Construction of Three High-Density Genetic Linkage Maps and Dynamic QTL Mapping of Growth Traits in Yellow River Carp (Cyprinus carpio haematopterus)

To provide the theoretical basis for researching growth, development, and molecular marker-assisted breeding of the economically important Yellow River carp (Cyprinus carpio haematopterus) using dynamic quantitative trait locus (QTL) mapping, we constructed three genetic linkage maps from 207 progeny using a new modified genotyping-by-sequencing method. The three maps contained 16,886, 16,548, and 7482 single nucleotide polymorphism markers, respectively, with an average interval of 0.36 cM, 0.45 cM, and 1.00 cM. We identified 148 QTLs related to four growth traits that were located on 25 chromosomes from three growth stages of Yellow River carp. A total of 32, 36, 43, and 37 QTLs were associated with body length, height, width, and weight, respectively. Among them, 47 QTLs were detected for only one growth trait in one stage, but all of the other QTLs were co-localized. Of the 14 main QTLs, 13 were located on chromosome 12, which suggests the presence of growth-related genes on this chromosome. We then detected 17 candidate genes within 50 K upstream and downstream of the 14 main QTLs. This is the first report of the dynamic QTL mapping of growth traits of Yellow River carp, and the results can be used in future studies of growth, development, and molecular-assisted breeding of this species.


Introduction
The construction of genetic linkage maps and quantitative trait locus (QTL) mapping are very important techniques for the study of fish genetics and breeding. Many genetic linkage maps have been constructed for common carp (Cyprinus carpio) since the first map was reported in 2000 [1]. In subsequent studies, markers such as simple sequence repeat (SSR) and single nucleotide polymorphism (SNP) were gradually applied to construct the genetic linkage map of common carp [2,3].
Yellow River carp (Cyprinus carpio haematopterus) is one of the most economically important aquatic fish in the world. It is a well-known endemic species of common carp in the Yellow River valley, and it is also an important freshwater cultured fish in China. This species is famous in China because of its golden scales and red tail, long body shape, and delicious tender meat. It also has strong cold resistance, a high feed conversion rate, and is easy to catch. Although there are some QTL mapping [4][5][6] and adaptive evolution [7] studies of Yellow River carp, there are only a few genetic linkage maps, and dynamic QTL mapping has not been reported.
Some scholars have conducted QTL mapping of the growth-related traits of common carp. For example, Laghari et al. selected 109 SSR markers, 31 expressed sequence tag SSR (EST-SSR) markers, and 167 SNP markers, and they identified 7 QTLs associated with body weight and growth rate [8]. With the completion of common carp genome sequencing [9] and the development of simplified genome sequencing technology, the density of the common carp genetic linkage map and the accuracy of QTL mapping are increasing. Peng et al. selected 28,194 SNP markers and detected 22 growth-related QTLs and 7 sex-related QTLs [5]. Chen et al. conducted QTL mapping and a genome-wide association analysis of the head type traits of Yellow River carp and identified 12 loci that were significantly related to head type [6]. Feng et al. mapped the growth traits of the F2 population of Yangtze River carp and found 21 growth-related QTLs that could explain the observed phenotypic variation (PVE) [10]. Wang et al. constructed a high-density genetic linkage map based on a common carp F2 family using a streamlined restriction site-associated DNA genotyping method based on sequencing the uniform fragments produced by type IIB restriction enzymes [4]. They then carried out QTL mapping studies and identified 12 chromosomal QTLs and 2 genome-wide growth-related QTLs.
Many QTL mapping studies of fish have focused on a certain stage of development, which is generally defined by the researcher based on the study question. For example, in aquaculture studies, QTL mapping of fish weight could be focused on the time when a fish grows to commercial size. This approach allows only a static study of the cumulative effect of these traits, and it cannot reflect the effect of genes on the whole growth process. To address this issue, Zhu proposed the conditional QTL and unconditional QTL analysis methods based on the composite interval mapping of a mixed linear model, which some scholars refer to as dynamic mapping [11]. Wu et al. developed software to locate the net genetic effect QTL from time "1" to "t" [12]. Other researchers have conducted dynamic QTL analyses of several species, including the potato (Solanum tuberosum) [13], wheat (Triticum aestivum) [14], rice (Oryza sativa) [15], and cotton (Gossypium hirsutum) [16].
Similar methods have been used for dynamic QTL mapping in fish. McClelland and Naish mapped 53 QTLs related to three traits by measuring the growth rate, body length, and body weight of coho salmon (Oncorhynchus kisutch) over eight time periods [17]. Laghari et al. studied the growth rate and body weight of 92 individuals of F1 common carp at 10, 11, and 12 months of age and identified 7 QTLs related to body weight [8]. Gutierrez et al. marked Atlantic salmon (Salmo salar) fry from five families with passive integrated transponders (PIT) when they reached 25 g, and then used SNP markers to locate the weight-related QTLs at four time points [18]. In the weight-related QTL detection of gilthead seabream (Sparus aurata) at four time points in a growth cycle, the QTLs of linkage group 1 were significant at all time points, but the effect varied with time, and one QTL of linkage group 21 seemed to affect the later growth weight of the species [19].
In 2018, Qi et al. [20] invented a modified genotyping-by-sequencing (GBS) method, which is a simplified genome sequencing technique based on GBS technology that is better at constructing polyploid linkage maps than previously available methods [21]. In the modified GBS, the size of the recovered fragment is adjusted by adjusting the volume ratio of the magnetic bead solution to the connecting product, which avoids uncontrollable processes, such as glue cutting and recovery, and improves the repeatability of the technology. Qi et al. also used a new analysis process to greatly improve genotyping accuracy [21]. The modified GBS is suitable for outcrossed polyploid species, complex genome species, and species with high repetitive sequences [20]. In this study, we will use the modified GBS and the new analysis process to analyze the experimental results. By comparing the QTLs associated with growth traits at different time periods, we may provide a theoretical basis for studying the growth and development mechanisms at work in Yellow River carp. Our results may also be applicable to the molecular marker-assisted breeding of Yellow River carp.

Mapping Family and Measuring Phenotype
In May 2017, one female fish was selected from a fast-growing Yellow River carp family established in 2015 and one male fish was selected from a slow-growing Yellow River carp family. These parent fish were allowed to mate to establish the F1 family. When the average body length of the progeny reached about 10 cm at 5 months old, 300 fish were randomly selected and injected with PIT markers. Their body length, height, thickness, and weight were measured at 13, and 17 months of age using an electronic balance (Wuxin Weighing Instrument Co., Ltd., Ningbo, China) and Vernier calipers (Ningbo Deli Co., Ltd., Ningbo, China). During the last measurement, a piece of the caudal fin of each fish was cut and preserved with anhydrous ethanol. After dead fish and those whose markers could not be detected were removed from the sample, 207 F1 progeny remained as the mapping population. Their genomic DNA was isolated from fins using the traditional proteinase-K digestion and phenol-chloroform extraction method. The quality of the genomic DNA was checked using a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Waltham, MA, USA). The average value, standard deviation of each trait and Pearson correlation coefficient between traits were performed using SPSS24.0 (IBM Corp, Armonk, NY, USA).

GBS Sequencing and Creation of Maternal (HA), Paternal (AH), and Both Parent (HH) SNP Datasets
The DNA samples were sent to Shanghai Ouyi Biotechnology Co., Ltd. (Shanghai, China) for sequencing. The libraries of the parents and 207 F1 progeny were prepared following the methods of Qi et al. [21]. Read processing, SNP calling and filtering were conducted following the methods of Qi et al. [21]. Chi-square tests were used to determine the markers' deviation from the expected Mendelian segregation ratios. Highly distorted markers (p-value ≤ 1 × 10 −10 ) were discarded. The remaining SNPs were categorized according to their segregation ratio and parental genotype. Population type was set as "backcross" (BC) for the maternal and paternal data sets, and as "F2 intercross" for the HH data set [21]. SNPs that were heterozygous in both parents with progeny segregating as A:H:B (1:2:1) constituted the HH dataset. SNPs with genotypic score A or B in parent 1 and H in parent 2 and segregating as A:H or B:H (1:1) were combined, and the B scores were converted to A scores to form the paternal (AH) dataset. Similarly, the H:A and H:B markers were combined, and the B scores were converted to A scores to yield the maternal (HA) dataset, following the methods of Qi et al. [21]. The GBS reads were submitted to NCBI-SRA (Acc. PRJNA788161).

Map Construction and Estimation of Genome Size
The HA, AH, and HH genetic maps were constructed using a combination of MSTMap [22] and MAPMAKER (modified from Lander et al. [23] as described by Qi et al. [20]. The population type was set as "backcross" for the maternal and paternal datasets and as "F2 intercross" for the HH dataset. Co-segregating SNPs were added to the framework map to the same location as their representative marker using an in-house Python script. Linkage maps were drawn with MapChart [24].

QTL Location and Consensus QTL Analysis
The body length, height, thickness, and weight measured for the F1 progeny were used as inputs for the QTL analysis. QTL mapping was conducted with WinQTLCart 2.5 [25] using the composite interval mapping model. The population type for QTL analysis was set as "SF2" for the HH maps and as "B2" for the maternal and paternal maps. A walk speed of 0.5 cM was used for QTL identification. The logarithm of odds (LOD) threshold for significant QTLs (p ≤ 0.05) was determined by 1000 permutations. After the QTLs were located, all QTLs located within 1M base pairs of each other on the physical map were identified as being consensus QTLs.

Screening for Growth-Related Candidate Genes
Potential growth-related candidate genes within QTLs with LOD value > 5.5 and PVE > 10% were identified. Briefly, we screened the +50 kb genome regions surrounding the significant SNPs based on the common carp reference genome (BIGD Genome Warehouse, https://bigd.big.ac.cn/gwh/Assembly/497/show (accessed on 9 May 2019) and annotated the candidate genes by conducting BLAST analysis against the Swiss-Prot database.

Phenotype Data
It can be seen from Table 1 that the body weight in the T1T2 period (4 months) accounted for more than 66% of the body weight in the T2 period (17 months). Pearson correlation analysis showed that there is a significant correlation between each trait at the 0.01 level (Table 2). However, the correlation coefficients between the traits differed greatly, as the smallest was only 0.348, and the largest was 0.987. The correlation coefficients of BL, BH and BT in the T1T2 stage for all traits were <0.8, and the correlation coefficients among all other traits were >0.8. These results show that the increase of body mass in the T1T2 stage is caused by many traits, rather than being caused mainly by a certain trait.

GBS Tags
High-throughput sequencing of the 209 GBS libraries from the 2 parents and 207 progeny yielded 6.9 and 7.2 million clean reads for the dam and sire, respectively, and 1.1 billion reads for the offspring, with an average of 5.6 million reads per progeny. The sequencing values were converted into genotypes, and only the recombination events across the progeny were retained for linkage map construction. The simplified linkage map (one marker per 2 cM) was used for the permutation calculation. Finally, 16,886 SNPs that were heterozygous only in the female parent were used to construct the maternal map, 16,548 SNPs that were heterozygous only in the male parent were used to construct the paternal map, and 7482 SNPs that were heterozygous in both the male and female parents were used to construct the HH map.

Construction of the Genetic Linkage Map
Three genetic linkage maps (maternal, paternal, and HH) that included 50 linkage groups were constructed ( Figure 1). The maternal map contained 16,886 SNPs and spanned 6103.6 cM, with an average marker interval of 0.36 cM (Table S1, (Table S1).

QTLs for Growth Traits
QTL fine mapping based on the high-density genetic linkage maps identified 148 QTLs, 8 of which were not located on chromosomes (Table S2). Thirty-two QTLs associated with body length were distributed in 11 LGs and 8 chromosomes ( Figure S1), with logarithm of odds (LOD) scores ranging from 3.01 to 13.20 and PVEs ranging from 0.10% to 13.20%. The 36 QTLs associated with body height were distributed in 14 LGs and 11 chromosomes ( Figure S2), with LOD scores ranging from 3.06 to 7.91 and PVEs ranging from 0.38% to 12.43%. Forty-three QTLs associated with body thickness were distributed in 18 LGs and 17 chromosomes ( Figure S3), with LOD scores ranging from 3.00 to 7.60 and PVEs ranging from 0.01% to 12.20%. The 37 QTLs associated with body weight were distributed in 12 LGs and 10 chromosomes ( Figure S4), with LOD scores ranging from 3.01 to 6.97 and PVEs ranging from 2.26% to 13.06%. Table S2 and Figure 2 show that there were great differences in the number and positions of QTLs for different traits in different stages, and some QTLs were expressed only during a specific time period.

Candidate Growth-Related Genes
After BLAST searching against the genome of Yellow River carp using extended SNP marker sequences, 17 potential candidate genes were identified based on the SNP position in the genome ( Table 4). The UPF0609 protein C4orf27 homolog (cd027), ankyrin-2 (ank2), and calcium/calmodulin-dependent protein kinase type II delta (camk2d) genes located on chromosome 1 were associated with body height at stage T1T2, and the other 14 candidate genes were located on chromosome 12. Synaptotagmin-1 (syt1) was associated with five QTLs covering BW, BT, BH and three stages. Zinc finger NFX1-type containing 1 (znfx1) was associated with body height at stage T1. Small nuclear ribonucleoprotein-associated protein B (snrpb) was associated with six QTLs covering all four traits and three stages. Itchy E3 ubiquitin protein ligase (itch) was associated with five QTLs covering BH, BT, BW and three stages. N-terminal EF-hand calcium-binding protein 1 (neca1) and protein CBFA2T2 (cbfa2t2) were associated with five QTLs covering all four traits and three stages. The remaining eight candidate genes were associated with five QTLs covering BH, BT, and BW at different times. Most of these candidate genes were related to cell proliferation, energy metabolism, immunity, and growth.
The modified GBS method is a simplified genome sequencing technique based on GBS technology that is better at constructing polyploid linkage maps than previously available methods [20,21]. Using this technique, we constructed three independent genetic linkage maps of Yellow River carp. Prior to our study (Table S3), the largest sample size of carp used to create a linkage map was 190 (627 markers) [37], and the highest density was 0.38 cM [5]. Thus, the mapping population used in this study (207 progeny) is the largest to date, and the map density is the highest to date. We constructed three independent genetic linkage maps because (i) we used the F1 population, (ii) the construction methods for the different markers differed, and (iii) the three maps were more in line with the genetic characteristics of the mapping population. The population type was set as "backcross" for the maternal and paternal datasets and as "F2 intercross" for the HH dataset. The mapping method used in this paper is different from the previous mapping methods [4,5], and the map made by this new method provides a new choice for the QTL location of Yellow River carp.
Dynamic QTL mapping can detect the expression of trait-related genes in different growth stages. We detected 148 QTLs located on chromosomes in different growth periods, and 47 of them were only detected for a certain trait in a certain time period. Seven QTLs were detected for different traits in the same time period, 37 QTLs were detected in two time periods, and 57 QTLs were detected in three time periods. These results show that some genes are continuously expressed throughout the growth and development period, whereas others are only expressed at a certain stage of development. This is consistent with the results of the correlation analysis of growth traits, which show that there is a certain correlation among the four traits in the three periods. Among them, the correlation between body length, body height, body width and the other traits in the T1T2 period was relatively low (0.3-0.8), while the correlation coefficients among other traits were >0.8. Gutierrez et al. reported similar results for Atlantic salmon. They analyzed the QTLs associated with net weight gain between time periods and found that the QTLs associated with weight appeared in different time periods in different families [18]. Similar results were found in the analysis of weight-related QTLs at four time points in one growth cycle of gilthead seabream. Some QTLs can be detected in multiple periods, while others can only be detected in a certain period [19]. We detected 47 QTLs in the T1T2 period of Yellow River carp growth, and 22 of them only existed in this growth stage. Most of the QTLs were different, so we then located conditional QTLs that were associated with growth traits. Compared with previous research about the growth-related QTL mapping of Yellow River carp only for a specific time [4,5], in this study we located the growth-related QTLs in three different time periods, which included not only two specific time points, but also the growth traits in a specific time period. This approach allowed us to explore the changes to the growth-related QTLs of the Yellow River carp.
Although conditional QTL mapping is rarely used in aquatic animals, it can play a very important role in breeding. For example, after injecting fish with PIT tags, their body mass can be measured before and after winter and then the QTLs associated with cold tolerance can be mapped [38]. The early growth traits of fish are greatly affected by motherhood [39,40], and dynamic QTL mapping can reduce the influence of the maternal effect and locate growth-related QTLs more exactly. Additionally, different types of fish have different growth characteristics, and measuring growth traits from the early growth stages to adulthood allows for more accurate QTL mapping. Conditional QTL can also be used to study gonadal development by detecting the condition before and after gonadal development in fish [41].
Many QTLs exhibit pleiotropy, meaning that a QTL is associated with different traits. Sometimes, one gene may affect multiple growth traits at the same time. For example, Cqtl-7, Cqtl-19, and Cqtl-21 in this study were detected for different traits in all three time periods. Similar results were reported previously for barley (Hordeum vulgare) [42] and Yellow River carp [4]. Through our QTL mapping analysis of the growth traits of Yellow River carp, we found that the QTL distribution showed a regionalization trend, with some growth traits having co-located QTLs. The QTLs with LOD > 5.5 and PVE > 10% were localized on chromosome 1 (BH-MM-T1T2-2) and chromosome 12 (the other 14 QTLs), which indicated that the major genes related to growth were localized on chromosomes 1 and 12. However, these results differed from those reported by Wang et al. [4] and Peng et al. [5], possibly due to the different genetic backgrounds of the families in the different studies. We also identified 47 QTLs associated with growth traits that were detected only in a single time period and for a single trait, which indicated that these QTLs were greatly affected by the growth period and environment. The LOD and PVE values of these QTLs were relatively small, which suggested that the QTLs detected in only one period and trait had little effect on growth traits.
In this study, the QTLs with LOD value > 5.5 and PVE > 10% were mainly concentrated on chromosome 12 and contained 14 candidate genes. Most of these genes are associated with cell proliferation, energy metabolism, immunity, and growth, but there may also be some as-yet unidentified functions. In future studies we will explore the relationship between alleles of these genes and growth, and we will screen for a number of molecular markers and alleles closely related to growth in order to provide a basis for the molecular marker-assisted selection of Yellow River carp.