Linkage mapping and QTL analysis of growth traits in Rhopilema esculentum

R. esculentum is a popular seafood in Asian countries and an economic marine fishery resource in China. However, the genetic linkage map and growth-related molecular markers are still lacking, hindering marker assisted selection (MAS) for genetic improvement of R. esculentum. Therefore, we firstly used 2b-restriction site-associated DNA (2b-RAD) method to sequence 152 R. esculentum specimens and obtained 9100 single nucleotide polymorphism (SNP) markers. A 1456.34 cM linkage map was constructed using 2508 SNP markers with an average interval of 0.58 cM. Then, six quantitative trait loci (QTLs) for umbrella diameter and body weight were detected by QTL analysis based on the new linkage map. The six QTLs are located on four linkage groups (LGs), LG4, LG13, LG14 and LG15, explaining 9.4% to 13.4% of the phenotypic variation. Finally, 27 candidate genes in QTLs regions of LG 14 and 15 were found associated with growth and one gene named RE13670 (sushi, von Willebrand factor type A, EGF and pentraxin domain-containing protein 1-like) may play an important role in controlling the growth of R. esculentum. This study provides valuable information for investigating the growth mechanism and MAS breeding in R. esculentum.

Edible jellyfish R. esculentum distributes in the northwest Pacific Ocean and is a popular seafood in Asian countries, especially in China 1,2 . R. esculentum is rich in protein and minerals while low in calories and fats, making it an ideal nutritive ingredient for developing oral formulations, functional food and nutricosmetics 3 . In addition, the collagen peptides from R. esculentum can accelerate the wound healing process of mice 4 and have antihypertensive activity 5 , suggesting it can be applied in pharmaceutical industry. The edibility, nutritional value and medicinal properties of R. esculentum make it an economical fish resource and widely farmed in aquaculture systems of China 2 . However, overfishing decreased the population number of R. esculentum and caused the scarcity of natural resource 6 . To sustainable development and utilization of R. esculentum resources, it's necessary to investigate the genetic background and genetic improvement.
Genetic linkage maps play an important role in studies of genome and genetic 7 . The advance of the 2b-RAD method attracts many researchers' attention and has accelerated the identification of SNP markers for constructing genetic linkage maps in farmed fishes, such as Carassius auratus 8 , Cyprinus carpio haematopterus 9 , Hypophthalmichthys nobilis 10 and Hemibagrus wyckioides 11 . The average marker interval of these linkage maps was between 0.44 and 0.87 cM and helped mapping QTLs of interested growth traits for genetic breeding [8][9][10][11] . Based on the constructed linkage maps, many QTLs related to economic traits, such as body weight, body length, nutritional metabolisms and sex have been identified 7,10-12 . For example, one QTL related to body weight in Apostichopus japonicas, explaining 11.8% of the phenotypic variation 7 .
Although the R. esculentum genome has been released 1,13 , the genetic linkage map and QTL for growth in R. esculentum have not been reported yet. In previous studies, the researchers only identified some markers in R. esculentum, such as microsatellite for detection of genetic diversity and conservation of germplasm resources 14 and SNPs and simple sequence repeats (SSRs) for assisting MAS breeding 15 . To improve the growth of R. esculentum through MAS breeding, we identified SNP markers by 2b-RAD method and constructed the first genetic linkage map of R. esculentum. Furthermore, we detected some QTLs and genes related to growth traits of R. esculentum based on the linkage map. This study will provide some insights for investigating the growth mechanism and genetic research of R. esculentum in the future.

Results
2b-RAD sequencing. A total of 1539 million raw reads were generated by 2b-RAD sequencing, including 19 million from the male and female parents and 1520 million from the offspring ( Table 1). The ratio of clean reads to raw reads is higher than 93.63% and the alignment rate of the clean reads to the genome reference sequence of R. esculentum was close to 63% ( Table 1). The tag number of female and male parents was both higher than that of F1 offspring while the tag depth was both lower than that of F1 offspring.

SNP markers and linkage mapping.
A total of 9100 SNP markers were generated with the transition to transversion (TS/TV) ratio of 1.48. After filtering, 6674 SNP markers can be used for constructing the linkage map. Through the neighbor-joining tree and principal component analysis, we excluded the abnormal individuals and confirmed the parents and offspring as a whole for constructing the linkage map ( Supplementary Figures). For female and male-specific linkage maps, 1427 and 1460 SNP markers were separately selected ( Table 2).
The SNP markers were grouped to 21 LGs, corresponding to the chromosome number of R. esculentu 13 (Fig. 1). The genetic lengths of the female and male-specific maps were 1360.29 and 1200.62 cM with an average marker interval of 0.95 cM and 0.82 cM, respectively (Fig. 1).
The consensus linkage map of R. esculentum was constructed using 2508 SNP markers with map coverage of 98.68% (  Table S1). The longest LG is LG14 with a genetic length of 88.9 cM, 1.75-fold higher than that of the shortest LG20 (Fig. 2). QTL mapping. Six QTLs were detected for the growth of R. esculentum, including three for umbrella diameter and three for body weight (Table 3). These six QTLs located on four LGs, explaining 9.4 to 13.0% variation ( Table 2). Two QTL regions for umbrella diameter that located on LG 14 and 15 are close to the QTL regions for body weight. In addition, the LOD value of the identified QTLs floated between 3.0 and 4.3 (Table 3).
Candidate genes related to growth traits. A total of 28 and 35 candidate genes in the QTL regions were identified associated with umbrella diameter and body weight, respectively (data not shown). Of these candidate genes, 27 genes were overlapped, associating with both umbrella diameter and body weight ( Table 4). As we can see from Table 4, 22 genes have functional annotation and ten genes have metabolic roles, such as thermogenesis, mTOR signaling pathway, glutathione metabolism, biosynthesis of nucleotide sugars according to the KEGG database (Table 4). Due to the calcium-binding EGF-like (EGF_CA) domain, RE13670 (sushi, von Willebrand factor type A, EGF and pentraxin domain-containing protein 1-like) was considered to play a key role in controlling the growth of R. esculentum (Table 4). Using online software NCBI-blast, we found RE13670 was most homologous to sushi, von Willebrand factor type A, EGF and pentraxin domain-containing protein 1-like of Acropora digitifera with a percent identity of 30.36%, supporting the accuracy of gene annotation. However, the metabolic roles of RE13670 were not found (Table 4).

Discussion
The molecular marker shows potential for investigating the growth of aquaculture animals in the genetic breeding industry 16 . SNP markers as one of the molecular markers were genotyped easily to construct the genetic linkage maps for guiding the genetic breeding of aquaculture animals 17 . Researchers have identified 1,034,708 SNPs in R. esculentum by transcriptome sequencing 15 , yet we only identified 9100 SNPs by 2b-RAD sequencing, much lower than that in the previous study 15 . The difference in SNP numbers between the two studies may be caused by the SNP analysis method. In the previous study, SNPs were detected by GATK2 software without reference genome while by RADtyping software according to R. esculentum genome in this study. Due to the simplicity and  www.nature.com/scientificreports/    www.nature.com/scientificreports/ identified SNPs in R. esculentum are very important markers for genetic breeding and this is the first report of linkage map in R. esculentum. The high-density linkage map plays an important role in performing QTL mapping and finding genes related to the growth traits in aquatic animals 7,8,10,11 . Numerous QTLs about growth traits, such as body weight, body length, sex as well as disease resistance were identified based on the high-density linkage map 8,[19][20][21][22][23][24][25] . For Nibea albiflora, 15 QTLs were detected associated with body weight, explaining 14.7-35.7% of the phenotypic variations 21 . For C. auratus at 2 months, eight QTLs in eight chromosomes were discovered associated with the body weight, explaining 10.1-13.2% of the phenotypic variations 8

. In this study, three QTLs distributed in three
LGs were detected associated with body weight of R. esculentum and explained 9.4-10.5% of the phenotypic variations, following the result of C. auratus 8 . In addition, our studies showed that growth-related traits body weight and umbrella diameter in R. esculentum are positively correlated with the two close QTLs (Supplementary  Table S2; Table 3).
Positional cloning of candidate genes with the help of QTL mapping may provide an efficient method for selective breeding in aquaculture animals 8,21,26 . We identified 27 candidate genes in LG 14 and 15 corresponding to the growth of R. esculentum and RE13670 showed the most possibility in controlling the growth due to the EGF_CA domain, following the growth-related genes reported in C. auratus 8 . For C. auratus, five candidate genes show potential for body weight and two genes (TGF-beta and EGF-like domain) may be the most promising according to their role in early growth and development of vertebrates 8 . EGF_CA domain is a calcium-binding EGF-like domain and needs calcium for performing biological function 27 . EGF_CA domain presented in extracellular (mostly animal) and membrane-bound 27 and has three main roles, including protein-protein interactions, as a spacer unit and structural stabilization 28 . Although the functional significance of EGF_CA domain in aquaculture animals is unclear, it is worth confirming if EGF_CA domain is associated with the growth of R. esculentum. With the release of genomic data of R. esculentum 1,13 and the development of biotechnology, the gene function studies of R. esculentum will be improved for studying genetic breeding.

Conclusion
In this study, a total of 9100 SNP were identified and a high-density linkage map was constructed with a marker interval of 0.58 mM in of R. esculentum using the 2b-RAD method. Based on the linkage map, six QTLs were identified associated with the growth of R. esculentum and one candidate gene RE13670 containing EGF_CA domain in LG14 may play the key role in controlling the growth of R. esculentum. Although one full-sib family of R. esculentum is limited, the identified SNPs and genes for growth will accelerate the MAS breeding of R. esculentum.

Materials and methods
Mapping family and DNA extraction. A full-sib family of R. esculentum was constructed in Yingkou City, Liaoning province, China. We selected the two parents with a large difference in growth from the breeding pond and cultured their offspring in the nursery pond. The two parents and random 150 offspring at sevenmonth-old (juvenile jellyfish) were chosen for sequencing. The body weight and umbrella diameter of the offspring were measured (Supplementary Table S1). The genomic DNA was extracted referred to the previous method 13 . The DNA quality was measured by Qingdao OE Biotech Co., Ltd.
2b-RAD sequencing. The 2b-RAD libraries of R. esculentum specimens were constructed at Qingdao OE Biotech Co., Ltd., following the published method 29 . Firstly, 100-200 ng genomic DNA was digested by 1 U BsaXI (New England Biolabs). Secondly, the ligation reaction was conducted to add specific adaptors to the digested genomic DNA. Thirdly, the ligation products were amplified in MyCycler thermal cyclers (Bio-Rad). Fourthly, the PCR products were purified using a MinElute PCR Purification Kit and digested using SapI (New England Biolabs). Fifthly, the digested products were transferred to the tube containing magnetic beads for incubation and then transferred the supernatant to a new tube for ligation using T4 DNA ligase (New England Biolabs). After that, the ligation products were purified and barcodes were introduced by PCR using barcodebearing primers. Finally, PCR products were purified and pooled for sequencing using the Illumina Novaseq 6000 PE150 sequencing platform.
The raw data of 2b-RAD sequencing were trimmed for getting the high-quality data, and then the high-quality data were aligned to reads. The reads with the BsaXI site were extracted and aligned to the reference reads using SOAP (version 2.21) 30 . The reference reads were extracted from R. esculentum genome (NCBI Genome ID: 56778) after electronic digestion using the BsaXI enzyme.
SNP detection and filter. The aligned data were used to detect SNP by RADtyping software 31 . The maximum-likelihood (ML) algorithm was used to detect homozygote or heterozygote in co-dominant markers. SNPs were filtered with the following criteria: (1) Segregating markers that could be genotyped over 80% of the progenies were retained; (2) Markers with a minor allele frequency (MAF) less than 0.01 were discarded; (3) Polymorphic loci that contain more than two alleles were excluded; (4) The aligned reads with more than two SNPs were discarded. Based on the detected SNPs, we constructed neighbor-joining tree using treebest (version 4.1) 32 and performed principal component analysis using ADMIXTURE (version 1.3.0) 33 for confirming the species as a full-sibling family. SNP markers were annotated using the software SnpEff (version 4.1) 34 .
Linkage map construction. After filtering, the linkage group was divided with the LOD value of 2-15.
Marker distances were calculated using Kosambi's mapping function 35 . The male and female-specific linkage maps were constructed by software JoinMap (v5.0) 36 . The male-specific linkage map was constructed using