A High-Density Genetic Linkage Map and QTL Fine Mapping for Body Weight in Crucian Carp (Carassius auratus) Using 2b-RAD Sequencing

A high-resolution genetic linkage map is essential for a wide range of genetics and genomics studies such as comparative genomics analysis and QTL fine mapping. Crucian carp (Carassius auratus) is widely distributed in Eurasia, and is an important aquaculture fish worldwide. In this study, a high-density genetic linkage map was constructed for crucian carp using 2b-RAD technology. The consensus map contains 8487 SNP markers, assigning to 50 linkage groups (LGs) and spanning 3762.88 cM, with an average marker interval of 0.44 cM and genome coverage of 98.8%. The female map had 4410 SNPs, and spanned 3500.42 cM (0.79 cM/marker), while the male map had 4625 SNPs and spanned 3346.33 cM (0.72 cM/marker). The average recombination ratio of female to male was 2.13:1, and significant male-biased recombination suppressions were observed in LG47 and LG49. Comparative genomics analysis revealed a clear 2:1 syntenic relationship between crucian carp LGs and chromosomes of zebrafish and grass carp, and a 1:1 correspondence, but extensive chromosomal rearrangement, between crucian carp and common carp, providing evidence that crucian carp has experienced a fourth round of whole genome duplication (4R-WGD). Eight chromosome-wide QTL for body weight at 2 months after hatch were detected on five LGs, explaining 10.1–13.2% of the phenotypic variations. Potential candidate growth-related genes, such as an EGF-like domain and TGF-β, were identified within the QTL intervals. This high-density genetic map and QTL analysis supplies a basis for genome evolutionary studies in cyprinid fishes, genome assembly, and QTL fine mapping for complex traits in crucian carp.

(S. . Among these methods, 2b-RAD is a simple and flexible method invented by S.  that adopts type IIB restriction enzymes for genome-wide genotyping. Compared with other GBS methods, it allows screening of almost every restriction site in the genome, regulates genome coverage more flexibly, and has a simpler protocol for library preparation (S. . 2b-RAD has been used successfully for constructing high-density maps and QTL fine mapping in several aquaculture species (Jiao et al. 2014;Shi et al. 2014;Cui et al. 2015;Tian et al. 2015;Dou et al. 2016;Fu et al. 2016;Palaiokostas et al. 2016).
Crucian carp (Carassius auratus) is widely distributed in Eurasia and is one of the most important freshwater fish species for Chinese aquaculture. Previous karyotype studies demonstrated that crucian carp has a diploid chromosome number of 2n = 100, a number twice as much most of other cyprinid fishes (2n = 50 or 48) (Knytl et al. 2013;Ojima et al. 1979;Kobayasi et al. 1970). Recent studies have revealed diploid (2n = 100), triploid (3n = 150), and tetraploid (4n = 200) crucian carp in natural populations, and these three ploidy populations often coexist with each other in natural waters (Liu et al. 2001;Xiao et al. 2011;Gui and Zhou 2010). Thereby, it was believed that crucian carp may have experienced multiple successive rounds of chromosome doubling, which provided an excellent material for genome duplication and evolution studies Gui and Zhou 2010;Zhou et al. 2002). Moreover, as an important farmed fish, the global aquaculture production of crucian carp reached 2.91 million tons in 2015 (FAO, 2016). Growth rate is one of the most important traits for breeding programs in crucian carp (Zhou et al. 2000a;Hulata 1995). Several strains of triploid crucian carp have been developed in China via gynogenesis or interspecific hybridization between crucian carp and common carp, which has allowed great progress in breeding, and significantly increased the production of crucian carp (Gui and Zhou 2010;Zhou et al. 2000aZhou et al. ,b, 2001Liu et al. 2001Hulata 1995). However, the genetic basis and architecture for growth modulation in crucian carp are still poorly understood because few genetic and genomic resources are available . It is well known that growth is a typical quantitative trait controlled by multiple genes known as quantitative trait loci (QTL), and may be influenced by environmental factors (Lynch and Walsh 1998;Mackay 2001;Mackay et al. 2009). Traditional selective breeding methods have encountered some difficulties such as uncertainty, extensive workload, being time-consuming, and being slow to take effect. Hence, molecular breeding methods, such as marker assisted selection (MAS) and genomic selection, are needed to accelerate breeding process in fish (Yue 2014;Liu and Cordes 2004;Tong and Sun 2015).
The objectives of this study include: (i) construction of a high-density SNP-based linkage map using 2b-RAD technology in diploid crucian carp; (ii) comparative genomics analysis between crucian carp and zebrafish (D. rerio), grass carp (Ctenopharyngodon idellus), and common carp (Cyprinus carpio). (iii) QTL fine mapping for body weight, and identification of candidate genes that may involve early growth of crucian carp. This study will provide a framework for further studies on genome evolution, comparative genomics, and fine-scale QTL mapping for economic traits in crucian carp.

MATERIALS AND METHODS
Mapping family and DNA extraction A total of 102 adult crucian carp was collected from a wild population in the Zhangdu Lake (Wuhan, China) and used as candidate broodfish in the production of mapping families. Genetic distances among these fish were evaluated using 20 previously reported microsatellite loci (Zheng et al. 2010). In April 2015, 17 candidate F1 full-sib families were established by crossing 14 sires and 17 dams via artificial fertilization. Larval fish of each family were raised in a plastic tank, and fed with brine shrimp (nauplii of the Artemia) twice a day. Body weights were measured and recorded for all families at the age of 2 months posthatch. One family with the highest genetic distance and largest within-family phenotypic variations in body weight was selected as the mapping panel for linkage map construction and QTL analysis for early growth in this study. A total of 160 progenies was randomly selected from this family, and a small piece of clipped fin was sampled from both parents and progenies and stored in 95% ethanol for DNA extraction. Genomic DNA was extracted from preserved fin tissues following a standard phenol-chloroform DNA extraction procedure (Sambrook and Russell 2001). The quality of DNA was checked by a NanoDrop 2000 spectrophotometer (Thermo Scientific), and 1% agarose gel electrophoresis. The concentrations of all DNA samples were adjusted to 50 ng/ml. All experimental animal programs involved in this study were approved by the Animal Care and Use Committee at the Institute of Hydrobiology, Chinese Academy of Sciences.
2b-RAD library construction and de novo genotyping Before library preparation, the number of possible restriction sites was calculated based on crucian carp genome draft assembly (unpublished data). The 2b-RAD library was prepared following the protocol originally described by S.  with minor modifications . Two parents and 160 offspring were used for the construction of 2b-RAD libraries. Briefly, 250 ng of genomic DNA was digested with 1 unit BcgI restriction enzyme (New England Biolabs) at 37°for 4 hr. The digested DNA was ligated at 16°for 8 hr with a 25 ml total volume reaction consisting of 1 unit T4 DNA ligase (New England Biolabs), 0.5 mM adapter 1 and adapter 2, 0.5 mM ATP (New England Biolabs). The ligation fragments were amplified in a 25 ml total volume consisting of 5 ml of ligated DNA, 0.5 mM P5 and P7 primer, 0.5 mM P4 and P6-BC primer, 0.75 mM dNTP (Shanghai Sangon, China), 5 ml 5· Phusion HF buffer, and 0.5 unit Phusion High-Fidelity DNA Polymerase (Thermo Scientific). The cycling conditions were: 98°f or 30 sec; 98°for 20 sec, 63°for 50 sec, 72°for 30 sec for 15 cycles, 72°for 5 min. The amplification products were purified via retrieval from 8% polyacrylamide gels. All libraries were pooled with equal amount of products from each library to make a final library which was sequenced on a lane at the HiSeq2500 platform (Illumina). The raw read data were archived at the NCBI Sequence Read Archive (SRA) under accession number PRJNA327320.
Raw reads were first trimmed to remove adapter sequences, and the terminal 2-bp positions. Reads without restriction sites or containing long homopolymers (.10 bp), ambiguous bases (N), low-quality sequences (more than five positions with a quality ,20) or mitochondrial origins were removed. The remaining trimmed reads with 32 bp in length were used for subsequent analysis. Filtered reads were analyzed with the software RADtyping program v1.0 (Fu et al. 2013) for de novo 2b-RAD genotyping. This software used stringent criteria in filtering candidate markers, and only those loci with at least four reads supporting were kept in the following analysis.

Linkage map construction
Only those markers that segregated in parents and could be genotyped in at least 80% of the offspring were used for further analysis. Markers that present significant segregation distortion in the x 2 goodness-of-fit tests (P , 0.05) were eliminated in the linkage analysis. Slightly distorted markers (P . 0.05) were also used for linkage analysis. A consensus linkage map was constructed by JoinMap 4.1 program (Van Ooijen 2006) with a threshold LOD score of 15.0. Male-and female-specific linkage map calculations were performed using the function of "Create Maternal and Paternal Population Nodes" in JoinMap 4.1, with a threshold LOD score of 10.0. The visualized linkage maps were drawn using MapChart v2.2 (Voorrips 2002). The linkage groups (LGs) of crucian carp were named according to their homologous groups Figure 1 Genetic distances and marker distribution of 50 linkage groups in the consensus linkage map of crucian carp. Within each linkage group, blue, green, and red lines represent maternal heterozygous SNPs, paternal heterozygous SNPs, and SNPs heterozygous in both parents, respectively.
n of common carp and zebrafish based on the results of comparative genome analyses. The expected genetic map length was calculated in two ways: G e1 (Fishman et al. 2001) and G e2 (Chakravarti et al. 1991), and the average of these two indexes was used as the predicted total genetic map length (G e ). The recombination ratio of female to male was calculated by the ratios of mean distances of each LG in the female and male maps.
Comparative genome analysis All mapped 2b-RAD sequences (32 bp) were first aligned against the crucian carp genome draft assembly (unpublished data). Those markers mapped at a single genome position were then extended by adding 100 nucleotide sequences from each side. A total of 5734 extended 2b-RAD sequences was searched against the genomes of zebrafish (Danio rerio, GRCz10), grass carp (Ctenopharyngodon idella), and common carp (C. carpio) using the basic local alignment search tool (BLASTN) with an e-value cutoff of 1e210. If a single marker sequence aligned multiple targets at different positions, only the top hit (lowest e-value) alignment was retained. The genomic synteny was visualized using the software Circos v0.67 (Krzywinski et al. 2009).
QTL analysis for body weight QTL mapping analysis was performed using the MapQTL 6.0 (Van Ooijen and Kyazma 2009) software program with a Multiple QTL Mapping (MQM) model. A mapping step size of 1 cM and five neighboring markers were used in QTL analysis. The genome-wide LOD threshold (significance level) or group-wide LOD threshold (suggestive level) were estimated using the permutation test (10,000 replicates) in MapQTL 6.0 with a confidence interval of 95%.
Identification of potential candidate genes Potential candidate genes within each QTL region were hunted through comparative genomics. We performed sequence similarity searches (BLASTN) for all QTL-associated SNP markers against the whole genome sequences of crucian carp (unpublished data) and common carp (Xu et al. 2014). Only annotated genes closest to the peak of corresponding QTL region were regarded as candidate genes.

Data availability
The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

RESULTS
2b-RAD genotyping A total of 6.56, 3.98 and 177.44 million reads were produced by 2b-RAD sequencing for the female parents, and male parents and 160 progenies (1.11 million reads per progeny), respectively. High-density linkage map All polymorphic SNP markers were grouped into 50 linkage groups, which is consistent with the haploid chromosome number of crucian carp (2n = 100) (Knytl et al. 2013). The consensus map consisted of 8487 markers and spanned 3762.88 cM, with an average interval of 0.44 cM (Figure 1 and Table 1). The number of markers per group distortion markers (7.15%) were distributed across the linkage maps, and they mainly gathered in some LGs such as LG1, LG6, LG12, LG22, and LG50 (Table 1 and Supplemental Material, Table S1). The expected map length was calculated to be 3807.22 cM (G e1 ) and 3809.784 cM (G e2 ), and the average G e of 3808.5 cM from these two estimation methods is taken as the expected genome length of crucian carp. The genome coverage of this consensus genetic map is then 98.8%. The female map comprised 4410 SNPs and spanned 3500.42 cM, with an average marker interval of 0.79 cM (Table 1 and Figure S1). Marker orders are largely conserved between female map and male map, although some LGs showed minor intrachromosomal rearrangements ( Figure S3). The average recombination ratio of femaleto-male is 2.13:1 (Table 1). Of the 50 LGs, 43 showed higher recombination rate in males than females, and, in contrast, seven LGs showed lower recombination rate in males than females. Interestingly, significant differences in recombination ratios between the female and male maps were observed in LG47 (13.76:1) and LG49 (9.17:1) (Table 1 and Figure S3).

Comparative genome analysis Comparative genomic analysis was performed between crucian carp
LGs and zebrafish, grass carp, and common carp chromosomes. For 5734 extended sequences of SNP markers, a total of 1072 markers was uniquely aligned to the genome of zebrafish. The results showed that two LGs of crucian carp were homologous to one particular chromosome of zebrafish, suggesting a clear 2:1 relationship of crucian carp LGs and zebrafish chromosomes (Figure 2A and Figure 3). Of the 1072 orthologous pairs, 966 pairs (90.1%) were located into the conserved syntenic blocks, revealing highly conserved synteny between crucian carp and zebrafish genome ( Figure 2B). Of note, 49% regions (37-75 Mb) on chromosome 4 of zebrafish had no homologous segments in crucian carp (Figure 2, A and B). A total of 1441 orthologous pairs between crucian carp and grass carp were identified, with 1267 (87.9%) mapped on those paired collinear blocks ( Figure 2C). Comparison analysis showed a 2:1 synteny between crucian carp and grass carp LGs, except for LG 24 in grass carp, which was aligned to four LGs in crucian carp (LG19, LG20, LG43, and LG44) ( Figure 2D and Figure 4). The comparative genomics results show that intensive chromosomal rearrangements were present between crucian carp and common carp ( Figure 2E and Figure 5). Of the 2201 markers uniquely anchored to the chromosomes of common carp, only 574 (26.1%) 1:1 orthologous pairs were identified ( Figure 2F). QTL for body weight QTL fine mapping based on above high-density genetic linkage map showed that eight chromosome-wide QTL associated with body weight were identified. These QTL distribute on five LGs (LG23, LG33, LG39, LG46, and LG49), with LOD scores from 3.71 to 4.91, and the phenotypic variance explained (PVE) ranging from 10.1 to 13.2% (Figure 6, Figure 7, and Table 2). The highest PVE value was detected for qBW46-a on LG46 with the LG positions of 32.579-32.808 cM ( Table 2). All QTL intervals in this study are ,1 cM except qBW23-a on LG23 (3.765 cM), with an overall average 0.87 cM for eight QTL intervals. Each QTL interval harbors 1-3 SNP markers (Table 2).
Potential candidate genes BLASTN searches of QTL-associated SNP sequences against the genome sequences of crucian carp and common carp gave five potential candidate genes ( Table 2). Markers in three QTL (qBW23-c, qBW39-a, and qBW49-a) failed to reveal any potential candidate genes because they were unable to align to any of the two genomes. Among those annotated genes with known biological functions (Table 2), epidermal growth factor-like domain (EGF-like domain) and transforming growth factor b (TGF-b) are two representatives.

DISCUSSION
With the development of high-throughput sequencing technology and the advancement of GBS methods, constructing well-defined genetic linkage maps using thousands of SNP markers is now available in many nonmodel organisms, including aquaculture species. High-density SNP linkage maps have been constructed using 2b-RAD sequencing in several aquaculture animals, including Chlamys farreri (3806 SNPs) (Jiao et al. 2014), pearl oyster (3117 SNPs) (Shi et al. 2014), Chinese mitten crab (10,358 SNPs) (Cui et al. 2015), sea cucumber (7839 SNPs) (Tian et al. 2015), bighead carp (3121 SNPs) , and gilthead sea bream (Palaiokostas et al. 2016) (12,085 SNPs). The number of markers on a linkage map is usually determined by the choice of restriction enzymes, the number of restriction enzymes, the total number of enzyme cut sites, and the polymorphism rate across the genome (Andrews et al. 2016). In this work, taking advantage of 2b-RAD technology, we constructed a high-density linkage map containing 8487 SNP markers with a resolution of 0.44 cM for crucian carp. This robust genetic linkage map will contribute to a better understanding of the genome structure, function, and evolution of crucian carp. In addition, the linkage map will be highly valuable to fine mapping for more complex traits and chromosome assembly of WGS in crucian carp.
Sex differences in recombination rates have been reported in a number of teleosts (Zhu et al. 2014). In many fishes, females usually have a higher recombination frequency than males . Recombination ratio differences between female and male have been observed in zebrafish (2.74:1) (Singer et al. 2002), rainbow trout (1.68:1 and 3.25:1) (Rexroad et al. 2008;Sakamoto et al. 2000), Atlantic halibut (1.89-2.53:1) (Reid et al. 2007), turbot (1.6:1) (Bouza et al. 2007), barramundi (2.06:1) (Wang et al. 2007), and grass carp (2.03:1) (Xia et al. 2010). In this study, the average female to male recombination ratio of crucian carp was 2.13:1, which was similar to many other fish species. In the majority of female LGs, the recombination ratios were obviously higher than those of male LGs. However, significant malebiased recombination suppression was observed in LG47 and LG49. This is a very interesting phenomenon worthy of attention. Previous studies have suggested many factors, such as sex chromosomes, regions around centromeres and/or telomeres, large areas of repetitive DNA, and heterochromatin could influence recombination rate (Ninwichian et al. 2012;Jones et al. 2013). Recombination suppression between sex chromosomes is a common phenomenon in vertebrates, and is important in maintaining the stability of the sex-determining regions and leads to the degeneration of Y or W chromosomes (Charlesworth et al. 2005;Bergero and Charlesworth 2009). In three spine sticklebacks and medaka, recombination in the sex determination region is reduced in the male linkage map relative to the female linkage map (Peichel et al. 2004;Roesti et al. 2013;Matsuda et al. 1999;Naruse et al. 2000). It was reported that crucian carp has an XX/XY sex chromosome system (Yamamoto 1974), and has a pair of heteromorphic chromosomes that were taken as X and Y chromosomes (Ruiguang 1982). In our study, significant male-biased recombination suppression on LG47 and LG49 may suggest that these two LGs are potential sex chromosomes of the crucian carp genome. Further studies are required to confirm this hypothesis, and elucidate the genetic mechanism for sex determination in crucian carp and other cyprinid fishes.
A high-density genetic map constructed with sequence-based markers makes it possible for genome evolution studies in nonmodel species (Z. . In this study, a detailed syntenic relationship was established between crucian carp LGs and zebrafish, grass carp, and common carp genomes via genetic maps and assembled genomes. In zebrafish, all chromosomal regions were covered by homologous loci of crucian carp except for half of Chromosome 4. The annotation information of zebrafish reference genome sequence indicated that this unique region harbored high gene duplication, high density of small nuclear RNAs (snRNAs), and may be related to sex chromosomes (Howe et al. 2013). This unique region presented obscure synteny with human, mouse, chicken, and common carp genomes (Howe et al. 2013;Xu et al. 2014), and our results were in accordance with those reports. The results of comparative mapping with zebrafish genome demonstrated a high level integrity of our linkage map. The conserved synteny with a clear 2:1 relationship between crucian carp LGs and zebrafish chromosomes was similar to that observed between common carp and goldfish (Peng et al. 2016;Xu et al. 2014;Zhang et al. 2013;Zhao et al. 2013;Kuang et al. 2016). In addition, a number of minor chromosome rearrangements were detected between crucian carp and zebrafish, which was similar to that reported in common carp (Xu et al. 2014). These findings suggested that crucian carp, similar to common carp, had undergone the 4R-WGD. Similarly, a clear 2:1 relationship was also observed between LGs of crucian carp and grass carp, except for LG24 of grass carp, which showed synteny to LG19, LG20, LG43, and LG44 of crucian carp. Previous studies indicated that LG24 of grass carp was orthologous to chromosome 10 and chromosome 22 of zebrafish (Z. , so, after the 4R-WGD, grass carp LG24 is syntenic to four LGs of crucian carp. Comparative analysis between crucian carp LGs, and common carp chromosomes demonstrated a 1:1 relationship, but with extensive chromosomal rearrangements. Crucian carp and common carp are two closely related species in the Cyprinidae, as they have the same number of chromosomes (2n = 100) and can produce hybrid offspring (Hulata 1995;S.J. Liu et al. 2001S.J. Liu et al. , 2016. Previous studies indicated that the common ancestor of crucian carp and common carp experienced the 4R-WGD 10.9-13.2 MYA, and the speciation occurred ranging from 8.1 to 12.9 MYA (Yuan et al. 2010;David et al. 2003;Chistiakov and Voronova 2009). A recent study confirmed that the 4R-WGD happened 8.2 MYA in common carp (Xu et al. 2014). Therefore, it was believed that crucian carp and common carp diverged Figure 7 The distribution of eight QTL for body weight on five genetic linkage groups of crucian carp.
n from each other after the 4R-WGD event (Somamoto et al. 2005;Chistiakov and Voronova 2009). Then, the duplicated genes faced different destinations, such as nonfunctionalization, subfunctionalization, neofunctionalization, and gene dosage effects, which is very important for biological evolution, adaptation, speciation, and diversification (Glasauer and Neuhauss 2014;Lien et al. 2016). In this case, it can be speculated that chromosome structural difference between crucian carp and common carp could have occurred since they evolved in different directions. A similar genome structure has been seen in salmonids, in which common ancestor undergone the 4R-WGD 80 MYA, and then the genome experienced extensive chromosomal rearrangements (Berthelot et al. 2014;Lien et al. 2016;Brieuc et al. 2014;Guyomard et al. 2012;Gharbi et al. 2006;Naish et al. 2013;Timusk et al. 2011). Furthermore, the rediploidization process followed by the 4R-WGD would also resulted in significant genome rearrangements (Lien et al. 2016). Therefore, large genomic reorganizations between crucian carp and common carp may be due to the independent genome evolution and rediploidization process after genome duplication. Growth is an important trait of interest in aquaculture species, and great efforts have been devoted to promoting the growth rate of crucian carp (Liu et al. 2001;Gui and Zhou 2010;Zhou et al. 2000aZhou et al. ,b, 2001). However, for genetic improvement of quantitative traits, traditional breeding methods have encountered some bottlenecks and problems (Mackay et al. 2009). QTL fine mapping and positional cloning of candidate genes may have been an efficient approach for breeding programs in aquaculture animals, especially for quantitative traits (Peng et al. 2016;Yu et al. 2015;Xiao et al. 2015;Tian et al. 2015;Shao et al. 2015;Li et al. 2015;Shi et al. 2014;). In the present study, the high-density linkage map allowed us to perform QTL fine mapping for body weight, and discover potential candidate genes for early growth stage in crucian carp. As a result, eight chromosome-wide QTL associated with body weight were located in five different LGs, which reflected the complexity of this polygenic trait. However, it is worth noting that the average map length of eight QTL was 0.87 cM, and this relatively narrow genomic region would facilitate further validating and positional cloning of potential major genes for growth in crucian carp. Among five potential candidate genes identified by genomic synteny analysis, EGF-like domain and TGF-beta may be most promising because their biological functions are likely associated with early growth and development in vertebrates. For example, EGF regulates cell proliferation and growth in human, which play an important role in cellular organization and membrane repair when the tissue is destroyed (Engel 1989). TGF-b modulates cell proliferation, differentiation, apoptosis, and immune regulation (Sporn et al. 1986). Our QTL results and uncovered candidate growth genes would lay a foundation for genetic improvement for growth in crucian carp. Nevertheless, we must recognize that our QTL results could suffer from the limitation of the use of a single family, which may have higher power to detect family-specific or rare QTL, but could also trade off the detection power of common QTL shared across families. In the future, joint multiple population analysis will be needed to detect QTL shared among multiple families with a wider scope of inference. In addition, taking into account the Beavis effect, QTL effects tend to be overestimated as sample size is relatively small, and the genetic architecture of the character highly polygenic, as is the case with body weight.

Conclusion
In summary, a high-density linkage map of crucian carp was constructed by 2b-RAD method with 8487 SNPs (0.44 cM/marker). Comparative genomics among four cyprinid fishes (crucian carp, zebrafish, grass carp, and common carp) provides new insights into genome duplication and chromosomal rearrangements in crucian carp. Eight chromosome-wide QTL for body weight were detected at quite narrow regions of five LGs with the PVE values of 10.1-13.2%. A few potential candidate genes, e.g., EGF-like domain and TGF-b, whose biological functions are likely involved in genetic regulation of early growth, were identified from those QTL intervals. Our present study provides valuable genetic and genomic resources for further studies of comparative genomics, genome evolution, chromosome assembly, and QTL fine mapping in crucian carp.