Construction of a high-density genetic map and fine QTL mapping for growth and nutritional traits of Crassostrea gigas
BMC Genomics volume 19, Article number: 626 (2018)
Both growth and nutritional traits are important economic traits of Crassostrea gigas (C. gigas) in industry. But few work has been done to study the genetic architecture of nutritional traits of the oyster. In this study, we constructed a high-density genetic map of C. gigas to help assemble the genome sequence onto chromosomes, meanwhile explore the genetic basis for nutritional traits via quantitative trait loci (QTL) mapping.
The constructed genetic map contained 5024 evenly distributed markers, with an average marker interval of 0.68 cM, thus representing the densest genetic map produced for the oyster. According to the high collinearity between the consensus map and the oyster genome, 1574 scaffold (about 70%) of the genome sequence of C. gigas were successfully anchored to 10 linkage groups (LGs) of the consensus map. Using this high-qualified genetic map, we then conducted QTL analysis for growth and nutritional traits, the latter of which includes glycogen, amino acid (AA), and fatty acid (FA). Overall, 41 QTLs were detected for 17 traits. In addition, six candidate genes identified in the QTL interval showed significant correlation with the traits on transcriptional levels. These genes include growth-related genes AMY and BMP1, AA metabolism related genes PLSCR and GR, and FA metabolism regulation genes DYRK and ADAMTS.
Using the constructed high-qualified linkage map, this study not only assembled nearly 70% of the oyster genome sequence onto chromosomes, but also identified valuable markers and candidate genes for growth and nutritional traits, especially for AA and FA that undergone few studies before. These findings will facilitate genome assembly and molecular breeding of important economic traits in C. gigas.
The Pacific oyster, Crassostrea gigas (C. gigas) originated in Asia, and has become an important cultured mollusc in many countries as a result of its rapid growth rate and high yield [1, 2]. In 2016, the global aquaculture production of C. gigas was estimated to be 625,925 metric tons with a value of approximately 1.34 billion dollars (FAO, 2016). As an edible aquatic product, quality and flavour plays a key role in the commercial value of C. gigas. Growth and nutritional traits are not only associated with the oyster production but also important indicators of oyster quality. Oysters are rich in glycogen, essential amino acids (EAAs), unsaturated fatty acids (UFAs), and vitamins . High levels of nutrients give the oyster high physiological and medicinal value. For example, glycogen from the oyster meat can be adsorbed directly by the human body, thus lighten the burden of the pancreas . As a sulphur-containing neurotransmitter, taurine (Tau) plays a key role in the development of the nervous and cardiovascular systems [5, 6], and is involved in many physiological and metabolic processes [7, 8]. The protein content of the oyster dry meat weight was reported as high as 57%, and the quality and composition of its EAAs were found to be superior to the milk [3, 4]. In addition, the UFAs of the oyster meat, such as C20:5ω3 (EPA) and C22:6ω3 (DHA), were reported to prevent cardiovascular and cerebrovascular diseases [9, 10]. Because of its commercial and industrial importance, breeding work has been performed to improve the quality of C. gigas. Compared with the growth traits, which have got great attention in previous phenotype-dependent selective breeding programs [11,12,13,14], few nutritional traits of C. gigas, other than glycogen , have been subjected to molecular breeding. Though nutritional traits of the oyster were reported to vary with the environment, like the season and sea area [16, 17], the oyster family cultured in the same environment exhibited great phenotype variation of the nutritional traits (Additional file 2, ), which indicated obvious genetic effect for nutritional traits.
Linkage analysis is an important tool for quantitative trait loci (QTL) mapping, marker-assisted selection (MAS) programs, and genome assembly. For the oyster, a total of 20 articles have reported studies concerning genetic map construction or QTL mapping [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38]. To date, QTLs for growth [31, 36], glycogen [34, 37], disease resistance [24, 29], and shell color [34, 38] traits have been detected in previous studies. However, no QTLs for AA and FA traits have been detected before and reported maps did little for the assembly of the oyster genome because of the absence of the genome or enough markers. In this study, we selected the double-digest genotyping by sequencing (ddGBS) method to genotype our full-sib mapping family. With thousands of polymorphic loci detected, we constructed a high-density genetic map. We not only evaluated the collinearity between the genetic map and the genome of C. gigas, but also conducted QTL mapping for many economic traits, including AAs and FAs that had never undergone QTL mapping or association analysis in this species. This work assists to assemble the genome sequence of C. gigas onto chromosomes, meanwhile help illuminate the genetic mechanism of important economic traits so as to promote MAS programs of C. gigas.
Mapping population, DNA and RNA extraction, and cDNA cloning
In July 2013, the mapping family was constructed by mating one maternal oyster (Changli, China) with another paternal oyster (Qingdao, China). All oysters were pasted with different labels on the shell to distinguish from natural spat before culturing in the sea in September 2013. In January 2015, a total of 169 F1 offsprings were randomly selected from this family as the mapping population. Mantles of the oysters were stored at − 80 °C for DNA and RNA extraction and the other tissues were sampled and stored at − 20 °C for nutrients detection. DNA was extracted with a DNeasy® Blood & Tissue Kit (QIAGEN, Germany) according to the user’s instructions. RNA was extracted using TRIzol reagent (Invitrogen, USA) according to manufacturer’s protocol. cDNA was reverse-transcribed from 1 μg of total RNA in a 20-ul reaction mixture using PrimeScript RT reagent kit with gDNA Eraser (TaKaRa, Japan) following the manufacturer’s instructions. All experiments were conducted with approval from the Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China.
Phenotype detection and analysis
Growth and nutritional traits were analysed for all 169 full-sib offsprings. Growth traits include shell height (SH), shell length (SL), shell width (SW), body weight (BW) and weight of soft tissue (STW), while nutritional traits include glycogen, 18 AAs, and 22 FAs. All nutritional traits were detected using the dry soft tissue that was dried with the freeze dryer (Boyikang, China). Glycogen content was detected using a liver/muscle glycogen assay kit (Nanjingjiancheng, China) according to the instructions. AAs were detected using an automatic amino-acid analyser LC-3000 (Eppendorf-Bio-Tronik, Germany) according to Chinese national standard GB/T 18246–2000. FAs were detected and analysed using a gas chromatograph (GC, Agilent 7890A) method  with slight modifications. To better understand the characteristics of detected traits, histogram distribution, and Pearson’s correlation analyses were performed using R 3.2.3.
GBS library preparation and sequencing
Before carrying out the experiment, Perl package RestrictionDigest  was used to evaluate the number of restriction sites, and the AciI-BfaI restriction enzyme-cut combination was selected. After the DNA concentrations of 171 individuals (two parents and 169 full-sib progenies) were quantified to 20 ng/μL by Qubit 2.0 (Thermo Fisher Scientific, Waltham, MA), the GBS library was constructed with slight modifications to the reported method . In short, the DNA quantities of two parents were three-times that of the progeny. Adapter 1 and 2 containing the enzyme-cut site, barcode, and primer sequence matches AciI and BfaI, respectively. The length of the barcode ranged from 4 to 13 bp. A total of 9 kinds of adapter 1 and 10 kinds of adapter 2, were used (Additional file 1). Primers used for PCR amplification are also listed in Additional file 1. Pair-end 125-bp (125PE) sequencing of the library was performed on the Illumina HiSeq 2500 system (Illumina, San Diego, CA) according to the manufacturer’s recommendations.
Quality control and genotyping of sequencing data
Raw reads were first filtered to clean reads with the following criteria: reads with adapter sequences, unreliable reads (> 10% ‘N’ in the sequence), and low-quality reads (> 50% positions with quality < 5) were removed. After all reads were trimmed to 112 bp using a Perl script, Stacks software was used to distinguish and genotype the mapping population . First, the command “process_radtags -c -q -r --inline_inline --renz_1 aciI --renz_2 bfaI -i gzfastq” was run to demultiplex lane data to one fastq file per individual. In this step, reads of low quality, reads with wrong enzyme-cut sites, or with incorrect barcodes were discarded. Second, pair-end reads of each individual were aligned to an updated Pacific oyster genome (V9.2, in preparation) using the BWA aln parameter to produce a “sam” file . Based on the sam file, reads aligned to multiple genome positions were deleted by the “grep -v ‘XA:Z’” command, and badly aligned reads (align score < 20) were removed with a Perl script. After that, pstacks, cstacks, sstacks, and genotype procedure of Stacks were run in turn to build and genotype the haplotype marker across the full-sib family. In summary, only markers with a sequencing depth over 5× were reserved. We further filtered for the missing value rate (> 90%, that means genotyped in more than 152 individuals) for each of the marker. The output type was set as the CP population towards JoinMap4.1 software.
Linkage map construction and genome collinearity analysis
Prior to map construction, a segregation distortion test was performed. There are four categories of markers according to the segregation patterns: lm x ll (markers from the male parent with a 1:1 segregation ratio in the family), nn x np (markers from the female parent with a 1:1 segregation ratio in the family), and ef x eg and ab x cd (markers from both parents with 1:1:1:1 segregation ratio in the family). Markers showing significant segregation distortion with P < 0.01 in Chi-square goodness-of-fit tests were excluded. Sex-specific and averaged linkage maps were constructed using JoinMap4.1 . Markers were assigned to different LGs by a LOD threshold of 5.0. Marker distance was calculated using haldane function within the maximum likelihood (ML) algorithm. Sex averaged map was drawn using Mapchart 2.3 .
To evaluate the collinearity between the genetic map and genome assembly, the sequences of the markers, which were extracted from the tsv file produced by cstacks, were mapped onto the genome of C. gigas. We quantified the number of scaffold that contains at least two markers, and then analysed the distribution of the markers on the same scaffold.
QTL mapping of growth and nutritional traits
MapQTL6.0 software was employed to detect QTLs for all detected traits using interval mapping (IM) and restricted multiple QTL model (rMQM) . The LOD threshold was set by 3000-times-permutation tests at P < 0.05 for each trait. Mapping step size was set to 0.1 cM, and only QTLs with a LOD score exceeding the threshold were regarded as significant. To further confirm the QTL mapping results, we conducted association analysis between traits and genotypes using Tassel 5.0 with the MLM (PCA + K) method . Considering the results from both QTL mapping and association analysis, we identified markers that were significantly correlated with the traits, and identified candidate genes around these markers utilising the reference genome and annotation file of C. gigas.
Quantitative gene expression analysis of candidate genes
For each gene identified in the QTL, we firstly genotyped the full-family by the nearest marker of the gene. Then we compared the phenotypic values of the individuals that with different genotypes to ensure that the marker was indeed associated with the trait. Samples with the genotype that showed higher and lower phenotypic values were then named ‘high group’ and ‘low group’, respectively. Based on the genotyping results, we selected eight samples with the highest phenotypic values from ‘high group’ and eight samples with the lowest phenotypic values from ‘low group’, respectively. Further, we compared the gene mRNA expression levels between the two groups. The real-time reverse transcript polymerase chain reaction (qRT-PCR) experiment was performed using SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa, Japan) on ABI 7500 Fast Real time Thermal Cycler (Applied Biosystems, USA) according to the instructions. The expression levels of the gene were calculated with the 2-△△CT method. Besides, elongation factor (EF) primers were used as the internal control primers , and a random sample from the ‘high group’ was taken as the reference sample. All the primers used in this study were listed in Additional file 1.
Additional file 2 shows phenotypic values and variations for growth and glycogen, AA, and FA traits of 169 full-sib progenies. The distribution of different traits displayed various patterns with abundant variations, e.g. glycogen content ranged from 26.5 to 217.24 mg/g, representing a difference of more than 8 fold; the total AA content (All_AA) ranged from 27.38 to 56.78%; and the total FA content (All_FA) ranged from 16.9 to 58.9 mg/g (Additional file 3). Among the 18 detected AAs, glutamic acid (Glu) and Tau exhibited the highest levels, accounting for 15.21 and 11.94% of the total AA content, respectively (Additional file 4 A). Twenty-two fatty acids were detected, C16:0 exhibited the highest level, accounting for 47.82% of the total FA content (Additional file 4 B). We also conducted phenotypic correlation analysis for all traits. The results indicated that growth and glycogen were positively correlated with each other (Additional file 5 A), and the same result was also shown by AAs (Additional file 5 B) and FAs (Additional file 5 C) correlation pattern. Furthermore, by selecting representative traits from each major category, we found that glycogen, growth, and FA traits were positively correlated with each other, and they were all negatively correlated with AA traits (Additional file 5 D).
GBS library sequencing and marker genotyping
After reads were demultiplexed and cleaned using Perl scripts and process_radtags program, a total of 860.7 million (M) clean reads were generated. There were, on average, 7.54 and 4.3 M clean reads for parents and progeny, respectively, which corresponded to an actual sequencing depth of 27× on average. Using data from the parents, a catalogue of 296,187 markers was built. Based on this catalogue, 55,792 polymorphic markers with more than 5× sequencing depth were identified across all the samples. We calculated the missing rate for all the polymorphic markers, and found that most of the markers (59.63%) were missing in more than half of the population (Additional file 6). To insure the quality of the markers, we only used 7075 markers with the missing rate of less than 10%, and further filtered the markers with the segregation distortion test (P < 0.01). In total, 5094 markers were retained for the genetic map construction.
High-resolution genetic map construction and analysis of recombination rate
Overall, 5024 markers were successfully assigned to 10 LGs by LOD 5.0 using JoinMap4.1 (Table 1). Detailed information of the markers can be seen in Additional file 7. For sex-specific maps, the paternal map consisted of 2993 markers with a total length of 1765.49 cM, whereas the maternal map consisted of 2963 markers with a total genetic map length of 1931.26 cM (Additional file 8). Using 932 shared markers (ef x eg and ab x cd, which are heterozygous in both parents), a consensus genetic map was constructed. Different types of markers were evenly distributed across the map, with the exception of LG6 and LG10. LG6 contained more paternal markers, while LG10 consisted mainly of maternal markers (Fig. 1). The genetic length of each LG ranged from 148.32 (LG5) to 312.43 (LG4) cM, with the total map length of 1982.07 cM. All 5024 markers of the consensus map occupied 2931 different positions, and the unique marker interval (interval only for markers occupied different positions on the map) ranged from 0.53 (LG1) to 1.22 cM (LG10), with an average interval of 0.68 cM (Table 1).
Marker distribution pattern on each LG showed a similar slope among different LGs, which indicated similar recombination rates of different LGs and evenly spaced markers on the consensus map (Fig. 2). Similar recombination rates were also confirmed by comparisons between sex-specific maps. First, the maternal map showed slightly higher values than the paternal map concerning both average map length and average marker interval (the average female-to-male ratio was 1.14 and 1.18, respectively), but these differences were not significant (t-test, P > 0.05) (Additional file 8). Second, as marker interval was closely linked to recombination rate, we used the interval of 932 shared markers to evaluate recombination rates between the sex-specific maps. Recombination rates of shared markers between sex-specific maps and the consensus map of C. gigas revealed similar recombination rates and high consistency between sex-specific maps (Fig. 3). The average female-to-male ratio for recombination rates of shared markers among different LGs was 1.07, ranging from 0.43 (LG6) to 1.41 (LG2) (Additional file 8). Recombination rates of the female-specific map were slightly higher than those of the male, although no significance was found (t-test, P > 0.05). High consistency between sex-specific maps laid the basis for generating a high-qualified consensus map.
Genome collinearity analysis
Based on the high-qualified consensus map, we used the genome position of the markers to evaluate the collinearity between the map and genome of C. gigas. Results showed that all 5024 markers were aligned to 1574 scaffold, covering 68.81% of the oyster genome. Overall, 3809 markers with a unique genome position and an alignment length of over 100 bp were used for further analysis. These 3809 markers were distributed on 1303 scaffold of the genome, and 761 of those scaffold had more than one marker. Overall, 86.47% of the scaffold with at least two markers (658 of 761) were located on the same LG of the consensus map (Table 2). Collinearity analysis of the scaffold position on the genome towards the marker position on the consensus map showed that most plots were on the diagonal or adjacent position, which indicated that markers of the same scaffold were distributed on the same or adjacent position of the consensus map (Fig. 4).
QTL mapping of growth and nutritional traits and candidate genes expression levels analyses
We detected a total of 41 QTLs for 17 of the 46 traits (Table 3). For growth traits, 15 QTLs were detected for SL, BW, and STW, with the explained phenotype variation (PVE) ranging from 8.9 to 15.2% (Table 3, Fig. 5a, Additional file 9). The correlation coefficient between BW and STW was 0.8 (P < 0.01, Additional file 5 A), and all of the QTLs were on LG4 and LG9 (Table 3). Two protein-coding genes were identified in qBW_2 (Table 4, Fig. 5a). Bone morphogenetic protein 1 (BMP1) is a metalloproteinase that can regulate the formation of the extracellular matrix and dorsal-ventral patterning in early embryos [46, 47]; While glucose-dependent insulinotropic receptor (GIPR), which is responsible for insulin secretion, has been reported to be associated with obesity and weight gain [48,49,50]. Two genes associated with growth (Tables 3-4) were also identified in on QTL of STW on LG4 (qSTW_4). Alpha-amylase (AMY) has been widely reported to affect growth, e.g. AMY inhibits the growth of Porphyromonas gingivalis , and polymorphisms in AMY was correlated with growth traits in Litopenaeus vannamei ; While fibroblast growth factor receptor (FGFR), was found to be significantly correlated with four growth traits in common carp .
For glycogen, only one QTL on LG2 with a PVE of 8.6% was detected (Table 3, Fig. 5b). This QTL interval identified many zinc finger proteins (ZFPs), such as ZFP729, ZFP668, and ZFP37 (Table 4). ZFP is a transcription factor that functions in many important physiological processes by regulating gene transcription . Of the three types of identified ZFPs, ZFP37 might be the key gene for its role in the regulation of spermiogenesis . Glycogen metabolism was widely reported to regulate sperm motility and capacitation [56,57,58,59]. Therefore, we suggest that ZFP37 might be an upstream regulatory gene in glycogen metabolism, and may regulate spermiogenesis by influencing the glycogen content.
Seventeen QTLs were detected for 7 AAs, including Glu, threonine (Thr), glycine (Gly), Tau, histidine (His), proline (Pro), and valine (Val) (Table 3, Fig. 5c, Additional file 9). The correlation coefficient between Glu, Thr, and Val exceeded 0.94 (P < 0.01, Additional file 5 B), and similar QTLs ranging from 25.152 to 34.413 cM were detected for these AAs on LG2 (Table 3, Fig. 5c, Additional file 9). One key gene identified in the QTL was delta-1-pyrroline-5-carboxylate synthase (P5CS) (Table 4). P5CS catalyses the reduction of Glu to delta1-pyrroline-5-carboxylate (P5C), a critical step in the de novo biosynthesis of Pro, ornithine, and arginine , and is thus a vital enzyme in Glu and other AA metabolism. Two genes related to glycometabolism were identified in qTau_2 and qTau_4, respectively (Table 4). Succinate dehydrogenase (SDH) can catalyse the oxidation of succinate to fumarate, with the reduction of ubiquinone to ubiquinol during the citric acid cycle ; While glycogen phosphorylase (GP) catalyses the rate-limiting step in glycogenolysis in animals by releasing glucose-1-phosphate from the terminal alpha-1,4-glycosidic bond . Tau was reported to have hypoglycaemic effect in previous studies [8, 63, 64]. Therefore, we suggest that Tau may affect glucose levels by regulating enzymes in the metabolic pathway. Furthermore, we identified a neutral and basic amino acid transport protein (rBAT) in qTau_4 with a PVE of 9% (Tables 3-4). rBAT functions as cysteine transporter [65,66,67]. As a sulfur-containing amino acid, cysteine is a substrate of tau synthesis. Therefore, rBAT may influence Tau levels by regulating cysteine transport.
A total of eight QTLs were detected for six FAs (C17:1, C20:1ω9, C20:2, C20:3ω6, C20:4ω6, and C22:0). Eight QTLs detected for FAs were distributed on six LGs, which indicated that different genes might be associated with different kinds of FAs (Table 3, Fig. 5d, Additional file 9). For C20:4ω6, one QTL located on LG4 and two QTLs located on LG9 were detected (Table 3, Fig. 5d). Two genes identified in qC20:4ω6_3 were shown to be associated with body weight or fat mass (Table 4). Dual-specificity tyrosine-phosphorylation-regulated kinase (DYRK), can prevent metabolic disorders like obesity by reducing fat mass ; while a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS), was widely reported to be associated with adipogenesis, fat content, and obesity [69,70,71,72]. qC20:1ω9, located on LG6, could explain 10.2% of the PVE for C20:1ω9 (Table 3, Additional file 9). One candidate gene identified in this QTL was growth hormone secretagogue receptor (GHSR). GHSR polymorphisms are associated with carcass traits in sheep, and individuals with a mutant ‘TC’ genotype exhibit greater levels of abdominal fat . For C17:1, the only QTL qC17:1 was located on LG5 with a PVE of 8.9% (Table 3, Additional file 9). The gene identified in this QTL was phospholipid scramblase (PLSCR). PLSCR proteins are responsible for the translocation of phospholipids between the two monolayers of a membrane lipid bilayer, and are critical to the normal regulation of lipid metabolism [74,75,76].
To further study the relationship between the candidate genes and the traits, we analysed the mRNA expression levels among samples with different phenotypic values. Of all 13 identified genes, we designed 9 primers successfully. The results showed that six genes exhibited significant expression difference (P < 0.05) between groups with different genotypes and phenotypic values (Fig. 6). Boxplot from Fig. 6 also showed that phenotypic values between different genotypes were partly overlapped with the other, which may reflect the minor effect of single gene to quantitative traits. Further, for the six differently expressed genes, BMP1 and AMY were identified for growth traits, P5CS and GP were for AA, while ADAMTS and DYRK were for FA. P5CS and GP genes showed high mRNA expression levels in ‘low group’, but low mRNA expression levels in ‘high group’; While all the other four genes showed higher mRNA expression levels in ‘high group’.
High-quality genetic map of C. gigas
After applying rigorous filter criteria (missing rate < 10%, segregation distortion P < 0.01), 5094 markers were used to construct the genetic map, of which, 5024 (98.62%) were successfully assigned to 10 LGs with the unique marker density of 0.68 cM (Table 1). Furthermore, even marker distribution (Table 1, Fig. 1) and similar patterns of recombination rates between sex-specific maps (Fig. 3) all confirmed the high quality of the consensus map. The Pacific oyster was reported to be of high heterozygosity and polymorphism , which may be the main reason for high missing rate for the enzyme-cut markers. It’s intriguing to see that LG6 and LG10 contain almost exclusively paternal and maternal markers, respectively. The Pacific oyster showing widespread segregation distortion in inbred F2 or F3 families, carries a large load of highly deleterious, primarily recessive mutation [78,79,80]. To some extent, this inbreeding depression was lineage- and environment-specific . Based on which, we speculate that female (sampled from Changli) LG6 and male (sampled from Qingdao) LG10 may bear some structural variation towards the breeding environment, which will cause lethal zygotes and further result in paternal-marker dependent LG6 and maternal-marker dependent LG10. The length of the consensus map was 1982.07 cM, which was longer than previously reported maps with the length ranging from 558.2 to 1084.3 cM [22, 23, 29, 31, 35, 40]. Different mapping methods in JoinMap4.1 may be the major reason causing this difference. Multi-point ML algorithm with haldnane function would increase the map length without considering double recombination events . However, ML algorithm runs more quickly than regression algorithm when there’re more than 400 markers on one LG (a few hours versus more than 1 month) and ML yields better fitting maps than regression method in most cases . Longer map length was also observed in other aquatic species, for which the genetic map was constructed with thousands of markers mainly from NGS [82, 83]. In summary, we have constructed a high-quality genetic map for C. gigas with the highest number of markers, the greatest marker density, and the largest genome coverage so far. This highly qualified genetic map will provide a basis for genome assembly and fine mapping of vital economic traits.
Recently, high-density genetic maps have been an important tool for genome scaffold assembly of aquatic species. For example, 246 genome scaffold, which represented 95.78% of the assembled genomic sequences of Paralichthys olivaceus, were accurately anchored onto 24 LGs by a GBS-based high-resolution genetic map . Similarly, 2839 of 2963 (95.8%) genome scaffold of Cyprinus carpio were mapped to unique LGs by a micro-array based high-density genetic map . In this study, all 5024 markers were mapped to 1524 genome scaffold, covering 68.81% of the genome. Further analysis showed that only 13.53% of the scaffold with more than one marker distributed on different LGs of the genetic map (Table 2). Our results not only showed a good collinearity between the consensus map and the oyster genome but also indicated a better genome assembly of C. gigas. Accordingly, we successfully assembled a chromosome-level genome framework covering about 70% of the genome sequence of C. gigas, providing a key basis for further comparative genomics and fine QTL mapping analysis.
QTL mapping and identification of candidate genes for growth and nutritional traits
In this study, we utilised a ddGBS-based high-quality linkage map to identify genes controlling both growth and nutritional traits of C. gigas. We not only identified known growth-related genes like AMY and BMP1, but also identified new genes participating in AA and FA metabolism and transportation, like P5CS, rBAT, PLSCR, and DYRK. Candidate genes identified for growth were different from the results before . Variation in growth may be influenced by many loci with small effects . Compared with before, different mapping family and enzyme-cut experimental design in this study may produce diverse segregation loci and result in disparate candidate genes. Polymorphism of glycogen phosphorylase was reported to be associated with the glycogen content in C. gigas , while it was identified in the QTL of Tau in this study. Considering the significant negative correlation between Tau and glycogen (Additional file 5 D), this gene may play a key role in the metabolism regulation of both Tau and glycogen. Genes identified in FA QTLs were reported either to regulate lipid metabolism or to be associated with body weight (Table 4), which was reasonable considering the correlation between fat and body weight. Further, we conducted gene expression analysis for the identified genes. Six genes showed gene expression difference between groups with different genotypes and phenotypic values, which indicated that the gene may influence the trait by transcriptional regulation. Meanwhile, genes identified for AA traits exhibited different expression patterns compared to growth and FA traits, which reflected diverse regulatory mechanism underlying different candidate genes and traits. These six genes will be key targets for molecular breeding of important economic traits. One underlying concern for this study is the overestimation of QTL mapping resulted from low sampling size (Beavis Effect ). In this study, we not only conducted strict data filter for each step towards QTL mapping but also verified gene function on transcriptional levels, both of which will convince our findings strongly. Further studies based on the construction of bigger populations or families with more divergent phenotypic difference will also improve the power of QTL mapping.
In this study, we used ddGBS method to generate a high-density genetic map with the most markers for C. gigas. Based on the high collinearity between the linkage map and the oyster genome, we successfully assembled nearly 70% of the oyster genome sequence onto chromosomes for the first time. Besides, we provided valuable markers and candidate genes for growth, and firstly for AA and FA traits of C. gigas by QTL mapping. Six genes that can influence the traits by transcriptional regulation may be causative genes for further studies concerning qualified traits improvement. In a whole, our findings will promote genetic dissection and molecular breeding of important economic traits in C. gigas.
Total amino acids
Total fatty acids
Quantitative trait loci
Weight of the soft tissue
Mann R. Exotic species in aquaculture. Oceanus. 1979;22(1):29–35.
Guo X. Use and exchange of genetic resources in molluscan aquaculture. Rev Aquac. 2009;1(3–4):251–9.
Wu Y, Sun H. Progress of functional products from oyster. Hebei Fisheries. 2007;8:6–9.
Shen S, Xu C, Tong D. Nutrition and health function of oyster and its exploitation and utilization. J Heibei Agric Sci. 2009;10:79–81.
El Idrissi A, Shen CH, L'Amoreaux WJ. Neuroprotective role of taurine during aging. Amino Acids. 2013;45(4):735–50.
Menzie J, Pan C, Prentice H, Wu J-Y. Taurine and central nervous system disorders. Amino Acids. 2014;46(1):31–46.
Marcinkiewicz J, Kontny E. Taurine and inflammatory diseases. Amino Acids. 2014;46(1):7–20.
Chen W, Guo J, Zhang Y, Zhang J. The beneficial effects of taurine in preventing metabolic syndrome. Food Funct. 2016;7(4):1849–63.
Yamaguchi R, Shirai T, Umeda M, Shiba T, Igarashi M, Asada T. The relationship between thrombotic diseases and the activities of EPA and DHA. Thromb Haemost. 1983;50(1):136.
Dickinson KM, Delaney CL, Allan R, Spark I, Miller MD. Validation of a brief dietary assessment tool for estimating dietary EPA and DHA intake in Australian adults at risk of cardiovascular disease. J Am Coll Nutr. 2015;34(4):333–9.
Ward RD, English LJ, McGoldrick DJ, Maguire GB, Nell JA, Thompson PA. Genetic improvement of the Pacific oyster Crassostrea gigas (Thunberg) in Australia. Aquac Res. 2000;31(1):35–44.
Degremont L, Bedier E, Boudry P. Summer mortality of hatchery-produced Pacific oyster spat (Crassostrea gigas). II. Response to selection for survival and its influence on growth and yield. Aquaculture. 2010;299(1–4):21–9.
Li Q, Wang Q, Liu S, Kong L. Selection response and realized heritability for growth in three stocks of the Pacific oyster Crassostrea gigas. Fish Sci. 2011;77(4):643–8.
Wang QZ, Li Q, Kong LF, Yu RH. Response to selection for fast growth in the second generation of Pacific oyster (Crassostrea gigas). J Ocean Univ China. 2012;11(3):413–8.
Zhong X, Li Q, Yu H, Kong L. Development and validation of single-nucleotide polymorphism markers in the Pacific oyster, Crassostrea gigas, using high-resolution melting analysis. J World Aquacult Soc. 2013;44(3):455–65.
Martino RC, da Cruz GM. Proximate composition and fatty acid content of the mangrove oyster Crassostrea rhizophorae along the year seasons. Braz Arch Biol Technol. 2004;47(6):955–60.
Saito H. Lipid and FA composition of the pearl oyster Pinctada fucata martensii: influence of season and maturation. Lipids. 2004;39(10):997–1005.
She Z, Li L, Qi H, Song K, Que H, Zhang G. Candidate gene polymorphisms and their association with glycogen content in the Pacific oyster Crassostrea Gigas. PLoS One. 2015;10(5):e0124401.
Hedgecock D, Hubert S, Li G, Bucklin K. A genetic linkage map of 100 microsatellite markers for the Pacific oyster Crassostrea gigas. J Shellfish Res. 2002;21(1):381.
Yu ZN, Guo XM. Genetic linkage map of the eastern oyster Crassostrea virginica Gmelin. Biol Bull. 2003;204(3):327–38.
Guo XM, Wang Y, Yu ZN. Physical and linkage mapping in the eastern oyster (Crassostrea Virginica Gmlin). J Shellfish Res. 2004;23(1):294–5.
Hubert S, Hedgecock D. Linkage maps of microsatellite DNA markers for the pacific oyster Crassostrea gigas. Genetics. 2004;168(1):351–62.
Li L, Guo XM. AFLP-based genetic linkage maps of the Pacific oyster Crassostrea gigas Thunberg. Mar Biotechnol. 2004;6(1):26–36.
Yu ZN, Guo XM. Identification and mapping of disease-resistance QTLs in the eastern oyster, Crassostrea virginica Gmelin. Aquaculture. 2006;254(1–4):160–70.
Hedgecock D, Perry GML, Voigt ML. Mapping heterosis QTL in the Pacific oyster Crassostrea gigas. Aquaculture. 2007;272:S267–8.
Lallias D, Beaumont A, Haley C, Boudry P, Heurtebise S, Lapegue S. A first-generation genetic linkage map of the European flat oyster Ostrea edulis (L.) based on AFLP and microsatellite markers. Anim Genet. 2007;38(6):560–8.
Perryt GNL, Voigt M-L, Hedgecock D. Mapping QTL controlling growth and body size in the Pacific oyster. J Shellfish Res. 2008;27(4):1040.
Lallias D, Gomez-Raya L, Haley CS, Arzul I, Heurtebise S, Beaumont AR, Boudry P, Lapegue S. Combining two-stage testing and interval mapping strategies to detect qtl for resistance to bonamiosis in the european flat oyster Ostrea edulis. Mar Biotechnol. 2009;11(5):570–84.
Sauvage C, Boudry P, de Koning DJ, Haley CS, Heurtebise S, Lapegue S. QTL for resistance to summer mortality and OsHV-1 load in the Pacific oyster (Crassostrea gigas). Anim Genet. 2010;41(4):390–9.
Plough LV, Hedgecock D. Quantitative trait locus analysis of stage-specific inbreeding depression in the Pacific oyster Crassostrea gigas. Genetics. 2011;189(4):1473.
Guo X, Li Q, Wang QZ, Kong LF. Genetic mapping and QTL analysis of growth-related traits in the Pacific oyster. Mar Biotechnol. 2012;14(2):218–26.
Stick DA, Camara MD. Identification and mapping of growth-related QTL using microsatellite and AFLP markers for the Pacific oyster. J Shellfish Res. 2012;31(1):349.
Ge JL, Li Q, Yu H, Kong LF. Identification and mapping of a SCAR marker linked to a locus involved in shell pigmentation of the Pacific oyster (Crassostrea gigas). Aquaculture. 2014;434:249–53.
Zhong X, Li Q, Guo X, Yu H, Kong L. QTL mapping for glycogen content and shell pigmentation in the Pacific oyster Crassostrea gigas using microsatellites and SNPs. Aquac Int. 2014;22(6):1877–89.
Hedgecock D, Shin G, Gracey AY, Van Den Berg D, Samanta MP. Second-generation linkage maps for the Pacific oyster Crassostrea gigas reveal errors in assembly of genome scaffolds. G3-Genes Genomes Genetics. 2015;5(10):2007–19.
Wang J, Li L, Zhang G. A high-density snp genetic linkage map and QTL analysis of growth-related traits in a hybrid family of oysters (Crassostrea gigas x Crassostrea angulata) using genotyping-by-sequencing. G3-Genes Genomes Genetics. 2016;6(5):1417–26.
Liu S, Qi LI, Hong YU, Kong L. Relationship between single nucleotide polymorphism of glycogen synthase gene of Pacific oyster Crassostrea gigas and its glycogen content. J Ocean Univ China. 2017;16(1):168–74.
Wang J, Li Q, Zhong X, Song J, Kong L, Yu H. An integrated genetic map based on EST-SNPs and QTL analysis of shell color traits in Pacific oyster Crassostrea gigas. Aquaculture. 2018;492:226–36.
Yang XH, Guo YQ, Yan JB, Zhang J, Song TM, Rocheford T, Li JS. Major and minor QTL and epistasis contribute to fatty acid compositions and oil concentration in high-oil maize. Theor Appl Genet. 2010;120(3):665–78.
Wang J, Li L, Qi H, Dua X, Zhang G. RestrictionDigest: a powerful Perl module for simulating genomic restriction digests. Electron J Biotechnol. 2016;21:36–42.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22(11):3124–40.
Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Van Ooijen JW. Multipoint maximum likelihood mapping in a full-sib family of an outbreeding species. Genet Res. 2011;93(5):343–9.
Voorrips RE. MapChart: software for the graphical presentation of linkage maps and QTLs. J Hered. 2002;93(1):77–8.
Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.
Wardle FC, Welch JV, Dale L. Bone morphogenetic protein 1 regulates dorsal-ventral patterning in early Xenopus embryos by degrading chordin, a BMP4 antagonist. Mech Dev. 1999;86(1–2):75–85.
Hopkins D, Keles S. Ds: the bone morphogenetic protein 1/Tolloid-like metalloproteinases. Matrix Biol. 2007;26(7):508–23.
Saxena R, Hivert M-F, Langenberg C, Tanaka T, Pankow JS, Vollenweider P, Lyssenko V, Bouatia-Naji N, Dupuis J, Jackson AU, et al. Genetic variation in GIPR influences the glucose and insulin responses to an oral glucose challenge. Nat Genet. 2010;42(2):142–U175.
Ono S, Suzuki Y, Fukui N, Sawamura K, Sugai T, Watanabe J, Tsuneyama N, Someya T. GIPR gene polymorphism and weight gain in patients with schizophrenia treated with olanzapine. J Neuropsychiatry Clin Neurosci. 2015;27(2):162–4.
Ceperuelo-Mallafre V, Duran X, Pachon G, Roche K, Garrido-Sanchez L, Vilarrasa N, Tinahones FJ, Vicente V, Pujol J, Vendrell J, et al. Disruption of GIP/GIPR axis in human adipose tissue is linked to obesity and insulin resistance. J Clinical Endocrinol Metab. 2014;99(5):E908–19.
Ochiai A, Harada K, Hashimoto K, Shibata K, Ishiyama Y, Mitsui T, Tanaka T, Taniguchi M. Alpha-amylase is a potential growth inhibitor of Porphyromonas gingivalis, a periodontal pathogenic bacterium. J Periodontal Res. 2014;49(1):62–8.
Jingjing XIN, Xiaolin LIU, Xilian LI, Xianzong W, Hao H, Jianhai X. PCR-RFLP polymorphism of alpha-amylase gene and its association with growth traits of Litopenaeus vannamei. Acta Oceanol Sin. 2011;33(3):124–30.
Xiao T, Lu C, Li C, Cao D, Cheng L, Sun X. Correlation analysis of microsatellite markers derived from scale genes (ant, eda, edar, fgfr) with four growth traits in common carp (Cyprinus carpio). J Fishery Sci of China. 2014;21(5):883–92.
Li WT, He M, Wang J, Wang YP. Zinc finger protein (ZFP) in plants-a review. Plant Omics. 2013;6(6):474–80.
Burke PS, Wolgemuth DJ. ZFP-37, a new murine zinc finger encoding gene, is expressed in a developmentally regulated pattern in the male germ line. Nucleic Acids Res. 1992;20(11):2827–34.
Albarracin JL, Fernandez-Novell JM, Ballester J, Rauch MC, Quintero-Moreno A, Pena A, Mogas T, Rigau T, Yanez A, Guinovart JJ, et al. Gluconeogenesis-linked glycogen metabolism is important in the achievement of in vitro capacitation of dog spermatozoa in a medium without glucose. Biol Reprod. 2004;71(5):1437–45.
Aparicio IM, Bragado MJ, Gil MC, Garcia-Herreros M, Gonzalez-Fernandez L, Tapia JA, Garcia-Marin LJ. Porcine sperm motility is regulated by serine phosphorylation of the glycogen synthase kinase-3 alpha. Reproduction. 2007;134(3):435–44.
Guo TB, Chan KC, Hakovirta H, Xiao Y, Toppari J, Mitchell AP, Salameh WA. Evidence for a role of glycogen synthase kinase-3 beta in rodent spermatogenesis. J Androl. 2003;24(3):332–42.
Jack S, Acharya R, Kline D, Olson GE, Vijayaraghavan S. Glycogen synthase kinase-3 phosphorylation - evidence for a role for the phosphoinositide 3-kinase signaling pathway in sperm motility regulation. Biol Reprod. 2000;62:206–7.
Zhuang GQ, Li B, Guo HY, Liu JL, Chen F. Molecular cloning and characterization of P5CS gene from Jatropha curcas L. Afr J Biotechnol. 2011;10(66):14803–11.
Oyedotun KS, Lemire BD. The quaternary structure of the Saccharomyces cerevisiae succinate dehydrogenase - homology modeling, cofactor docking, and molecular dynamics simulation studies. J Biol Chem. 2004;279(10):9424–31.
Bacca H, Huvet A, Fabioux C, Daniel JY, Delaporte A, Pouvreau S, Van Wormhoudt A, Moal J. Molecular cloning and seasonal expression of oyster glycogen phosphorylase and glycogen synthase genes. Comp Biochem Physiol B-Biochem Mol Biol. 2005;140(4):635–46.
Gao Y, Guo J, Zhang Y, Chen W. Research progress of hypoglycaemic effect of taurine. J Chin Inst Food Sci Technol. 2016;16(1):202–10.
Murakami S. Role of taurine in the pathogenesis of obesity. Mol Nutr Food Res. 2015;59(7):1353–63.
Bisceglia L, Stanziale P. Novel human pathological mutations. Gene symbol: SLC3A1. Disease: cystinuria. Hum Genet. 2010;127(4):473.
Markazi S, Kheirollahi M, Doosti A, Mohammadi M, Koulivand L. A novel mutation in SLC3A1 gene in patients with cystinuria. Iran J Kidney Dis. 2016;10(1):44–7.
Mizukami K, Raj K, Giger U. Feline Cystinuria caused by a missense mutation in the SLC3A1 gene. J Vet Intern Med. 2015;29(1):120–5.
Song W-J, Song E-AC, Jung M-S, Choi S-H, Baik H-H, Jin BK, Kim JH, Chung S-H. Phosphorylation and inactivation of glycogen synthase kinase 3 beta (GSK3 beta) by dual-specificity tyrosine phosphorylation-regulated kinase 1A (Dyrk1A). J Biol Chem. 2015;290(4):2321–33.
Bauters D, Scroyen I, Deprez-Poulain R, Lijnen HR. ADAMTS5 promotes murine adipogenesis and visceral adipose tissue expansion. Thromb Haemost. 2016;116(4):694–704.
Bauters D, Spincemaille P, Geys L, Cassiman D, Vermeersch P, Bedossa P, Scroyen I, Lijnen HR. ADAMTS5 deficiency protects against non-alcoholic steatohepatitis in obesity. Liver Int. 2016;36(12):1848–59.
Kumar S, Chen M, Li Y, Wong FHS, Thiam CW, Hossain MZ, Poh KK, Hirohata S, Ogawa H, Angeli V, et al. Loss of ADAMTS4 reduces high fat diet-induced atherosclerosis and enhances plaque stability in ApoE(−/−) mice. Sci Rep. 2016;6:31130.
Ql Y, Zhang X, Wang Y, Wang D, Guo Z, Liu P. The expression of ADAMTS2 and collagen genes in muscle tissue and its relationship with meat quality characters in cattle. J Yunnan Agric Univ. 2014;29(2):173–8.
Bahrami A, Miraei-Ashtiani SR, Mehrabani-Yeganeh H. Associations of growth hormone secretagogue receptor (GHSR) genes polymorphisms and protein structure changes with carcass traits in sheep. Gene. 2012;505(2):379–83.
Wiedmer T, Zhao J, Li LL, Zhou QS, Hevener A, Olefsky JM, Curtiss LK, Sims PJ. Adiposity, dyslipidemia, and insulin resistance in mice with targeted deletion of phospholipid scramblase 3 (PLSCR3). Proc Natl Acad Sci U S A. 2004;101(36):13296–301.
Sahu SK, Aradhyam GK, Gummadi SN. Calcium binding studies of peptides of human phospholipid scramblases 1 to 4 suggest that scramblases are new class of calcium binding proteins in the cell. Biochimica Et Biophysica Acta-General Subjects. 2009;1790(10):1274–81.
Sahu SK, Gummadi SN, Manoi N, Aradhyam GK. Phospholipid scramblases: an overview. Arch Biochem Biophys. 2007;462(1):103–14.
Zhang G, Fang X, Guo X, Li L, Luo R, Xu F, Yang P, Zhang L, Wang X, Qi H, et al. The oyster genome reveals stress adaptation and complexity of shell formation. Nature. 2012;490(7418):49–54.
Launey S, Hedgecock D. High genetic load in the Pacific oyster Crassostrea gigas. Genetics. 2001;159(1):255–65.
Plough LV, Hedgecock D. Quantitative trait locus analysis of stage-specific inbreeding depression in the Pacific oyster Crassostrea gigas. Genetics. 2011;189(4):1473–86.
Bucklin KA. Analysis of the genetic basis of inbreeding depression in the Pacific oyster Crassostrea gigas. Davis: PhD Dissertation, University of California; 2003. p. 140.
Plough LV. Environmental stress increases selection against and dominance of deleterious mutations in inbred families of the Pacific oyster Crassostrea gigas. Mol Ecol. 2012;21(16):3974–87.
Peng W, Xu J, Zhang Y, Feng J, Dong C, Jiang L, Feng J, Chen B, Gong Y, Chen L, et al. An ultra-high density linkage map and QTL mapping for sex and growth-related traits of common carp (Cyprinus carpio). Sci Rep. 2016;6:26693.
Shao C, Niu Y, Rastas P, Liu Y, Xie Z, Li H, Wang L, Jiang Y, Tai S, Tian Y, et al. Genome-wide SNP identification for the construction of a high-resolution genetic map of Japanese flounder (Paralichthys olivaceus): applications to QTL mapping of Vibrio anguillarum disease resistance and comparative genomic analysis. DNA Res. 2015;22(2):161–70.
Xu S. Theoretical basis of the Beavis effect. Genetics. 2003;165(4):2259–68.
The authors thank all members of the laboratory for valuable discussions.
This research was supported by the National Natural Science Foundation of China (31530079 to GFZ and 31572620 to LL, in oyster cultivation, sample collection and experiments), Strategic Priority Research Program of “Western Pacific Ocean System: Structure, Dynamics and Consequences” (XDA586 11020305, in data analysis and interpretation), and the Earmarked Fund for Modern Agro-industry Technology Research System (CARS-49, in manuscript drafting and publication).
Availability of data and materials
Sequencing raw data were deposited at NCBI under the accession number PRJNA389216.
Ethics approval and consent to participate
The parent oysters were sampled from the intertidal zone of Changli and Qingdao, respectively. The F1 progenies were cultured in a farm in Qingdao before attachment, and cultured in the sea after attachment of the larvae. There’s no specific permissions for the sampling of the oysters. And all the experiments were carried out with approval from Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Adaptors and primers used in this study. (XLS 34 kb)
Histogram of growth and nutritional traits. Individual panels show (A) growth and glycogen, (B) amino acid, and (C) fatty acid. Full names for the abbreviation are as follows: shell height (SH), shell length (SL), shell width (SW), body weight (BW), weight of the soft tissue (STW), total amino acids (All_AA), glutamic acid (Glu), taurine (Tau), aspartic acid (Asp), lysine (Lys), alanine (Ala), leucine (Leu), glycine (Gly), proline (Pro), arginine (Arg), valine (Val), serine (Ser), isoleucine (Ile), threonine (Thr), phenylalanine (Phe), tyrosine (Tyr), histidine (His), cysteine (Cys), methionine (Met), total fatty acids (All_FA), C20:5ω3 (EPA), and C22:6ω3 (DHA). (TIF 552 kb)
Phenotype values and variance. Full names for the abbreviation are shown in Additional file 2. (XLS 112 kb)
Amino acid and fatty acid composition. (A) Amino acid composition, (B) Fatty acid composition. Full names for the abbreviation are shown in Additional file 2. (TIF 625 kb)
Correlation analysis of growth and nutritional traits. (A) Correlation analysis of (A) growth and glycogen, (B) amino acid, (C) fatty acid, and (D) representative traits. Digits above the diagonal show correlation coefficient, while patterns below the diagonal show scatter plots of correlation. Significance correlation (P < 0.05) and (P < 0.01) was indicated by * and **, respectively. Full names for the abbreviation are shown in Additional file 2. (TIF 2998 kb)
Missing rate statistics for all polymorphic markers of the mapping family. (XLS 25 kb)
Detailed information of 5024 haplotype markers used for genetic map construction. (XLS 1260 kb)
Summary of the sex-specific maps. (XLS 35 kb)
About this article
Cite this article
Li, C., Wang, J., Song, K. et al. Construction of a high-density genetic map and fine QTL mapping for growth and nutritional traits of Crassostrea gigas. BMC Genomics 19, 626 (2018). https://doi.org/10.1186/s12864-018-4996-z
- Crassostrea gigas
- Genetic map
- Quantitative trait loci (QTL) mapping