Construction of a Genetic Linkage Map Based on SNP Markers, QTL Mapping and Detection of Candidate Genes of Growth-Related Traits in Pacific Abalone Using Genotyping-by-Sequencing

Pacific abalone (Haliotis discus hannai) is a commercially important high valued molluscan species. Its wild population has decreased in recent years. It is widely cultured in Korea. Traditional breeding programs have been implemented for hatchery production of abalone seeds. To obtain more genetic information for its molecular breeding program, a high-density linkage map and quantitative trait locus (QTL) for three growth-related traits was constructed for Pacific abalone. F1 cross population with two parents were sampled to construct the linkage map using genotyping by sequencing (GBS). A total of 664,630,534 clean reads and 56,686 SNPs were generated and 3,345 segregating SNPs were used to construct a consensus linkage map. The map spanned 1,747.023 cM with 18 linkage groups and an average interval of 0.55 cM. QTL analysis revealed two significant QTL in LG10 on the consensus linkage map of each growth-related trait. Both QTLs were located in the telomere region of the chromosome. Moreover, four potential candidate genes for growth-related traits were identified in the QTL region. Expression analysis revealed that these identified genes are involved in growth regulation of abalone. The newly constructed genetic linkage map, growth-related QTLs and potential candidate genes identified in the present study can be used as valuable genetic resources for marker-assisted selection (MAS) of Pacific abalone in molecular breeding program.


INTRODUCTION
Abalones are marine molluscan species distributed worldwide along coastal waters in tropical and temperate regions (Gordon and Cook, 2013). Pacific abalone (Haliotis discus hannai) is the most commercially important invertebrate species in South Korea as well as in South-East Asian countries. Because of its extraordinary food value due to the presence of bioactive molecules and high price, commercial aquaculture of abalone has expanded widely (Park and Kim, 2013;Hossen et al., 2021). It is reared in sea cages and extensive case farming has distributed along small bays and islands of South, South-West and West coasts of Korea (Choi et al., 2015). Commercial-scale production of abalone in Korea began in the 2000s. As of 2015, abalone production reached about 10,090 metric tons. Its production was doubled in the next 4 years, reaching about 20,100 metric tons in 2019 with commercial value of about 536 million USD (KOSTAT, 2020). However, aquaculture production of abalone is facing difficulties due to its low growth rate, summer mortality and frequent outbreaks of infectious diseases. A total of 684 hatcheries are producing abalone seeds in Korea using conventional selective breeding methods (Choi et al., 2015). Recent studies have shown that successful aquaculture mostly depends on molecular breeding for rapid growth, bigger size, disease resistance and higher survival of target species . Hence, it is essential to ascertain molecular breeding for Pacific abalone.
Growth is the most preferred economic important trait for aquaculture species including abalone. Substantial improvement of abalone seed production has been achieved using conventional phenotype-based selection for growth. However, such method is labor-intensive, time consuming, error prone and influenced by environmental factors (Gjedrem and Baranski, 2009). Molecular breeding using marker-assisted selection (MAS) can increase the accuracy of selection and shorten the operation time (Xu et al., 2017). Economically related traits such as growth, meat quality and immunity are mainly controlled by quantitative trait loci (QTLs). QTLs are defined as chromosomal regions linking single gene or genes clusters (Yue, 2014). With progress and growing application of contemporary biotechnology, markerassisted selection (MAS) is becoming popular over phenotypebased selection in animal breeding program. Markers linked to QTLs are used in MAS to accelerate the molecular breeding technique that can increase the accuracy of selection and achieve genetic improvement through direct and prompt selection (Feng et al., 2018;Song et al., 2018).
Mapping of QTL is a principal step in molecular breeding program. Genetic linkage map is an essential tool for genomic research as it can provide a basis for QTL mapping (Tian et al., 2015). The availability of a bulk number of molecular markers is a prerequisite to construct a high-density genetic linkage map. Several molecular markers such as express sequence tag (EST), single sequence repeats (SSR), restriction fragment length polymorphism (RFLP), amplified fragment length polymorphism (AFLP), randomly amplified polymorphic DNA (RAPD) and single nucleotide polymorphism (SNP) can be used to construct a genetic linkage map. Among these, SNPs are the most abundant markers located in the genome, and appropriate for broadscale genotyping and construction of a high-density linkage map (Lindblad-Toh et al., 2000). A number of linkage maps have been reported for Pacific abalone using DNA-based low-density SSR markers (Sekino and Hara, 2007), AFLP markers (Qi et al., 2007) and a combination of AFLP, RAPD, and SSR markers (Liu et al., 2006(Liu et al., , 2007. Even so, these linkage maps are limited by the lack of dense polymorphic molecular markers. Recently introduced next-generation sequencing (NGS) technologies are cost-effective, easy, and rapid to detect genome-wide molecular markers (such as SNPs) for construction of a high-density linkage map. Among various NGS technologies, genotype-by-sequencing (GBS) is a new method to perform sequence-based genotyping by preparing libraries of thousands of DNA fragments. Thus, GBS technology provides a large quantity of quality data that can be used to construct a fine map using SNP markers from a reference genome to detect QTLs (Salazar et al., 2019). After the introduction of GBS, QTL mapping using GBS has been performed for several aquaculture species such as snapper (Ashton et al., 2019), Pacific oyster , Pacific bluefin tuna (Uchino et al., 2018) and scallop (Jiao et al., 2014).
In the present study, GBS approach was used for largescale identification of SNP markers and construction of a SNPbased high-density linkage map for Pacific abalone. Further QTL mapping was performed to detected SNPs associated with growth-related traits. Finally, candidate genes were identified in the QTL region that might regulate growth-related traits of Pacific abalone.

Mapping Family
An F1 family of Pacific abalone was produced in May 2019 for genetic mapping. Female and male parents were obtained from two geographically distant zones in Korea: a smallersized female ( in West coast of Korea. The parents were maintained in a commercial abalone hatchery, Tou-Jeong Soosan abalone hatchery, Dolsan-eup, Yeosu, Korea. Abalones were reared in recirculating sea water tanks. After reaching at full maturity, induced spawning was performed by heat-induction and UV-irradiation (Uki and Kikuchi, 1984;Leighton, 2008). Briefly, fully mature abalones were exposed to sunlight for one hour in an upside-down orientation and another half an hour in an upside-up orientation. Heat induced abalones were then kept in a spawning tank with UV-irradiated seawater and aeration until spawning. After fertilization and production of abalone larvae, 15,000 progenies were reared under commercial condition at the larval rearing unit with recirculating sea water in Tou-Jeong Soosan abalone hatchery. These larvae were ensured to have optimal grow-out densities.

Husbandry of Abalone Larvae and Monitoring of Water Quality
A 1.5 L min −1 sea water flow rate was maintained continuously in the rearing tanks. Water quality parameters such as temperature, dissolved oxygen, pH and salinity were measured daily with a YSI digital water quality meter (Pro 10102030; Xylem Inc., United States) (Supplementary Table 1). Abalone larvae were fed every day in spring-summer season and every other day in winter season to satiation with commercial abalone feed (Gold Ocean Premium Abalone Feed, JOEUN Corp., Wando, Korea). The rearing tank was cleaned every day in spring-summer season and at every other day interval in winter season by changing water and washing with sea water flash.

Ethics Statement
Animal experiments were conducted in accordance with guidelines of the Institutional Animal Care and Use Committee of Chonnam National University (CNU IACUC) with a permission number CNU IACUC-YS-2020-5. Animals were used according to Article 14th of the Korean Animal Protection Law of the Korean government. Animals were cared for in accordance with Guidelines for Animal Experiments of Chonnam National University.

Sample Collection and Trait Measurement
Twelve months after hatching in May 2020, 96 F1 offspring (48 largest and 48 smallest) were randomly selected and sampled as described previously (Baranski et al., 2008;Li and He, 2014). For linkage analysis, three quantitative growth traits of abalone, i.e., total weight (TW), shell length (SL) and shell width (SW), were measured using an electric balance (SPX223KR, OHAUS, United States) and a digimatic caliper (CD-P20S, Mitutoyo Corp., Japan) to the nearest mg and mm, respectively (Supplementary Table 2). A portion of the nondestructive cephalic tentacle was seized from each juvenile abalone for extraction of genomic DNA (gDNA), washed with phosphate buffered saline (PBS), immediately flash frozen in liquid nitrogen and stored at −80 • C until gDNA extraction. After taking cephalic tentacles, juvenile abalones were tagged and returned to a new grow-out tank. Cephalic tentacles of parent abalone were also collected previously during the breeding experiment and stored at -80 • C after flash frozen in liquid nitrogen.

Genotyping-by-Sequencing (GBS) Library Preparation, Sequencing and SNP Genotyping
Extraction of Genomic DNA Genomic DNA of each F1 progeny and parent abalone were extracted from cephalic tentacle using the CTAB (Cetyl trimethylammonium bromide) method as describe previously (Winnepenninckx, 1993;Arseneau et al., 2017) with slight modifications. Briefly, samples were cryogenically ground to fine power in a bead beater. After adding 200 µL of CTAB extraction buffer (0.1 m Tris−HCl pH 8.0, 20 mm EDTA pH 8.0, 1.4 M NaCl, 2% CTAB w/v) containing 73 µL 2-mercaptoethanol to each 100 mg homogenized tissue, the mixture was vortexed thoroughly. The homogenate was then incubated at 65 • C in a water bath for 1 h and then centrifuged at 14,000 g for 5 min. The supernatant was then transferred to a new tube. After adding 5 µL of RNase A, the mixture was incubated at 37 • C for 20 min. An equal volume of phenol/chloroform/isoamyl alcohol (25:24:1) was added to the mixture. After vortexing for 5 s, the mixture was incubated at room temperature for 5 min and then centrifuged at 13,000 g for 5 min to separate aqueous phases. The clear aqueous upper phase was then transferred to a new tube. After adding 200 µL of chloroform/isoamyl alcohol (24:1), previous steps were repeated and the upper clear aqueous phase was transferred to a new tube. Later, DNA samples were precipitated by adding 2.5 volumes of ice-cold 100% ethanol and incubated at −70 • C for 30 min. After the incubation, samples were centrifuged at 13,000 g for 30 min at −4 • C. The supernatant was then decanted without disturbing the pellet. The pellet was subsequently washed with 500 µL ice-cold 70% ethanol. Ethanol was decanted and residual ethanol was removed by air-drying. Finally, DNA was dissolved in 20 µL TE buffer (10 mM Tris, pH 8, 1 mM EDTA). The integrity of each DNA sample was evaluated by 1.2% agarose gel electrophoresis and concentration of the DNA was determined with a Nano Photometer spectrophotometer (IMPLEN, United States). The ratio of absorbance at 260/280 nm was used to assess the purity of DNA.

Construction of GBS Library and Sequencing
Construction of GBS library and Illumina HiSeq X sequencing were performed at SEEDERS, Korea using genomic DNA of each female and male parents and 90 progenies. A GBS library was constructed following the protocol described previously (Mascher et al., 2013) by double digestion with two restriction enzymes: PstI and MspI. Chronological steps of GBS library construction included adaptor annealing, DNA digestion with PstI and MspI, adaptor ligation, sample pooling, DNA purification and multiplexed PCR. The quality of the constructed library was analyzed through a quality control check and agarose gel electrophoresis. Pooled GBS library was then sequenced using a HiSeq X platform (Illumina, Inc., United States) by the pairedend read method.

Sequence Data Analysis
Raw sequences were demultiplexed into 96 samples according to barcode sequences. Barcode and adapter sequences were removed using Cutadapt (version 1.8.3) software (Martin, 2011). Low quality sequences were trimmed using DynamicTrim and LengthSort programs of SolexaQA (v.1.13) package (Cox et al., 2010). For DynamicTrim, a phred score ≥20 was used as the criterion. For LengthSort, a read length of ≥25 pb was applied. DynamicTrim removed low-quality bases at either ends of short reads based on the phred score to carry out a purification process for high-quality clean reads. LengthSort removed the reads, from which excess bases were then cut with DynamicTrim.

Raw SNP Detection and Generation of SNP Matrix
Cleaned sequences were aligned to the reference genome sequence of Pacific abalone consisting of 80,032 scaffolds (Supplementary Table 3) using a Burrows-Wheeler Aligner (BWA) program version 0.6.1-r104 (Li and Durbin, 2009). A BAM format file was generated to detect raw SNP (In/Del) using default values with the following options: a seed length (l) of 30, a maximum difference in the seed (-k) of 1, number of threads (-t) of 16, a mismatch penalty (-M) of six, a gap open penalty (-O) of 15, and a gap extension penalty (-E) of eight. Raw SNPs (In/Del) were detected and consensus sequences were extracted using SAMTools version 0.1.16 ). Prior to raw SNP detection, SNP validation was carried out using an in-house script of SEEDERS (Kim et al., 2014). To compare SNPs among analytic targets, an integrated SNP matrix was produced for samples. Based on respective coordinates, SNPs were classified as homozygous (SNP read rate ≥ 90%), heterozygous (40% ≤ SNP read rate ≤ 60%) and ambiguous (20% ≤ SNP read rate ≤ 40% or 60% < SNP read rate < 90%).

SNP Filtering, Genotyping, and Construction of Genetic Linkage Map
Genetic linkage maps were constructed using the JoinMap 4.1 program (Van Ooijen, 2006) with population type crosspollination (CP) (Grattapaglia and Sederoff, 1994) under the condition of logarithm of odds (LOD) of 6.0 or higher with a maximum distance of 30 cM. To construct an integrated map of male and female parents, all segregating markers that showed polymorphism in at least one parent were used in the JoinMap configuration for CP mode (lmxll, nnxnp, hkxhk). Ratios of marker segregation were calculated using Chi-square test. Markers that satisfied the expected Mendelian segregation ratio were included for mapping. Filtered markers contained different segregation types, including < nn x np >, < lm x ll > and < hk x hk >. Markers showing significantly distorted segregation (p < 0.05) were excluded from map construction. Markers were grouped with a minimum logarithm of odds (LOD) score of 6.0 and a recombination frequency of 0.45. A regression mapping algorithm was used to build the linkage map. Map distances were calculated in centiMorgans (cM) according to the Kosambi mapping function (Kosambi, 2016). For all other calculations, default settings were used.

QTL Mapping
QTL analysis was performed using R/qtl (Broman et al., 2003). Initially, an interval mapping (IM) was performed to identify potential QTLs. Later, a multiple QTL model (MQM) using the stepwiseqtl function was performed for forward and backward selections with a check for interactions to identify multiple QTLs. Based on IM results, the maximum number of QTLs to search for each trait was set to be three. Significance thresholds were determined using 1,000 permutations (p < 0.05). Results from QTL analysis were used to construct a QTL map, and their positions were used in a default model.

Identification of Potential Candidate Genes in the QTL Region
To identify potential candidate genes associated with growth was investigated in the QTL region. Genes encoded in the four scaffolds regions located in LG10 were BLAST searched against the reference genome sequence of Pacific abalone (Nam et al., 2017) and sequences of the NCBI non-redundant Molluscan database.

Quantitative Real Time PCR (qRT-PCR) Analysis of Candidate Genes
The mRNA expression of candidate genes in fast growth (49.38 ± 3.36 g) and slow growth (14.02 ± 3.27 g) abalone groups were investigated by qRT-PCR. The experimental abalones were two years old and obtained from an independent family having ten individuals in each group. Total weight (TW) was measured to the nearest mg and, shell length (SL) and shell width (SW) were measured to the nearest mm (Supplementary Table 4). From each abalone, tissue samples such as cephalic tentacle (CT), epipodial tentacle (ET), muscle (MUS), mantle (MNT), heart (HRT), gill (GIL), hepatopancreas/digestive gland (DG), ovary (Ov), and testis (Te) were collected, snap frozen in liquid nitrogen, and stored at -80 • C until total RNA extraction. A separate qRT-PCR analysis of candidate genes was also conducted for parent abalones and presented in the Supplementary Figure 1.

Total RNA Extraction and cDNA Synthesis
Total RNAs were extracted from all sampled tissues using an ISOSPIN Cell & Tissue RNA kit (Nippon Gene, Tokyo, Japan). First-strand cDNAs were synthesized from total RNAs using and oligo (dT) primer (OdT) (Sigma) and a Superscript R III First-Strand cDNA synthesis kit (Invitrogen, United States). All steps of RNA extraction and cDNA synthesis were conducted following the manufacturer's instructions.

qRT-PCR Analysis
Quantitative real-time PCR (qRT-PCR) assay was carried out using 2x qPCRBIO SyGreen Mix Lo-Rox kit (PCR Biosystems,     United Kingdom) on a LightCycler R 96 System (Roche, Germany) as described previously . Briefly, qRT-PCR was executed using 1 µL of cDNA template in a total reaction volume of 20 µL containing 1 µL (10 pmol) of each forward and reverse primer (Supplementary Table 5), 10 µL SyGreen Mix, and 7 µL ultra-pure water. Triplicate reactions were performed for target and reference genes for each tissue sample. PCR amplification conditions were: pre-incubation at 95 • C for 3 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 30 s. At the end of each cycle, a fluorescence reading was recorded for quantification. Specific amplification of each cDNA sample was verified by melting curve analysis and gel electrophoresis of the qRT-PCR product. Relative gene expression was quantified based on cycle threshold using the 2 − CT method (Livak and Schmittgen, 2001). Housekeeping β-actin gene (GenBank accession no. AY380809) was used as an internal reference to normalize the gene expression. Relative mRNA expression levels were analyzed statistically and expressed as mean ± SEM (standard error of the mean). Changes in relative mRNA expression levels between slow and fast growth groups of abalone were analyzed by Student's t-test. Statistical significance was set at p < 0.05. All data were analyzed, and graphs were generated using GraphPad Prism 5.1 software 1 .

Values of Phenotypic Traits
A total of 90 Pacific abalone F1 progenies were used to generate the linkage map. Average values of total weight, shell length and shell width were 15.58 ± 6.80 g, 49.30 ± 9.50 mm and 33.36 ± 6.21 mm, respectively. These growth traits were strongly correlated with each other (r = 0.968-0.990; p < 0.01 for all traits). Shell length and shell width showed the highest correlation (r = 0.990) (Supplementary Table 6).

GBS Sequencing
A total of 752,414,780 reads were obtained from 90 progenies (among 96 sample, 6 progenies were discarded because of poor DNA quality) and parents with an average length of 151 bp. The total length was 113,614,631,780 bp with a GC content of 45.25% and an average Q30 ratio of 91.69% (Table 1). A total of 664,630,534 (91.69%) reads passed the quality filter with a total length of 80,519,620,461 bp ( Table 2). There were 12,651,462 reads from the female parent and 43,000946 reads from the male parent. Offspring per capita had an average of 7,498,108 reads (Table 3).

SNP Discovery and Genotyping
A total number of 815,970 putative SNPs were detected. Among these SNPs, 144,819 and 200,695 putative SNPs were found in the female parent and the male parent, respectively ( Table 3). There were 73,025 and 104,357 homozygous SNPs, and 32,525 and 49,280 heterozygous SNPs in the female and male parents, respectively. In the offspring, there were 58,985 homozygous SNPs and 24,120 heterozygous SNPs ( Table 4). Numbers of homozygous and heterozygous SNPs in all samples are presented in Figure 1, with parents having more SNPs than the offspring.  All sequences containing SNPs are presented in Supplementary  Table 7. After initial genotyping, the CP model in JoinMap 4.0 was used to select segregation types nn x np (26,487 SNP markers), lm x ll (14,602 SNP markers) and hk x hk (15,597 SNP markers), resulting in a total of 56,686 SNPs that were polymorphic in at least one parent and one progeny ( Table 5).

Characteristics of Genetic Mapping
A total of 3,345 markers (1,251 female-specific, 944 malespecific and 1,150 anchor) and phenotypic data of F1 population developed from two geographically distant parents were used to construct a linkage map. A total of 18 linkage groups (LGs) were generated from 3,314 markers (Figure 2). The consensus linkage map had a map distance of 1747.023 cM. The average distance between adjacent markers was 0.55 cM, with the lowest distance of 0.32 cM and the highest distance of 0.85 cM ( Table 6). The number of markers per LG varied from 98 to 297, with a mean of 184.11 markers per LG.

QTL Analysis
The integrated linkage map with 3,314 effective SNPs in 18 LGs of Pacific abalone was used to detect QTLs for three growthrelated traits. Results are presented in Table 7. Two suggestive overlapped QTLs for three growth-related traits were detected in one linkage group (LG 10, Figure 2) at the position of 0 and 8.104 cM in the telomere region of the chromosome. Two SNPs were located on each growth-related trait, with LOD values ranging from 13.490 to 16.513 for shell length, from 30.886 to 33.883 for shell width, and from 25.962 to 29.006 for total weight (Table 7 and Figure 3).

Identification of Potential Candidate Gene in the QTL Region and Expression Analysis
Four candidate genes were identified encoded on four scaffolds (HDSC01365CG00050, HDSC01365CG00060, HDSC01365CG00070 and HDSC01726CG00010) in the QTL region of LG10. Associated GO categories of these genes included activating signal cointegrator 1 complex subunit (ascc1), ecdysoneless homolog (Ecd), galectin-4 (Gal-4) and calmodulin-B (CaM-B). Details of these identified candidate genes are presented in Table 8.
The mRNA Expression analysis of identified candidate genes in slow and fast growth groups of abalones showed that mRNA expression levels of these genes were significantly upregulated or downregulated in experimental groups of abalones (Figure 4). Expression levels of ascc1 were found to be significantly higher in CT, ET, and MNT of abalones with fast growth than in abalone with slow growth (Figure 4A), whereas expression of Ecd was downregulated in abalone with fast growth, i.e., its expression levels were significantly lower in CT, ET, and MNT of abalone with fast growth than in abalones with slow growth (Figure 4B). mRNA expression levels of Gal-4 gene were significantly higher in ET, MNT, and HRT of abalones with fast growth than in abalones with slow growth (Figure 4C). Expression levels of CaM-B in ET and MNT of abalones with fast growth were also significantly higher than in abalones with slow growth (Figure 4D).

DISCUSSION
The present study reported a SNP marker based genetic linkage map that was constructed for Pacific abalone using GBS approach. GBS is an efficient and cost-effective method that could be used to develop SNP markers, construct genetic linkage map and select marker-based complex traits of many species with or without a reference genome (Elshire et al., 2011;Nie et al., 2017;Liu et al., 2020). This is the first report of QTL mapping using GBS to screen high quality and abundant SNP markers for constructing a genetic linkage map for Pacific abalone. Previously, several genetic linkage maps have been constructed for Pacific abalone using SSR markers (Sekino and Hara, 2007), AFLP markers (Qi et al., 2007), and a combination of AFLP, RAPD, and SSR markers (Liu et al., 2006). A growth-related QTL has also been reported for Pacific abalone using a combination of RAPD, AFLP and SSR markers (Liu et al., 2007). SNP markers are most abundant in an organism. They more appropriate than AFLP, RAPD, or SSR markers for genetic studies, especially for constructing a genetic linkage map, identification of QTLs and population genetic analysis (Vignal et al., 2002;Yu et al., 2011;Ren et al., 2016;Yáñez et al., 2016). To date, the development of SNP markers using GBS or RAD sequencing of an organism is simple and quick to construct a genetic linkage map and perform QTL analysis for molluscan aquaculture species with large genomes, such as pearl oyster (Li and He, 2014;Shi et al., 2014) small abalone (Ren et al., 2016), Pacific oyster Shi et al., 2020), scallop (Jiao et al., 2014) and Manila clam (Nie et al., 2017). In this study, a genetic linkage map was constructed using 3,345 SNP markers that were filtered    (Arai et al., 1982). The total length of the consensus map was 1747.023 cM with an average inter-location interval of 0.55 cM. This value of inter-location interval was higher than that of Pacific oyster (0.31 cM) (Shi et al., 2020), scallop (0.41 cM) (Jiao et al., 2014), or pearl oyster, Pinctada fucata martensii (0.39 cM) (Shi et al., 2014), but lower than that of small abalone (0.59 cM) (Ren et al., 2016) or pearl oyster, Pinctada fucata (1.14 cM) (Li and He, 2014). In this study, the target quantitative trait was growth. It was measured as total weight, shell length and shell width Asterisks indicated significant differences between slow and fast growth groups (p < 0.05). CT, cephalic tentacle; ET, epipodial tentacle; MUS, muscle; MNT, mantle; HRT, heart; GIL, gill; DG, digestive gland; Ov, ovary; Te, testis.
( Table 2). As growth is directly related to the production performance of aquaculture species, it is one of primary targets for marker assisted selection in genetic breeding programs. Two outbred geographically distant parents were used to produce F1 population in this experiment. Due to a longer reproductive cycle, it is difficult to produce F2 progenies for molluscan species. Therefore, most genetical mapping and QTL studies have used F1 cross population for molluscan species (Baranski et al., 2006;Liu et al., 2006;Sekino and Hara, 2007;Shi et al., 2010Shi et al., , 2014Zhan et al., 2012;Vervalle et al., 2013;Jiao et al., 2014;Ren et al., 2016;Li et al., 2018). Two suggestive growth-related QTLs in each quantitative growth trait were detected in the linkage group LG10 of Pacific abalone. Fewer numbers of growth-related QTLs have also been reported for aquaculture species, such as one QTL for Pacific bluefin tuna (Uchino et al., 2018), two QTLs for scallop (Jiao et al., 2014) and three QTLs for kelp grouper (Kessuwan et al., 2016). The two QTLs detected in the present study overlapped for three growth-related traits, suggesting that these growth traits might be controlled by a common region. Overlapping QTL region has also been reported for small abalone (Ren et al., 2016) and Asian seabass (Xia et al., 2014). In the present study, both detected QTLs were located in the telomere of chromosome 10. Since recombination rates are considered to be higher in telomeric region than in centrometric regions of chromosomes (Sakamoto et al., 2000;Anderson et al., 2012), a higher residual heterozygosity level is expected for telomeric regions (Fraslin et al., 2018). Therefore, the present result of the location of QTL in telomere region of chromosome suggests more advanced segregating population and it should be considered as a more robust result.
Beside for detecting genetic markers associated with targeted traits, QTL mapping is also useful for identifying candidate genes that regulate targeted traits (Mackay et al., 2009;Liu et al., 2020). Using BLAST search against the reference genome and sequences in the NCBI non-redundant Molluscan database, four potential candidate genes were identified in scaffolds of the QTL region (Table 8). These genes were mainly involved in signal transduction, regulation of transcription factor and growth. Although whether these genes are directly or indirectly responsible for growth of abalone remains unknown, there are some evidences showing that these genes are involved in the growth of several animal species. For example, Ecd is an evolutionarily conserved protein that plays a crucial role in embryonic development and cell growth (Garen et al., 1997;Kainou et al., 2006). Galectin gene plays an important role in diverse fundamental cellular processes such as growth regulation, lipid raft stabilization, protein apical trafficking, adhesion/migration and mediator release, cell adhesion, and wound healing (López-Lucendo et al., 2004;Cao and Guo, 2016). It has also been reported that galectins contribute to larval development of Pacific abalone (Sandamalika and Lee, 2020), muscle development in Drosophila (Pace et al., 2002) and chick muscle (Cooper and Barondes, 1990). The Gal-3 gene in Siniperca chuatsi was significantly upregulated in groups of fish with higher growth performance, suggesting that protein glycosylation might contribute to the growth performance of S. chuatsi (He et al., 2020). In general, CaM genes and their associated kinases are involved in increased growth and anabolic functions. Thus, their expression levels are expected to be higher in faster growing organisms (Kocmarek et al., 2014). It is evident that CaM genes are involved in regulatory mechanisms related to growth, calcium metabolism, biomineralization and shell formation of Mollusca (Li et al., 2004;Huang et al., 2007;Lim et al., 2018). Genes encoding insulin-like growth factor binding proteins, molluscan insulin related peptide, epidermal growth factor and myostatin are considered to be directly associated with growths of abalone and molluscan species. However, these genes were not detected in the QTL region of Pacific abalone. Results of the present study indicate that mapping of new SNP markers and assemblage of gene information in the QTL region can facilitate the identification of specific regions in the genome that contain growth-related genes. Furthermore, the differential expression patterns of ascc1, Ecd, Gal-4, and Cal-B between slow and fast growth abalones suggest that these genes might be associated with growth regulation of abalone.
In conclusion, QTL mapping using GBS in Pacific abalone revealed two new SNP markers on one linkage group that were significantly associated with growth traits. In addition, four candidate genes were identified on the same LG of these QTLs; which were found to be upregulated or downregulated in abalones with fast growth or slow growth. These genomic resources can be used in selective breeding program of Pacific abalone. These data can also be used as basic data for further genomic investigation of Pacific abalone.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10. 5061/dryad.msbcc2fz3.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Chonnam National University (CNU IACUC).

ACKNOWLEDGMENTS
The authors gratefully acknowledge the support of Dongik Choi (President, Tou-Jeong Soosan abalone hatchery, Dolsan-eup, Yeosu-si, Jeollanam-do, South Korea) for conducting the field experiment in his hatchery facilities.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.

2021.713783/full#supplementary-material
Supplementary Figure 1 | Relative mRNA expression levels of (A) ascc1, (B) Ecd, (C) Gal-4 and (D) CaM-B in different tissues of female and male parents of abalone used for generating of F1 progeny of QTL. Asterisks indicated significant differences between slow and fast growth group (p < 0.05). CT, cephalic tentacle; ET, epipodial tentacle; MUS, muscle; MNT, mantle; HRT, heart; GIL, gill; DG, digestive gland; GND, gonad.  Table 4 | Trait data of slow growth and fast growth abalones used in qRT-PCR expression analysis.

Supplementary
Supplementary Table 5 | List of primers used for quantitative real-time polymerase chain reaction (qRT-PCR) for expression analysis of QTL candidate genes in Pacific abalone with slow and fast growth.