Genome-Wide Mapping of Growth-Related Quantitative Trait Loci in Orange-Spotted Grouper (Epinephelus coioides) Using Double Digest Restriction-Site Associated DNA Sequencing (ddRADseq)

Mapping of quantitative trait loci (QTL) is essential for the discovery of genetic structures that related to complex quantitative traits. In this study, we identified 264,072 raw SNPs (single-nucleotide polymorphisms) by double digest restriction site associated DNA sequencing (ddRADseq), and utilized 3029 of these SNPs to construct a genetic linkage map in orange-spotted grouper (Epinephelus coioides) using a regression mapping algorithm. The genetic map contained 24 linkage groups (LGs) spanning a total genetic distance of 1231.98 cM. Twenty-seven significant growth-related QTLs were identified. Furthermore, we identified 17 genes (fez2, alg3, ece2, arvcf, sla27a4, sgk223, camk2, prrc2b, mchr1, sardh, pappa, syk, tert, wdrcp91, ftz-f1, mate1 and notch1) including three (tert, ftz-f1 and notch1) that have been reported to be involved in fish growth. To summarize, we mapped growth-related QTLs in the orange-spotted grouper. These QTLs will be useful in marker-assisted selection (MAS) efforts to improve growth-related traits in this economically important fish.


Introduction
Orange-spotted grouper, Epinephelus coioides (Epinephelinae, Serranidae), a protogynous hermaphrodite species, is mainly distributed in the Indo-West-Pacific region [1]. It is an economically important aquaculture species in many Asian countries, especially in China, due to its desirable taste and as a source of nutrition [2]. Orange-spotted grouper has become an important edible fish species in live fish markets and is an important cultured fish for sale in markets in southeast China. Market demand for the orange-spotted grouper has prompted the development of fish families and populations characterized by lower food consumption and higher growth rates. The application of marker-assisted selection (MAS) to orange-spotted grouper will be useful to improve important traits, such as feed conversion rate, meat quality and disease resistance. These traits have a major effect on productivity and profitability.
High-density genetic linkage maps and quantitative trait loci (QTL) mapping provide a framework for the MAS program. At present, genetic linkage maps have been constructed in over 28 fish species and economic traits have been mapped in at least 13 fish species [3]. These include the Nile tilapia (orechromis niloticus) [4], Atlantic salmon (Salmo salar) [5], rainbow trout (Oncorhynchus mykiss) [6], and channel catfish (Ictalurus punctatus) [7]. Microsatellite-based linkage maps have been reported for the white grouper (Epinephelus aeneus) [8] and kelp grouper (Epinephelus bruneus) [9]. The density of the genetic linkage maps is typically determined by the marker type. Single-nucleotide polymorphisms (SNPs), the most stable and abundant form of genetic marker, are the ideal marker type for construction of high-density genetic linkage maps [10][11][12]. The advent of next-generation sequencing (NGS) has made it possible to discover thousands of SNPs dispersed throughout the genome in a single procedure, even if little genetic information is available on the species [13]. A SNP-based high-density genetic linkage map has been constructed for orange-spotted grouper by our lab [14], but QTL mapping was not reported. This is the first report of a SNP-based high-resolution genetic linkage map of the orange-spotted grouper developed using double digest restriction site associated DNA sequencing (ddRADseq). We also report the identification of 27 significant growth-related QTLs and 17 corresponding genes in this economically important fish.

Sequencing and Genotyping
In total, 8.2ˆ10 8 raw 90 bp reads were generated and 7.3ˆ10 8 clean reads were retained after removal of low-quality raw reads. These high-quality reads were partitioned into ddRAD tags. The average number of ddRAD tags in each individual was 10,428,571. For the analysis of the F 1 mapping population, a total of 264,072 candidate ddRAD loci were inferred from 68 individuals. If both parents were homozygous, the SNPs were eliminated. The SNPs that were significant segregation distortion based on Mendelism by Chi-square analysis were also eliminated. Ultimately, 26,661 SNPs were obtained (see detailed information in Table S1).

Genetic Linkage Map
After construction of the genetic linkage map using the 26,661 SNPs, we obtained 24 linkage groups (LGs), which is consistent with the haploid chromosome number of the orange-spotted grouper [15]. A total of 3029 SNPs were successfully mapped to the 24 LGs (Table S2 and Figure S1) after discarding contradictory SNPs. These LGs spanned a total genetic distance of 1231.  Table S3.

Discussion
In this study, we applied the technology of ddRADseq, an extension of RADseq [16]. This technique has the advantage of providing improved efficiency and robustness by utilizing two restriction enzymes (a frequently-cutting enzyme and a rare-cutting enzyme) [16][17][18][19]. Under-sampling in read count was a constant question resulting from biased read representation in pooled sequencing experiments among individual samples. ddRADseq increases its sturdiness compared to RADseq [20][21][22]. ddRAD-based genetic maps of Nile tilapia (Oreochromis niloticus) [23], Japanese eel (Anguilla japonica) [17], Midas cichlids (Amphilophus spp.) [24,25], and Eurasian perch (Perca fluviatilis L.) [19] have been successfully constructed. Using ddRADseq, 264,072 candidate RAD loci were inferred, and 3029 high-quality SNPs were retained after a series of filtering. The number of mapped genetic markers in the present study was less than that in our previous report [14]. However, ddRAD-seq was better than MSG (multiplexed shotgun genotyping) in the detection of SNPs because MSG with a stringent methodology may overlook some significant loci [26]. There were fewer markers mainly due to the use of a small F 1 full-sib population.
QTL mapping is considered an efficient strategy for analyzing complex quantitative traits in a variety of fish species. For example, 17 QTLs were detected for the traits of body size (body length and weight) of Chinook salmon (Oncorhynchus tshawytscha) [27] and 11 QTLs were found to be related to the body shape of lake cichlid fishes [25]. Additionally, growth-related QTLs had been identified in Atlantic salmon (Salmo salar) [5] and Arctic charr (Saivelinus alpinus) [28,29]. Recently, a genetic linkage map of kelp grouper had been constructed using simple sequence repeat (SSR) markers, and growth-related (body weight and total length) QTL analysis was performed [30]. The high-resolution genetic map consisted of 714 SSR markers. One major growth-related QTL and several putative QTLs were detected.
The previous data were from the family without phenotype information and were used to construct a SNP-based high-density genetic linkage map. However, our current mapping data from ddRADseq were applied to construct a genetic map and perform QTL analysis based on the available phenotypic information. Twenty-seven QTLs associated with growth (body weight or body length) were identified and found to be distributed on six LGs (LG1, LG2, LG5, LG7, LG21 and LG24). Interestingly, 21 out of the 27 QTLs were concentrated on LG5 within a narrow region (15.1-33.3 cM). These had the highest LOD value (4.7), accounting for 14.5%-29.1% of the phenotypic variation. The small physical and genetic distance among QTLs within the cluster suggested that the cluster would be highly effective for future marker-assisted selection.
Target genes associated with the growth traits of orange-spotted grouper have been previously studied [31]. We compared our QTLs with the scaffold assembly and annotation of the orange-spotted grouper reference genome [32], and discovered a total of 17 genes within the QTL regions. Three of these genes (notch1, ftz-f1 and tert) have been reported to be involved in fish growth [33][34][35]. The notch1 signaling gene plays a role in the notochord development of zebrafish [32]. The development of the notochord is strongly linked with body length. Growth was correlated with gonadal development and sex change in groupers [36]. In tilapia, ftz-f1 is involved in the development of adrenal-gonadal and sex determination and its transcripts were only expressed in the gonads and kidneys [34]. Phylogenetic analysis of the tilapia ftz-f1 indicated that it was highly conserved among other teleosts. In addition, the orange-spotted grouper is a protogynous hermaphroditic and the expression of ftz-f1 in gonads is influenced by the development of the gonad [37]. The gene of tert is associated with the physiological aging of teleosts and plays a major part in the cell process of proliferation, differentiation and tumorigenesis [35]. Expression of tert was also significantly correlated with muscle telomerase activity (TA) in the skeletal muscle of many fish species [38]. TA has a major effect on the growth rate of orange-spotted grouper. Fourteen additional genes (Table S4) were also mapped to growth-related QTLs, but their exact roles with regard to fish growth remain to be determined.

Sample Preparation
Parent fishes (orange-spotted grouper) were captured from the South China Sea in Hainan. One female and one male, with desirable properties such as the coefficient of mature (stage IV) gonad, genetically diverse and highly heterozygous, high vigor, etc., were selected. An F 1 full-sib family was generated by crossing the female and male fish at the Daya Bay Seawater Fish Farm in Huizhou, Guangdong Province, China. The F 1 progeny were raised under a natural photoperiod, at the water temperature 28˝C. They were fed according to management practices of the fishery. Fin clips of the parents and 68 offspring (four months of age) were collected and stored in absolute ethanol at 20˝C before use. Phenotypic information on the F 1 progeny is summarized in Table S5. Phenotypic correlations were tested with the Pearson correlation coefficient, and the results indicated a significant correlation (p < 0.01) between body weight and body length (Table S6). Genomic DNA was extracted using a standard phenol-chloroform protocol [39]. The quality of all DNA samples was evaluated by Qubit Fluorometer (Invitron, Waltham, MA, USA). Electrophoresis was conducted on 0.6% agarose gels. All the experiments were carried out in accordance with the guidelines of the Animal Ethics Committee and were approved by the Institutional Review Board on Bioethics and Biosafety of BGI (No. FT14015, 12 March 2014).

ddRAD Library Construction and Sequencing
The ddRADseq library was prepared using a previously reported protocol [16] with some modifications. Briefly, the workflow included double-digestion, ligation reaction, pooling, purification, and amplification. The double-digestion and ligation reactions were prepared as 30-and 40-µL reactions (see more details in Tables S7 and S8), respectively. One µL of the ligation products from each sample was pooled 24 samples in a new centrifuge tube, respectively. The DNA fragments were purified (300-700 bp) were isolated on a 3% agarose gel, and then extracted with the QIA quick Gel Extraction Kit (Qiagen, Hilden, Germany). All the samples were amplified by 12 cycles of PCR in the 50-µL reaction volumes containing 2-µL purified of the PCR products, 20-µL Master Mix, 2-µL primers (1-µL each) and 26-µL ddH 2 O. The amplified products were purified using the QIA quick PCR Purification Kit (Qiagen) and checked the quality (concentration and purity) of DNA products via an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Finally, library sequencing was performed using a Hiseq 2000 platform (Illumina, San Diego, CA, USA) with 90-bp pair-end reads.

SNP (Single-Nucleotide Polymorphism) Calling and Genotyping
We discarded the Illumina short reads without sample-specific barcodes and restriction enzyme motifs. Subsequently, we used Soapsnp [40] to sort retained reads into loci and to genotype them. Lastly, we aligned the clean reads of the two parents to the previously assembled orange-spotted grouper genome [32] to eliminate monomorphic DNA sequences and then added the reads containing SNPs. SNPs were selected as described previously [41]. Regions with these putative SNPs were defined as reference SNP regions, and then the clean reads of progeny samples were aligned to the reference SNP regions. As noted earlier, the parental genotypes determined the genotypes of the individual offspring [42].

High-Density Genetic Map Construction
The SNPs that at least one parent was heterozygous and had high quality of genotype calls were retained for genetic linkage analysis. Those SNPs with significant segregation distortion (p < 0.01, x 2 test) were discarded. We used JoinMap 4.1 (Kyazma B.V., Wageningen, Gelderland, Netherlands) to conduct the genetic map construction with the genotypic data from the F 1 mapping population. The logarithm of odds (LOD) threshold was set as from two to 15. It was adopted as indicator to cluster analysis using regression mapping algorithm. Map distances were calculated using the Kosambi mapping function in centiMorgans (cM). All the SNPs were arranged and grouped into 24 linkage groups (LGs). Graphics of the linkage groups were generated using custom perl script.

QTL Mapping
We used WinQTLCart2.5 software (Wang S., Basten C.J. and Zeng Z.-B., Raleigh N.C., USA) to perform the QTL analyses with the method of composite interval mapping (CIM) [43]. The CIM method was executed using Model 6 with parameters of five control markers, a 10-cM window size and backward regression. We used 1000 permutations and the p-value (0.05) of whole genome-wide significance to determine the threshold level of likelihood ratio (LR) (11.5) or LOD (2.5). The LR peak location and surrounding district of every QTL was determined. Phenotypic variation resulting from growth-related QTLs was also calculated by the software.

Identification of the Genes on Growth Traits Associated QTLs
The genes on growth-associated QTLs were identified via following steps. First we identified the orange-spotted grouper genome scaffolds [32] from these QTL regions. Subsequently, we mapped the corresponding gene identities (IDs) based on their position on the scaffolds. Finally, we retrieved the homologous genes from the orange-spotted grouper gene annotation file [32].

Conclusions
A high-density genetic linkage map with a total genetic distance of 1231.98 cM was constructed for the orange-spotted grouper using the ddRADseq method. Twenty-seven growth-related QTLs were identified by QTL mapping and corresponding 17 genes were discovered. These results will provide a valuable resource in future marker-assisted molecular breeding studies of this economically important fish species.