Construction of a chromosome-level genome assembly for genome-wide identification of growth-related quantitative trait loci in Sinocyclocheilus grahami(Cypriniformes, Cyprinidae)

The Dianchi golden-line barbel (Sinocyclocheilus grahami,Cypriniformes, Cyprinidae) (Figure 1A) is a small-sized cyprinid with restricted distribution in the Dianchi basin and surrounding streams in Yunnan Province, China. Famed as a high-yield and nutrition-enriched fish, this species contains higher amount of protein, lipids, and essential amino acids than many common market fish, such asCtenopharyngodon idellus(Valenciennes, 1844) (grass carp) andHypophthalmichthys nobilis(Richardson, 1845) (bighead carp)(Zhao et al., 2013). However, theS. grahamipopulation has declined sharply since the 1960s due to habitat damage,water pollution, and alien species invasion (Yang et al., 2007).As such, it was listed as an animal of Second-Class National Protection in 1989, and as an endangered fish in the “China Red Book of Endangered Animals” in 1998 (Le & Chen, 1998).

Since 2007, our team has successfully achieved the artificial breeding ofS. grahami(Yang et al., 2007), which has not only helped in avoiding its wild extinction, but also opened a new era for its utilization. Moreover, after four generations of artificial selection, a new national breed (“S. grahami, Bayou No. 1”) with accelerated growth and weakened inter-muscular bones was certified in 2018. This new breed highlights a successful paradigm of wildlife protection in China, i.e., from endangered species conservation to sustainable utilization. It is foreseeable that with “S. grahami,Bayou No. 1”, focus will shift to the improvement of growth traits, which are among the most important economic characteristics in aquaculture (Yin et al., 2020).

Molecular marker-assisted selection (MAS) is an efficient method for trait improvement as it can effectively shorten selective breeding time and minimize possible negative risks to the environment (Yue, 2014). If high-precision chromosomal genomes are available, QTLs and candidate genes can be identified and applied for MAS. Commonly used methods for identifying QTLs and candidate genes in fish species studies include genome-wide association analysis (GWAS), QTL analysis, restriction site-associated DNA sequencing (RADseq), and transcriptome sequencing (RNA-seq) (Li et al.,2018; Liu et al., 2017; Salem et al., 2012; Yu et al., 2016).QTL-seq, which combines the advantages of both bulked segregant analysis (BSA) and whole-genome resequencing, is an innovative, rapid, and effective way in which to locate QTLs(Gu et al., 2018; Li et al., 2017; Takagi et al., 2013).

In the current study, we constructed a high-density linkage map ofS. grahamiusing RAD-seq based on a full-sib family,including 169 offspring and their parents (F1 generation from domesticated Muyanghe population) cultivated in the Endangered Fish Conservation Center (EFCC) of the Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences (CAS). Based on the generated high-density genetic linkage map, we constructed a chromosome-level genome assembly with the assistance of our previous scaffold-level genome sequences (Yang et al., 2016). We also used the obtained chromosome-level genome as a reference and performed QTL-seq for growth based on sibling (～2 000 individuals) and random (～20 000 individuals) populations ofS. grahami. All animal experiments were approved by the internal review board of KIZ, CAS (approval ID: SMKX-SQ-20170308-073). Detailed descriptions are provided in the Supplementary Materials and Methods.

Following RAD-seq and quality control, a total of 286 Gb of raw reads were generated through an Illumina HiSeqX-Ten platform. In total, 25 646 single nucleotide polymorphism(SNP) markers were obtained across all samples in the full-sib family, with 18 874 SNPs detected in more than 70% of the progeny. Among them, 12 412 SNPs with normal Mendelian segregation patterns were retained to construct a linkage map,with classification into 48 linkage groups.

With the assistance of our previous scaffold-levelS.grahamigenome (Yang et al., 2016), we assigned the linkage groups into 48 pseudo-chromosomes. Our final chromosomelevel assembly contained 12 412 high-quality SNPs, with a genome size of 1.49 Gb (see Supplementary Tables S1, S2),accounting for 84.99% of the scaffold-level assembly(1.75 Gb). Figure 1B shows the distributions of SNPs, genes,and GC content at 100 kb genomic intervals, as well as the one-to-one syntenic relationships of orthologous pairs. In total,20 558 gene pairs (430 syntenic blocks) betweenS. grahamiandDanio rerio(Hamilton, 1822) (zebrafish) were obtained after genome-wide alignments, thus supporting a 2:1 syntenic relationship between the two species for almost all chromosomes (Figure 1C).

Based on QTL-seq analysis (Takagi et al., 2013), two bulks from the sibling population (～2 000 individuals) were constructed, i.e., fastest-growth bulk consisting of 30 extremely large individuals and slowest-growth bulk consisting of 30 extremely small individuals. Two parallel bulks from the random population (～20 000 individuals) were also constructed, with 30 extreme-sized individuals contained in each bulk. In the four constructed bulks, the mean body lengths (and weights) of the fastest-and slowest-growth bulks were 47.83±2.71 mm (1.91±0.32 g) and 22.52±1.58 mm(0.2±0.05 g) in the sibling population, and 45.5±4.48 mm(1.88±0.79 g) and 20.83±2.11 mm (0.15±0.037 g) in the random population, respectively. In each population, body length and weight showed significant differences (P＜0.01)between the two extreme bulks (Supplementary Figure S1).

After whole-genome resequencing and quality control of the four bulk-pooled samples from the two populations, a total of 1 584 million reads (150 bp in length) were retained for subsequent analyses. According to the SNP identification procedure in QTL-seq analysis, a total of 282 476 and 281 611 sliding windows were detected in the sibling and random populations, respectively, with more than 15 SNPs in each window. Correspondingly, theFst values of each window ranged from 0.001 to 0.138 and 0.001 to 0.188, and the lowestFst values for the top 1/1 000 windows (Fst cutoff) were 0.085 and 0.098 in each population, respectively. After filtering those chromosomes that contained sliding windows (＜4) withFst values higher than the cutoff, as well as those regions with only sporadic sliding window distributions, we retained 12 QTLs from the sibling population and 29 QTLs from the random population (Figure 1D, E).

第一，思想准备。各级政府及其水行政主管部门要达成共识，真正转变观念，切实将思想认识统一到使市场在水资源微观配置中起决定性作用上来。同时，在全社会也要达成共识。

Additionally, when we combined Fisher’s Exact test and bulk frequency ratio (BFR) parameter (BFR≥4), 152 and 310 significant SNPs were identified in the sibling and random populations, respectively. The major SNPs (≥8) screened from the sibling population were distributed in Chr3, Chr17, Chr33,Chr34, and Chr46 (Figure 1F), while the major SNPs screened from the random population showed a relatively scattered distribution in Chr3, Chr10, Chr11, Chr13, Chr14, Chr15,Chr17, Chr18, Chr19, Chr24, Chr28, Chr30, Chr36, Chr37,Chr39, and Chr48 (Figure 1G).

To identify more reliable QTLs, we combined the significant SNPs identified from both the sliding windows and BFR parameter, and observed that growth-related QTLs converged in Chr3, Chr17, Chr33, Chr34, and Chr46 in the sibling population, and in Chr3, Chr11, Chr14, Chr15, Chr17, Chr18,Chr28, Chr36, Chr37, and Chr39 in the random population.Furthermore, when considering both sibling and random populations together, the overlapping QTL intervals between the two populations were in Chr3 (14.9-39.1 Mb) and Chr17(4.1-27.4 Mb) (Figure 1H, I), suggesting that these QTLs represent the greatest possible growth-related regions inS.grahami.

Simultaneously, we identified several genes within the two major candidate QTL intervals that contained significant SNPs(BFR≥4), which were thus considered as candidate genes relevant toS. grahamigrowth. These genes includedmitogenactivated protein kinase kinase 5(map2k5), which had two SNPs with BFR values of 8.4 and 5, respectively;signal transducer and activator of transcription 1-alpha/beta(stat1)andPHD finger protein 21A(phf21a), which contained SNPs with BFR values ≥5; andSRY-box 6(sox6),SMAD family member6(smad6),tetraspanin-3(tspan3), andunconventional myosin-Va(myo5a), which all contained SNPs with BFR values ≥4 but ＜5 (see Supplementary Figure S2).

The chromosome-level genome with 48 pseudochromosomes constructed here is consistent with previous karyotypic studies ofS. grahami(Li et al., 1983; Xiao et al.,2002). In addition, the general 2:1 syntenic relationship betweenS. grahamiandD. reriodemonstrated that the ancestor ofS. grahamiunderwent another whole-genome duplication (WGD) after the teleost-specific genome duplication (TSGD), with the TSGD maintaining only 24-25 chromosomes in most extant fish, such as the model species zebrafish (Xu et al., 2014).

Figure 1 Chromosome-level genome and growth-related QTLs of S. grahami

In MAS practice, common QTLs between multiple populations/families are usually treated as valuable targets for containing relatively stable gene sites over long-term evolution(Barton & Keightley, 2002). Therefore, the two major candidate QTL intervals in Chr3 (14.9-39.1 Mb) and Chr17(4.1-27.4 Mb) should be given priority in future MAS studies onS. grahami.

Moreover, according to the hypothesis that SNPs with greater BFR values locate at loci with more dominant relationships to target traits (Wang et al., 2013; Yao et al.,2017), we identified several candidate genes related to growth within these two QTL intervals. For example, the top candidate genemap2k5is highly conserved among various species. It is an upstream kinase ofextracellular signal regulated kinase 5(erk5) with abundant expression in skeletal muscle (Dinev et al., 2001), and themap2k5-erk5pathway is critical in muscle cell differentiation (Dinev et al., 2001), myogenesis (Chen et al., 2017; Mauro et al., 2002), neural differentiation (Nishimoto et al., 2005), and growth factor mediation (Carter et al., 2009).Therefore, it is reasonable to infer thatmap2k5mutations may have affected themap2k5-erk5pathway and induced differences in growth between the fastest- and slowest-growth bulks. Of course, other candidate genes, such asstat1,phf21a,smad6,sox6,myo5a, andtspan3, may also have played important roles in the differences in growth between bulks as they are involved in functions related to neuronal development, endochondral bone formation, chondrogenesis,skeletal muscle development, myeloid development, cell differentiation, and immunity (Dinh et al., 2010; Estrada et al.,2011; Jackson et al., 2015; Kim et al., 2012; Lefebvre, 2019;Najjar & Fagard, 2010; Rancura et al., 2008; Song et al.,2011; Yeh & Klesius, 2012).

Nevertheless, we also identified several QTL differences between the two populations. In theory, these QLTs should not be excluded as candidate QTLs, but further research is required. Lv et al. (2016) identified 165 QTLs for growthrelated traits in eight distinct families ofCyprinus carpio(Linnaeus, 1758) (common carp), and found only 36 QTLs were common in more than two families, with all others found to be distinct among families. Thus, both previous study and our present research suggest that differences between populations/families would, to a great extent, represent population/family-specific QTLs as growth is a complex trait affected by many genes. From this point of view, to improve MAS efficiency in complex traits such as growth, more populations/families with both simple and complex genetic backgrounds should be investigated.

DATA AVAILABILITY

All relevant sequences were deposited in the National Center for Biotechnology Information (NCBI) under accession No.PRJNA702560 and China National GeneBank (CNGB) under accession No. CNA0019203.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

J.X.Y., W.S.J., and Q.S. designed the study. Y.H.Y., X.A.W.,Y.W.Z., A.L.W., X.D.H., M.W., and X.F.P. collected the samples and performed the laboratory work. Y.H.Y., X.H.Z.,R.H.L., X.X.S., X.X.Y., and C.B. performed data analysis.Y.H.Y., X.H.Z., X.A.W., and R.H.L. drafted the paper. W.S.J.,J.X.Y., and Q.S. revised the paper. All authors read and approved the final version of the manuscript.

ACKNOWLEDGMENTS

We would like to thank all staff from the Endangered Fish Conservation Center (EFCC) of the Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences (CAS), for assisting us in the cultivation of fish materials.