High-resolution linkage and quantitative trait locus mapping using an interspecifc cross between Argopecten irradians irradians (♀) and A. purpuratus (♂)

Scallops represent a conspicuous group of bivalve molluscs, comprising over 300 extant species worldwide (Brand 2006). They are valuable models for studying early animal evolution and bivalve adaptation as they possess remarkable preservation of ancestral bilaterian linkage groups, and novel genomic features for adaptations to benthic life as semi-sessile filter feeders (Li et al. 2017; Wang et al. 2017). Also, scallops are popular aquaculture species and provide significant economic benefits. The production of maricultural bivalves is more than 13 million tonnes per year from 2010 to 2015 (Wijsman et al. 2019). The bay scallop, Argopecten irradians irradians, and the Peruvian scallop, A. purpuratus, are economically important species naturally distributed along the Atlantic coasts from Nova Scotia to the Gulf of Mexico and Pacific coasts from Peru to Chile, respectively (Wolff and Mendo 2000). Geographic isolation and climatic variability promoted their divergence and allopatric speciation (Brand 2006; Marín et al. 2013). Despite that, Argopecten scallops have distinct size, life span, temperature and hypoxia tolerance, these two species share similar life history and karyotypes (Hu et al. 2013, 2015; Wang and Guo 2004). Interspecific hybridization between the two species exhibited extraordinary increase in growth rates in both shell height and whole-body weight in F₁ offspring (Wang et al. 2011). The exploitation of this heterosis could greatly benefit the scallop industry (Hu et al. 2013).

Genetic linkage mapping is an efficient method to investigate the molecular mechanisms of hybridization and introgression between species. On the basis of genetic mapping, suppression of gene flow could be revealed by limited recombination between rearranged chromosomes (Ostberg et al. 2013; Panithanarak et al. 2004). Selections against particular hybrid genotypes result in serious deviations of observed genotypic frequencies from the expected Mendelian inheritance (Fishman et al. 2001; Kim et al. 2014). Therefore, segregation distortion markers have been extensively studied with interspecific populations to reveal the intrinsic incompatibilities (Brennan et al. 2014; Phadnis and Orr 2009; Streiff et al 2014).

Moreover, genetic mapping serves as a chromosomal-level framework to detect quantitative trait loci (QTL) including that for heterosis and trace inheritance of various performance traits (Garcia et al. 2008; Shang et al. 2016). The dominant, over-dominant, epistatic and other effects for heterosis have been tested in plants and multiple livestock and aquaculture species (Gardner and Lonnquist 1959; Kaeppler 2012; Krieger et al. 2010; Kusterer et al. 2007; Moll et al. 1964; Schwartz and Laughner 1969; Semel et al. 2006; Wolf and Hallauer 1997; Xiao et al.1995; Yu et al. 1997). Limited by the number of developed markers and high-throughput genotyping platforms, genetic studies in scallops lack progress in uncovering the molecular mechanisms underlying various commercial traits. Most of the previously constructed genetic maps in scallops were with relatively low resolutions, limiting the sensitivity and accuracy in determining quantitative trait genes (Li et al. 2005, 2012; Wang et al. 2004, 2005, 2007; Zhan et al. 2009). Recent development of genotyping by sequencing (GBS) technologies has allowed high-throughput and efficient genotyping (Elshire et al. 2011; Narum et al. 2013; Poland and Rife 2012). Of these techniques, an optimized representative restriction site-associated DNA (RAD) method, 2b-RAD, has gained popularity owing to its simplicity in library preparation and ability to provide even and tunable genome coverage with a low requirement of sample quality (Wang et al. 2012, 2016). This method has been successfully applied for genetic and breeding studies, such as high-density linkage map construction and QTL mapping, in a wide range of species (Barfield et al. 2016; Cui et al. 2015; Dixon et al. 2015; Guo et al. 2014; Jiao et al. 2013; Lowry et al. 2015; Pauletto et al. 2016; Seetharam and Stuart 2013; Tian et al. 2015).

In the present study, we conducted genetic linkage analysis using the 2b-RAD genotyping method based on a F₁ interspecific family with the objectives of (1) construction of high-density genetic maps for the two scallops, (2) identification of segregation distortion loci to provide insights into genomic incompatibilities of the two species, and (3) genome scan of QTLs for growth-related traits to dissect the genetic architecture of interspecific heterosis.

The sequencing data of the mapping family including two parents and 150 progenies (Fig. 1) are summarized in Supplementary Table S1. A total of 656 million reads were generated for the whole family, and 86.7% of the total reads passed quality trimming. For the parents, standard libraries were constructed and sequenced to a sufficient depth. A total of 33.1 and 32.3 million high-quality reads were generated for the parental bay scallop and Peruvian scallop, corresponding to the sequencing depth of 128.2 × and 122.8 ×, separately. For the hybrid progenies, RR libraries were used to reduce the total sequencing cost. Approximately, 3.4 million high-quality reads for each progeny were generated with an average sequencing depth of 22.3 ×, which was sufficient for accurate genotyping (Fu et al. 2013).

Figure 1.  The interspecific cross between Argopecten irradians irradians (♀) and A. purpuratus (♂), and two shell color types of the hybrid progenies. Type Ⅰ: The left valves were orange and the right were purplish orange. Type Ⅱ: The left valves were purplish white with different degree coverage of dark blotches and the right were purplish white

Clustering of the high-quality reads of the two parents generated 199, 962 representative reference sites, including 86, 544 parent-shared (~ 43.3%) and 113, 418 parent-specific (~ 56.4%) sites. After filtering sites with insufficient depth or from repetitive genomic regions, 73, 027 parent-shared and 111, 696 parent-specific sites were retained. For the parent-specific sites, 71, 692 (~ 64.2%) sites were from the bay scallop, and the remaining 40, 004 (~ 35.8%) sites were from the Peruvian scallop. Sequencing reads were then remapped to the constructed reference for genotyping. In total, 7653 polymorphic markers were identified in at least 80% progenies, which consisted of 1696 codominant markers and the rest 5957 dominant markers. The genotyping results for each segregation type are summarized in Table 1. The χ² test was carried out to examine segregation distortion of the codominant markers and identified 924 markers that significantly deviated from the expected Mendelian ratio (P < 0.05), including 508 bay scallop-specific SNPs (Aa × aa type), 363 Peruvian scallop-specific SNPs (aa × Aa type), and 53 shared SNPs (Aa × Aa type) (Supplementary Figs. 1, 2). Notably, 93% of the shared codominant markers were identified with non-Mendelian inheritance.

Figure 2.  The high-resolution genetic linkage map for (a) A. irradians irradians containing 2994 markers in 16 linkage groups, and (b) A. purpuratus containing 1585 markers in 16 linkage groups. The white spaces are the gaps without markers mapped and black lines represent the positions of markers in the linkage maps

In total, 6729 makers in accordance with Mendelian segregation (P ≥ 0.05) were used for the linkage mapping. These markers were partitioned into bay scallop heterozygous (Aa × aa and Aa × Aa type; 4561 markers) and Peruvian scallop heterozygous (aa × Aa and Aa × Aa type; 2172 markers) for constructing species-specific linkage maps. Using a LOD threshold of 10.0, markers were assigned into 16 linkage groups (LGs), which was consistent with the number of chromosomes of Argopecten haploid genome (Hu et al. 2013).

The linkage map of the bay scallop contained 2994 markers and spanned 945.95 cM (Table 2, Fig. 2a). The average marker interval was 0.32 cM, ranging from 0.18 cM in LG1 to 0.63 cM in LG16. The expected map length was estimated to be 957.55 cM (G_e1) and 959.04 cM (G_e2) by two different methods. Genome coverage calculated using the average of the two estimates (958.30 cM) was as high as 98.9%. For the Peruvian scallop map, 1585 markers were mapped and spanned over 800.46 cM (Table 3, Fig. 2b). The average marker interval was 0.51 cM, ranging from 0.29 cM in LG1 to 1.35 cM in LG14. The constructed genome was estimated to be 819.01 cM in length by averaging G_e1 (816.78 cM) and G_e2 (821.24 cM) and covering about 97.7% of the whole genome. A few marker intervals were notably large in the genetic maps, which may represent recombination hotspots. Due to the large genetic differences between the bay and Peruvian scallops, insufficient shared markers were eligible to integrate the two maps. Therefore, consensus maps of the two species were not constructed.

The Kolmogorov–Smirnov test indicated that all the five growth-related traits followed normal distribution. The Pearson's correlation coefficients among the five growth traits ranged from 0.69 (shell width with shell weight) to 0.93 (shell height with shell length), indicating that they were highly correlated and possibly controlled by the same set of QTLs (Table 4, Supplementary Fig. S3). Therefore, we conducted a principle component analysis on the five traits. The first principle component (PC1), which explained over 80% of the total variation, was used as a composite trait for QTL mapping (Fig. 3). The full stepwise multiple QTL model identified a sole growth-related QTL (qgt) in LG2 (1.3–1.9 cM) of the Peruvian scallop, generated a LOD score of 4.7 and explained 13.1% of the variation in PC1 of the five traits (Fig. 4a). Evidence for epistatic effects has been found also among markers on LG1, LG2, and LG8, despite that the interactions do not meet the model selection significance threshold (Fig. 4b). Candidate genes, including SPR (sex peptide receptor), Trip13 (pachytene checkpoint protein 2 homolog), and KIF13B (kinesin-like protein KIF13B), were found to be located within the 95% credible interval of the QTL detected on LG2. These genes could be of causal importance that is worthy of further investigation. Notably, no LOD peaks above the threshold were detected in the bay scallop, which may be due to the non-polymorphism of candidate QTLs or associated markers in the bay scallop (Fig. 4c). Genomic scan of LOD profiles for five individual traits were also performed and revealed similar results (Supplementary Figs. S4, S5).

Figure 3.  Principal component analysis of five growth traits in the hybrid scallops. a The PC1 explained about 85% of the total variance. b SW shell width, BW body weight, AMW adductor muscle weight, SL shell length, SH shell height

Figure 4.  QTL mapping of growth traits and shell color in the hybrid scallops. The PC1 for five growth traits was used in QTL mapping with the linkage map of Peruvian scallop (a) and bay scallop (c). b LOD scores from a two-dimensional scan for growth traits. LOD_i measures the improvement in the fit of the full model over that of the additive model. It is displayed in the upper left triangle. LOD_f measures the improvement in the fit of the full two-locus model over the null model, indicating the evidence for at least one QTL, with allowance for interaction. It is displayed in the lower right triangle. QTL mapping for shell color with the linkage map of bay scallop (d) and Peruvian scallop (e). The solid and dashed lines indicate the LOD significance thresholds

The two shell color types were checked to be consistent with the segregated ratio of 1:1 (P ≤ 0.05). Despite that, a single significant QTL (qsc) on LG2 in the bay scallop was identified with a LOD value of 45.1, no significant epistatic interactions were detected by the stepwise QTL selection model (Fig. 4d). The Bayesian 95% credible interval of qsc was only 0.83 cM, ranging from 13.2 to 14.1 cM. It explained ~ 74.9% of the variation in shell color. By mapping the markers to the bay scallop genome assembly, potential causal genes were identified, including PECR (peroxisomal trans-2-enoyl-CoA reductase), MECOM (MDS1 and EVI1 complex locus), and RTTN (Rotatin), while for the Peruvian scallop, no significant QTL signal was identified (Fig. 4e).

Scallops are important fishery and aquaculture species, with a production of over 2 million metric tonnes in 2016, accounting for 13% of the global mollusc production (FAO STAT 2016). Much of the research effort has been on evaluation and selection of crossbreeding and interspecific hybridization (Hu et al. 2015; Lu et al. 2006; Wang et al. 2010, 2011; Zhang et al. 2014). Understanding the genetic bases underlying genomic incompatibility and heterosis has been an important task to facilitate scallop breeding and genetic improvement. In this study, we performed genetic mapping with an interspecific hybrid resource family. Taking advantage of 2b-RAD method, high-resolution genetic linkage maps were constructed for the bay and Peruvian scallops by cost-effectively genotyping a full-sibling family of 152 scallops. The mapping density of the bay scallop was 0.32 cM with a 97.7% coverage, surpassing the resolution of all the previously constructed genetic maps (Li et al. 2012; Qin et al. 2007a, b; Wang et al. 2007). The linkage map for the Peruvian scallops constructed in this study is the first one for this species, with a resolution of 0.51 cM and a mapping coverage of 98.9%. Large differences in the polymorphic sites have been observed between the two species, with ~ 60.5% markers were species-specific polymorphic. A few marker intervals (i.e. 25.9 cM in LG16) were notably large in the current linkage map, which may represent hotspots of recombination. Therefore, more markers in these regions could be used for further study. Only ~ 3.4% codominant markers were shared between the two species, resulting in insufficient markers for construction of consensus maps.

Markers deviating significantly from Mendelian inheritance were excluded from linkage mapping to avoid miscalculation. Despite that the level of segregation distortion varies greatly among species, ~ 54.5% markers showing segregation were identified in the hybrid progenies, which is much higher than the general segregation distortion ratio in marine bivalves, such as the Pacific oyster (Li and Guo 2004), the Zhikong scallop (Zhan et al. 2009), and the bay scallop (Li et al. 2012). In this study, over 93.0% of the SNPs shared by the two parents (i.e., Aa × Aa) deviated from the Mendelian segregation. Consistent with previous studies, the distorted markers tend to locate adjacently to form segregation distortion regions (SDRs), implying that inter-chromosomal rearrangements may occur within these regions (Douglas and Par 2003; Liu et al. 2010; Lyttle 1991). Chromosomal rearrangements may influence inheritance patterns by suppressing recombination in heterokaryotypes (Sturtevant 1917). Studies have reported that rearrangements, such as inversions, represent sets of coadapted alleles in a population that live in a constant environment. When two distinct populations are connected by migration, recombination loci involved in local adaptation regions are typically selected against (Dobzhansky 1947; Ortiz-Barrientos et al. 2002).

Elucidation of the genetic basis underlying performance traits of scallop has been a central task for genetic improvement. QTLs for growth- and reproduction-related traits have been detected in the bay scallop (Li et al. 2012), the Zhikong scallop (Jiao et al. 2013), and the Pacific lion-paw scallop (Petersen et al. 2012). As revealed by Wang et al. (2011), the hybrids exhibited higher growth rates in both shell height and whole-body weight than the parental scallops, with a mid-parental heterosis of 125.9%–138.9% for whole-body weight and 145.4%–156.2% for adductor muscle weight. In this study, growth-related traits were highly correlated (Table 4), suggesting that they were possibly regulated by a same set of QTLs. Therefore, a principle component analysis was conducted to generate a composite trait for QTL mapping. By applying the newly constructed high-resolution linkage maps, a major QTL for the composite growth traits was mapped to LG2 of the Peruvian scallop. Despite that the gaps in genomic assembly may limit identification of the potential causal genes within the QTL, growth-related genes, including SPR, Trip13, and KIF13B were recognized within the loose window of this QTL (Supplementary Table S2). These genes are involved in the control of reproductive behaviors (Bath et al. 2017; Yapici et al. 2008), gamete generation (Li and Schimenti 2007), and cytokinesis (Kurasawa et al. 2018; Venkateswarlu et al. 2005). The list of candidate genes and SNPs associated with traits supplies potential targets in the future selective breeding programs to improve scallop aquaculture production, but more validation work is still needed in a broader population.

In molluscs, shell coloration plays important roles in crypsis, mimicry and aposematism against visual predators including crabs, fish, and birds (reviewed in Williams 2017). Variations in shell colors and pigmentation patterns also associate with differences in growth, survival, and biological adaptations (Gary 1980; Wolff and Garrido 1991; Zheng et al. 2005). The mechanisms of shell coloration are complicated; they are affected both by environmental and genetic factors. Environmental effects including lights, climate, salinity and diets on shell color have been reported in previous studies (Lindberg and Pearse 1990; Liu et al. 2009; Miura et al. 2007). These modifications could be stable and long lasting once formed (Liu et al. 2009). Also, studies have revealed that shell colors are subject to strict genetic control. The general colorations and banding patterns of periostracum layer, raylike marks, and prismatic layer are controlled by single locus variations in the blue mussel (Newkirk 1980), hard clam (Chanley 1961), bay scallop (Adamkewicz and Castagna 1988), pearl oyster (Wada and Komaru 1990), and Yesso scallops (Zhao et al. 2017). Other coloration characteristics, such as pigmentation patterns, appear to be controlled by the epistatic effects of multiple nonallelic genes (Adamkewicz and Castagna 1988; Winkler et al. 2001). In this study, the hybrid scallop shells exhibited distinct background colors with a segregation ration of 1:1. A significant QTL signal on LG2 was detected in the bay scallop linkage map, explaining over 74.9% of the phenotypic variation. The markers located at the confidence intervals of this QTL constitute a group of valuable markers for further evaluation of their utility in marker-assisted selection. Candidate causal genes, and in particular, the MECOM gene, were identified within the QTL and focused our interest (Supplementary Table S3). MECOM is a transcriptional factor that is involved in regulation of myeloid differentiation of bone marrow progenitor cells to erythroid cells and granulocytes (Buonamici et al. 2003). Heme biosynthesis is induced during erythroid differentiation and is the main pathway for the formation of uroporphyrin pigments (Layer et al. 2010). Studies have revealed that uroporphyrin pigments are responsible for coloration of shells and soft tissues in several molluscs (reviewed in Williams 2017). In scallops and marine snails, several genes related to porphyrin metabolism were differentially expressed in mantles with different colors (Ding et al. 2015; Williams et al. 2017), whereas for MECOM, the difference is not significant in adult mantle tissues suggesting that it may express and function during early shell formation. In this study, the QTLs related to shell coloration and growth were revealed by mapping with one bi-parental family, which represents only a small cross section of the breeding germplasm. The expression of QTL depends on its presence and could be influenced by the genetic background (Blanc et al. 2006; Ogut et al. 2015). Therefore, further work would be needed to evaluate the QTL positions and effects with independent populations.

In conclusion, we generated the first high-resolution linkage maps for A. purpuratus and A. irradians irradians. The marker density and mapping coverage were comparable to the highest level of linkage maps in molluscs allowing for comparative genome analysis using both genetic and genomic approaches. QTLs and candidate genes related to growth and shell coloration were detected, indicating important information for further evaluation. These findings are an important addition to the genomic resources for scallop genetic studies and are especially useful for investigations on genomic incompatibility for hybridization, genome evolution of closely related species, and genetic enhancement programs in aquaculture.

The mapping family was constructed using bay scallops from a farmed strain in Jiaonan, Qingdao, China, and Peruvian scallops from the progenies of wild stocks in Peru. As both species are hermaphrodite with external fertilization, inter-specific hybridization of A. irradians irradians (♀) and A. purpuratus (♂) was performed following the procedures previously described by Wang et al. (2011) to avoid self-fertilization. Briefly, sexually mature scallops with gonadal condition at stage Ⅳ (Sastry 1963) were collected and subjected to an air exposure of 30 min followed by a temperature shock from 21 to 23 ℃. Spawning scallops were separated in 5 L capacity polyethylene buckets containing 3 L of 30 μm filtered seawater to prevent contamination. The collected sperm were filtered through a 25 μm mesh screen to mitigate mixed eggs. The eggs were checked under a microscope 20 min post-spawning to verify that they were not self-fertilized. The unfertilized eggs of one bay scallop were mixed with the sperm of one Peruvian scallop for insemination. The mapping family was reared following a previous study (Wang et al. 2011). Briefly, larvae were raised at 21–23 ℃ with aeration, and fed with microalgae, i.e. Isochrysis galbana, at concentrations of 20, 000–100, 000 cells/ml depending on their stages. The density of larvae in the tank was maintained at 10 larvae/ml by adjusting water volume. Collectors were placed into the tanks as substrates for settlement when the eyespots became apparent. After 10 days of metamorphosis, the settled spat were transferred to the field and subsequently sorted into grow-out cages in the open sea of Jiaonan.

After six-months post fertilization, a total of 150 hybrid progenies were randomly taken for genotyping and subsequent analysis. The adductor muscles from the progenies and their parents were collected and stored in 70% ethanol for DNA extraction. Genomic DNA samples were prepared following standard phenol–chloroform DNA extraction procedure (Maniatis et al. 1989). 2b-RAD sequencing libraries were constructed using the protocol developed by Wang et al. (2012). The adaptors and PCR primers were modified to fit the Illumina HiSeq sequencing platforms (Supplementary Table S4). For the parents, standard sequencing libraries were constructed using adaptors with 5′-NNN-3′ overhangs to target all the BsaXI fragments, whereas for the hybrid progenies, reduced representation (RR) libraries were constructed by adaptors with 5′-NNA-3′ and 5′-NNT-3′ overhangs to capture a subset of all the fragments. All the libraries were labeled with unique barcodes and pooled together for single end sequencing (1 × 36 bp) utilizing a HiSeq-2000 system. All the sequencing data are deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) platform with accession no. PRJNA488174.

RADtyping (version 1.0) was used to process the sequencing data and perform de novo dominant and codominant genotyping (Fu et al. 2013). Briefly, raw reads were first trimmed to remove adaptor sequences, low-quality reads (> 5 positions with quality < 10), uninformative reads (with no restriction sites, or containing ambiguous base calls), or reads with excessive homopolymer repeats (> 10 bp). Filtered reads of both parents were assembled into the exactly matching "allele" clusters, and then merged into "locus" clusters allowing two mismatches, which finally compromised the reference sites. Clusters with insufficient coverage (< 8 reads) or excessive depth (> 3000 reads) were excluded to avoid the interference of sequencing errors and repetitive elements. These clusters were further classified as parent-shared or parent-specific sites (sites present in one parent and absent in another due to the mutation of a recognition site) for subsequent dominant and codominant genotyping, respectively. Finally, filtered reads of the 150 progenies were mapped to the constructed reference for genotyping.

Markers were excluded from linkage analysis if they were not polymorphic in any of the parents or had low quality of genotype calls in more than 20% progenies. Mendelian segregation was examined through χ² goodness-of-fit test, markers with significant segregation distortion were also filtered out (P < 0.05). Segregation-distorted regions were identified by aligning these distorted markers to the genome assemblies of the bay scallop (BioProject No. PRJNA428789) and the Peruvian scallop (Li et al. 2018) using Blastn (BLAST 2.2.28 +). Filtered markers were then used to construct the linkage map using JoinMap (version 4.0) (Stam 1993). Markers were categorized into three groups according to their phase types (lm × ll and nn × np as species-specific markers, and hk × hk as species-shared markers). Only species-specific markers were used for map construction. Markers were assigned to different linkage groups with logarithm of odds (LOD) cut-off of 10. The recombination frequencies were converted into map distances (centimorgans) through the Kosambi mapping function (Kosambi 1944). The expected genetic map lengths (G_e) were estimated according to Chakravarti et al. (1991) with the following formulas: G_e1 = G_oa + 2S (G_oa: the observed genetic map length; S: average marker interval); \(G_{{{\text{e}}1}} = \sum\nolimits_{i}^{n} {L_{i} \times (m_{i} + 1)/(m_{i} - 1)}\) (L_i: the observed length of LGi; m_i: the marker number of LGi; n: the number of LGs). Map coverage was measured as the ratio of observed genetic map length and the average of expected lengths. Finally, the genetic linkage map was graphically presented through MapChart (version 2.2) (Voorrips 2001).

The distributions of the five growth-related traits, including body weight (BW), adductor muscle weight (AMW), shell length (SL), shell height (SH), and shell width (SW), were checked whether they followed the normal distributions by the Kolmogorov–Smirnov test using MATLAB. The phenotypic correlations among the five growth traits were estimated by the Pearson's correlation coefficients. The progenies were characterized with two distinct types of shell coloration patterns (Fig. 1). The shell color was analyzed as a binary trait, and its segregation ratio was tested by χ² test.

QTL mapping was conducted for both growth-related traits and shell color with R/qtl (Broman et al. 2003) using the following procedures. The best multiple QTL model was determined using Haley–Knott regression (Haley and Knott 1992). Subsequently, automated stepwise model selection procedures were used to scan additive and epistatic QTL at each step (Manichaikul et al. 2009). Significant QTLs were identified with the LOD thresholds (experiment-wise α = 0.05) determined by 1000 permutations. Bayesian 95% credible intervals were calculated for each QTL as a measure of location uncertainties (Broman and Sen 2009). Also, single QTL and pairwise marker interactions were checked to confirm the selected model. The markers were mapped to the genome assemblies of the bay scallop (BioProject No. PRJNA428789) and the Peruvian scallop (Li et al. 2018) using Blastn (BLAST 2.2.28 +), genes located within the QTLs were identified as candidates. The full linkage maps of the two species, the sequences of markers, the phenotypic data, and the codes for QTL mapping have been deposited in https://www.figshare.com with the Digital Object Identifier (DOI) https://doi.org/10.6084/m9.figshare.8246867.

We acknowledge the Grant support from National Natural Science Foundation of China (U1706203, 31172404 and 31572618), Taishan Scholar Project Fund of Shandong Province of China, and Youth Talent Program Supported by Laboratory for Marine Fisheries Science and Food Production Processes, Pilot National Laboratory for Marine Science and Technology (Qingdao) (2018-MFS-T07).

SW and CW: conceived and designed the study; JM, ZY, HP, and LY: prepared the materials and constructed the sequencing libraries; JM and QZ: conducted genetic mapping and other data analysis; ZB: advised and coordinated the study. All the authors contributed to manuscript writing, reviewing, and approved the fnal version for submission.

The authors declare that they have no confict of interest.

We declare that all applicable international, national, and/or institutional guidelines for sampling, care, and experimental use of organisms for the study have been followed and all necessary approvals have been obtained.