A Genome-Wide Association Study Identifies Candidate Genes Associated With Shell Color in Bay Scallop Argopecten irradians irradians

Molluscan shell color has consistently drawn attention for its abundant diversity and commercial use in shellfish breeding projects. Recently, two new strains of bay scallop (Argopecten irradians irradians) with different shell colors as marked phenotypic traits have been artificially bred to improve their economic values; however, the inheritance mechanism of their shell pigmentation is still unclear. In this study, a genome-wide association study (GWAS) was conducted to determine the genetic basis of shell color in bay scallops utilizing 29,036 high-quality single-nucleotide polymorphisms (SNPs) derived from 80 purple-red (PP) and 80 black-brown (BP) shell color individuals. The result of the GWAS showed that 469 SNPs (p <1.72E−6) significantly associated with shell color were mainly distributed in chromosome 7. The top three SNPs (i.e., chr7-12764003, chr7-13213864, and chr7-11899306) are located in the genic region of G-protein-coupled receptor-like 101 (GRL101), polyketide synthase 1 (PKS1), and phosphoinositide phospholipase C (PLC1), which have been widely reported to be involved in pigmentation. Successfully, the top three SNPs were verified in another non-breeding bay scallop population. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses obtained 38 GO terms covering 297 genes and aggregating pathways involving 252 annotated genes. Specifically, the expression profiles of the top three identified candidate genes were detected in mantles of PP and BP individuals by real-time quantitative reverse transcription PCR. The significantly higher expression levels of GRL101 (6.43-fold) and PLC1 (6.48-fold) in PP, and PKS1 (12.02-fold) in BP implied that GRL101 and PLC1 potentially functioned in PP shell coloration, and black pigmentation in BP might be principally regulated by PKS1. Our data provide valuable information for deciphering the phenotype differences of shell color in the bay scallop.

Molluscan shell color has consistently drawn attention for its abundant diversity and commercial use in shellfish breeding projects. Recently, two new strains of bay scallop (Argopecten irradians irradians) with different shell colors as marked phenotypic traits have been artificially bred to improve their economic values; however, the inheritance mechanism of their shell pigmentation is still unclear. In this study, a genome-wide association study (GWAS) was conducted to determine the genetic basis of shell color in bay scallops utilizing 29,036 high-quality single-nucleotide polymorphisms (SNPs) derived from 80 purple-red (PP) and 80 black-brown (BP) shell color individuals. The result of the GWAS showed that 469 SNPs (p < 1.72E−6) significantly associated with shell color were mainly distributed in chromosome 7. The top three SNPs (i.e., chr7-12764003, chr7-13213864, and chr7-11899306) are located in the genic region of G-protein-coupled receptor-like 101 (GRL101), polyketide synthase 1 (PKS1), and phosphoinositide phospholipase C (PLC1), which have been widely reported to be involved in pigmentation. Successfully, the top three SNPs were verified in another non-breeding bay scallop population. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses obtained 38 GO terms covering 297 genes and aggregating pathways involving 252 annotated genes. Specifically, the expression profiles of the top three identified candidate genes were detected in mantles of PP and BP individuals by real-time quantitative reverse transcription PCR. The significantly higher expression levels of GRL101 (6.43-fold) and PLC1 (6.48-fold) in PP, and PKS1 (12.02-fold) in BP implied that GRL101 and PLC1 potentially functioned in PP shell coloration, and black pigmentation in BP might be principally regulated by PKS1. Our data provide valuable information for deciphering the phenotype differences of shell color in the bay scallop.

INTRODUCTION
The shell of mollusks with diverse forms and alternative colors has always appealed to consumers and scientists (Comfort, 1950). Shell color of marine mollusks, an obvious and marked genetic-based phenotypic trait, has been proven to be closely related to production traits (Cong et al., 2014;Zhang et al., 2016). Therefore, shell color has been widely applied in the selective breeding project of economic shellfish, such as oysters (Xu et al., 2019;Han et al., 2020), mussels (Innes and Leslie, 1977;Li et al., 2014), scallops (Petersen et al., 2012;Ding et al., 2015), clams Nie et al., 2020), and abalones Hoang et al., 2017). To reveal the mechanism of shell color formation, several investigations have been conducted and found multiple factors, ranging from environmental factors to inner genetic factors, contributing to the higher variability of shell color in mollusks (Underwood and Creese, 1976;Kraeuter et al., 1984). For example, researches demonstrated that water temperature differences (Heller, 1992), salinity variations (Sokolova and Berger, 2000), and diet sources (Marchais et al., 2017) could affect shell color formation in mollusks. On the aspect of genetic basis, shell color has been demonstrated to be determined by one or a small number of major gene(s) in some shellfish species (Cole, 1975;Palmer, 1985). For instance, experimental crosses of Manila clam, Ruditapes philippinarum, verified that pigmented coloration in two valves is controlled by at least two genes (Peignon et al., 1995). In noble scallop Chlamys nobilis, a one-locusthree-allele model was proposed to elucidate the distribution of four different color variants, and the brown shell color is controlled by a recessive allele distinct from the other colors . Investigation in Pacific oyster, Crassostrea gigas, demonstrated that shell pigmentation is controlled by two genetic loci, with one responsible for the secretion of pigments and the other responsible for the distribution mode of pigmentation (Xu et al., 2019). Furthermore, transcriptomes and digital gene expression analysis of four different color clams, Meretrix meretrix, suggested that several potential genes and the Notch pathway played a crucial role in its shell color patterning (Yue et al., 2015). The same method was also exerted in two extreme color phenotypes of Yesso scallop, Patinopecten yessoensis, and 25 significantly differential expression genes were identified for unraveling shell color differences (Sun et al., 2016). These researches provided pivotal loci/genes responsible for the formation of shell coloration and served as candidate markers for mollusks breeding projects.
Biological pigments are considered as the key element of shell color diversity, and the distinct color is perceived because of different types of pigments, such as melanins, porphyrins, tetrapyrroles, and carotenoids, which had been identified in shellfish shells (Comfort, 1951;Stemmer and Nehrke, 2014). Melanins were first derived by oxidation and polymerization of tyrosine in animals or phenolic compounds in lower organisms (D'Iischia et al., 2013). Using the ultraviolet and infrared radiation spectral analysis, melanin was confirmed as the black pigment extracted from the scars where adductor muscle jointed shells in oysters . Eumelanin and pheomelanin, two forms of melanin, were verified to comprise the dark brown shell pigments with ∼76.6 and 23.4% components via spectrophotometry scanning and high-performance liquid chromatography in Yesso scallop (Sun et al., 2017). As for porphyrins, the cyclic structure of tetrapyrroles is often associated with red, brown, or purple shell coloration (Comfort, 1951). Investigation via modern chemical and multimodal spectroscopic techniques of shell pigments in marine snails Clanculus pharaonius and Clanculus margaritarius showed that pink-red dots and lines in shells were caused by porphyrins (Williams et al., 2016). Furthermore, carotenoids were regarded as the principal substance that presented orange adductor muscle and shell color in scallop since their first identification in the adductor muscle of rare orange Yesso scallop variants (Li et al., 2010). Subsequent analysis in reddish-orange and brown shell-color Yesso scallops identified candidate genes involved in carotenoid metabolism during shell color formation . A similar result was also reported in noble scallop populations, in which individuals with orange shell color contained dramatically higher carotenoid content compared with the brown ones (Zheng et al., 2010).
The bay scallop, Argopecten irradians irradians, is polymorphic for shell color (purple, orange, yellow, or white) and possesses advantages of fast growth and high production in aquaculture (Shumway and Parsons, 2006;Wang et al., 2007). Selective breeding of bay scallops for better growth and tolerance performance has been consecutively conducted for generations from the 1980s in China (Zheng et al., 2004(Zheng et al., , 2011Wang et al., 2017). Specifically, a potential phenotypic correlation between shell color and growth traits has been observed in artificially bred bay scallops (Wang et al., 2020). Consequently, understanding the genetic mechanism of shell color in bay scallops is of theoretical and economic importance. The first comprehensive classification of shell color based on experimental crosses in bay scallop showed that three background colors, six pattern colors, and the diversity of shell colors were under genetic control (Kraeuter et al., 1984;Elek and Adamkewicz, 1990). Subsequent study of shell color inheritance in selfing families found that orange or yellow parents produced both colorful and white progeny with a ratio of ∼3:1, implying that shell color in bay scallops was controlled by one locus being dominant. Furthermore, one amplified fragment length polymorphism (AFLP) marker was successfully identified in orange and white shell color bay scallops through genetic linkage analysis (Qin et al., 2007). Considering the large polymorphism of shell color in bay scallop, further genome-wide researches is still needed to identify the comprehensive candidate loci/genes and to explore the potential molecular mechanism.
Genome-wide association study (GWAS) is an accurate and powerful tool to identify novel genetic variants loci underlying complex traits and has been broadly applied in various fields (Zhao et al., 2011;Kim et al., 2012;Sukumaran and Yu, 2014;Ikeda et al., 2018;Ning et al., 2019). Investigations using GWAS have already successfully deciphered genetic variants of coloration traits, in which the identified specific markers could be potentially beneficial for selective breeding. For instance, a strong candidate locus, in arabidopsis pseudo-response regulator 2-like (APRR2) that regulates the accumulation of green pigments, was identified in melon and watermelon via GWAS (Oren et al., 2019). Information of two colors of single-nucleotide polymorphisms (SNPs) in French Saanen goats could be used to remove the undesired colored goats through molecular marker-assisted breeding (Martin et al., 2016). Furthermore, five potential loci related to carotenoid biosynthesis were identified in maize (Zeamays L.); therefore, a greater abundance of carotenoids accompanied by the deeper orange color of the kernel was potentially obtained (Owens et al., 2019). A similar GWAS result has also been concluded in Yesso scallop, in which three candidate genes (i.e., LDLR, FRIS, and FRIY) were detected around the two most significant potential loci and were proven to be involved in carotenoid metabolism in reddishorange shell individuals . These GWAS results not only provide a list of candidate loci/genes to be functionally validated but also offer an efficient approach for deciphering the genetic architecture of complex traits and selecting the breeding strategy.
This study aimed to identify significant SNP makers associated with shell color (purple-red and black-brown) in two new strains of bay scallops by performing GWAS, with the aid of whole-genome sequence databases of A. irradians irradians (unpublished data). The top three SNPs that dominantly control the shell color formation were further analyzed in their located candidate genes and then verified in another population. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to reveal the potential molecular mechanism in bay scallop shell coloration. Finally, we investigated the expression levels of the top three candidate genes in bay scallop mantles to detect their vital roles in shell coloration and to provide valuable information for the efficient selective breeding of shell color strains in bay scallop.

Scallop Materials and Miso-RAD Library Construction
Two strains of healthy bay scallops with nearly 9-month-old, namely, the purple-red shell color population (PP) (N > 500) and the black-brown shell color population (BP) (N > 500), were used as samples in 2016. The shell color trait of both strains is genetically stable due to the genomic selection of four consecutive generations from 2012 to 2015 by our group (Figure 1). Samples were collected from artificial scallop-rearing substrates at Qingdao Jinshatan Fishery Group Co. Qingdao (35 • 97 ′ 23 N, 120 • 27 ′ 28 E, Shandong Province). Subsequently, we randomly selected scallop samples (N = 80) from each group based on their shell color. Gills were removed from each scallop according to the procedure of Mao et al. (2013) and then preserved in 100% ethanol at −20 • C before genetic analysis. Our experiments were conducted according to the guidelines and regulations established by the Ocean University of China and the local government.
Genomic DNA was extracted from the gill using a standard phenol-trichloromethane method (Sambrook et al., 1989). Subsequently, Miso-RAD libraries were constructed by using 80 PP and 80 BP individuals, following the protocol developed by Wang et al. (2016). On summarizing, genomic DNA of the 160 selected bay scallop individuals was digested with BsaXI (New England BioLabs, MA, USA, Cat. no. R0609) at 37 • C for 45 min. Adaptor ligation, PCR amplification, digestion, and ligation were followed and after which every five DNA fragments were set as one group (barcoding and pooling). MinElute PCR Purification Kit (Qiagen, Hilden, Germany, Cat. no. 28004) was used to purify each product that concatenated five tags from five samples after performing amplification for two times. The 32 libraries were then pooled to run double-end sequencing on the Illumina Hiseq2000 system (San Diego, USA).

Sequence Data Processing and Genotyping
Raw reads were first preprocessed to eliminate paired-end (PE) reads with ambiguous base calls (N), long homopolymer regions (>10 bp), or excessive low-quality bases (>20% of bases with quality score <10) to obtain high-quality reads. PEAR software  was exerted to assemble the forward and reverse reads. The merged reads containing five tags were divided into single-tag datasets by using homemade Perl script, and via SOAP software (Li et al., 2009), each tag was subsequently aligned to the chromosome-level reference bay scallop genome (unpublished data) that was constructed by our group. The RADtyping software (Fu et al., 2013) was employed in the genotyping analysis of sequencing data. The selective rules of markers included polymorphism, sufficient genotype rate (with available genotype over 80% progenies), and minor allele frequency (MAF) higher than 0.05. Then, the qualified markers were used for downstream analysis.

Genome-Wide Association Study and Candidate Gene Searching
To reduce the false-positive rate, the population genomic structure analysis was performed to provide a population membership matrix through Admixture software 1.3 (Alexander et al., 2009), and a kinship matrix was estimated by TASSEL 5.0 (Bradbury et al., 2007). Association analysis between genotypes and shell colors was conducted using PLINK software 1.9  by employing the logistic regression model with these two covariates to detect the associate SNPs. The statistic pvalue for each SNP was calculated and summarized in ascending order. The threshold p-value for genome-wide significance was calculated using the Bonferroni correction based on the number of qualified makers (Johnson et al., 2010). When an SNP scored less than the significance cutoff (1.72E−6), the sequence of each maker was extracted and used as a query to exert BLASTN (Altschul et al., 1997) to search the draft genome of A. irradians irradians (unpublished data) that were annotated by the protein database NR. In addition, the Manhattan plot was drawn by an R package named "ggplot, " and these potential SNPs (p < 1.00E−5) were captured to conduct bioinformatics analysis, including GO function annotation and KEGG pathway enrichment analyses (http://www.omicshare.com/tools).

SNP Verification
To further evaluate the accuracy of GWAS, the top three SNPs (p-value < 1.00E−17) were selected to verify in the 30 purplered shell color samples and 30 black-brown shell color samples collected from another separate bay scallop population (N > 500) in the sea area of Huangdao District, Qingdao, Shandong Province, in 2018. Striated muscles of the above 60 individuals were separately sampled and preserved at −80 • C for subsequent analysis. The primer sequences are provided in Table 1. The PCR of fragments containing the above three SNPs was carried out as follows: 2 µl of template DNA (50 ng/µl), 15 µl of Q5 High-Fidelity 2× Master Mix (New England BioLabs, MA, USA, Cat. no. M0492S), 0.5 µl (10 pmol/µl) of each forward and reverse primer, and sterile water up to 30 µl volume was added for amplifying. The related parameters in the PCR system for three primers were pre-denaturing under 98 • C for 30 s, then denaturing under 98 • C for 10 s, annealing under 56 • C for 15 s, and extending under 72 • C for 10 s, respectively. After 28 cycles, the sample was extended under 72 • C for 2 min. The PCR products were sequenced using the Sanger method by Sangon Biotech (Shanghai, China). Sequences amplified by the primers were compared among the 60 scallops (30 PP and 30 BP, respectively) using the ClustalW2 multiple alignment program (http://www.ebi.ac.uk/Tools/msa/clustalw2/). When a mutation was detected, comparisons of genotype frequencies between the PP and BP were performed using Fisher's exact test to identify the shell color-associated mutation. p-values < 0.05 were considered statistically significant.

Candidate Gene Evaluation
To compare the differences in the expression levels of genes which are the top three SNPs annotated between PP and BP, 30 PP and 30 BP individuals from another non-breeding bay scallop population (N > 500) were randomly recruited, and the mantle was sampled in the sea area of Huangdao District, Qingdao, Shandong Province, in 2020. Total RNA of the sampled mantles was isolated following the method described by Hu et al. (2006) and then digested with DNaseI (Takara, Shiga, Japan). The concentration and purity of the RNA were determined using a Nanovue Plus spectrophotometer (GE Healthcare, Piscataway, NJ, USA), and the RNA integrity was assessed by agarose

SNP/Gene
Primers sequence Note gel electrophoresis. The first-strand cDNA was synthesized according to the protocol of the manufacturer using Moloney murine leukemia virus (MMLV) reverse transcriptase (Thermo, Wilmington, USA) in a 20 µl volume with 2 µg of each total RNA sample as the template and 0.5 µg of oligo (dT)18 (Takara Biotechnology, Liaoning, China) as the primer. The mixture was denatured at 65 • C for 5 min and then chilled immediately on ice. After adding the reverse transcriptase, reaction buffer, and dNTPs, cDNA was amplified under the following conditions: 42 • C for 90 min and 72 • C 10 min. The cDNA was stored at −20 • C and diluted to 5 ng/µl for use as the template in the real-time quantitative reverse transcription PCR (qRT-PCR). Data from the qRT-PCR were obtained and the gene expression level is shown as the fold change. The statistical analysis of the data was performed with SPSS software (version 21.0) via an independent sample t-test. Differences were considered significant at p < 0.05.

Sequencing Data
In total, sequencing of 32

Genomic Regions Associated with Shell Color
A total of 211,728 SNP markers obtained from 160 bay scallop individuals were exerted to quality control performed by PLINK software. After quality filtering (call rate > 80% and MAF > 5%), a set of 29,036 qualified SNPs were successfully assigned to the 16 chromosomes of bay scallop and used for  further genome-wide analysis (Figure 2). Genetic clustering illustrated that the admixed ancestry of PP and BP was shown for K = 2, with kinship coefficients varying from 0 to 0.02 between PP and BP individuals, then the corresponding components were brought as covariates. After conducting Bonferroni correction and 0.05 cutoff (p < 1.72E−6), we detected 469 SNPs significantly associated with shell color, and 285 of these 27 bp tags could be directly annotated in the bay scallop genome. Furthermore, the -log 10 (p-value) values of the top eight SNPs with a significant association of shell color were higher than 15, and these top eight SNPs were mainly distributed in chromosome 7 (Table 2). Specifically, the top eight loci were separately located at the genic regions of GRL101 (G-protein-coupled receptor-like 101), PKS1 (polyketide synthase 1), PLC1 (phosphoinositide phospholipase C), GLRA1 (glycine receptor subunit alpha-1-like), COL4 (zinc finger protein constants-like 4), CL1 (cathepsin L1-like), and gene LOC110449017 (not yet characterized), as well as at an intergenic region (where SNP Chr13-253787 anchored). Among these, the top three loci were distributed in the genic region of chromosome 7 (Figure 3). More specifically, SNP chr7-12764003 was mapped to location 41 bp in the coding sequence (CDS) 7 of GRL101, and SNP chr7-13213864 was cited in position 1,347 bp in CDS 5 of PKS1. For chr7-11899306, an SNP was located at position 14,136 bp in the first intron of PCL1. Subsequently, the top three SNPs were selected to be verified by sequencing the PCR products using the Sanger method by Sangon Biotech (Shanghai, China).

SNP Verification
All the top three SNPs provided by the GWAS results were successfully verified by Sanger sequencing in another separate bay scallop population. The allele and genotype frequency of these three SNP markers are shown in Table 3. For SNP chr7-12764003, dimorphic alleles could be observed, with AA and GG presenting in PP and BP, respectively. As for SNP chr7-13213864, the genotype of AA (100%) was only observed in PP, while another two genotypes of TT and AT were present in BP, with frequencies of 76.67 and 23.33%, respectively. Besides, the genotype is varied in SNP chr7-11899306 (for CC, TT, CT; 76.67, 0, 23.33% in PP and 6.67, 63.33, 30% in BP, respectively). Comparison between PP and BP individuals showed that genotype frequencies of the above three SNP markers were significantly different (p < 0.05). These results support that the verified SNPs

GO and KEGG Analyses of Candidate Genes
After a widely used and less-stringent threshold p-value < 1.00E−5 (Marigorta et al., 2018), we obtained 713 significant SNPs annotated in 374 unigenes. Subsequently, GO function annotation and KEGG pathway enrichment were applied to excavate the potential association between these genes and shell color. Results displayed that 38 GO terms were classified into the biological process (BP), cellular component (CC), and molecular function (MF) under level 2 covering 297 genes (Figure 4). In detail, genes related to the cellular process (105 genes), single-organism process (94 genes), and metabolic process (68 genes) showed remarkable terms in BP. Genes related to the membrane (53 genes) displayed the notable term in FIGURE 4 | The summary of Gene Ontology (GO) function annotation analysis for candidate genes that comprised significant SNPs (p < 1.00E−5). Each term describes the function of the gene cluster, and the length of colored bars represents the difference of the number of genes.
CC. Genes associated with binding (166 genes) and catalytic activity (80 genes) were highly presented in MF (Figure 4). In terms of KEGG pathway enrichment analysis, a total of 252 genes were acquired annotations and involved in 44 pathways ( Figure 5). The above candidate genes were chiefly distributed in metabolism (69 genes), genetic information processing (15 genes), environmental information process (55 genes), cellular process (38 genes), organismal systems (90 genes), and human disease (97 genes) ( Figure 5). Specifically, signal transduction (37 genes), global and overview maps (29 genes), infectious diseases (23 genes), cancers (21 genes), and endocrine system (19 genes) were mainly presented. Overall, the result indicated that the extracted GO terms and KEGG pathway are closely related to cellular process, membrane, signaling, and binding, and these possible molecular mechanisms may attribute to the shell color formation.

Expression Patterns of the Top Three Genes
Annotation results showed that the top three SNPs were located in three genes, including GRL101 (chr7-12764003), PKS1 (chr7-13213864), and PLC1 (chr7-11899306), which were further selected for comparing their expression differences in mantles between PP and BP individuals. As shown in Figure 6, the expression profiles of the three target genes were ubiquitous in both PP and BP individuals but exhibited distinct expression patterns. Specifically, an extremely significant higher (6.43fold, p < 0.01) expression level of GRL101 was detected in PP than that in BP. A similar phenomenon was also observed in PLC1, which displayed a prominently increasing expression level in PP (6.48-fold, p < 0.05) compared with BP. On the contrary, the expression of PSK1 was notably correlated with black-brown shell color rather than purplered shell color (12.02-fold of BP compared with PP). Overall, in the bay scallop, the significantly higher expression profiles of GRL101 and PLC1 were detected in PP, implying that they may play a pivotal role in purple-red shell formation. Besides, PSK1 is principally expressed in BP instead of PP, suggesting that it may involve black-brown pigmentation during shell formation.

DISCUSSION
Aquaculture breeding researches have confirmed that color trait in most aquaculture species is heritable, including red/white flesh color in salmon (Iwamoto et al., 1990), carapace/hepatopancreas color in crab , black/white shell color in oyster (Xu et al., 2019), white/orange in scallop adductor muscle , etc. Particularly, the bay scallop is polymorphic for the background color of the shell and presents a highly variable distribution of overlying pigments (Clarke, 1965). Although previous investigations showed one or several major genes contributed to shell color diversity of bay scallop FIGURE 5 | The summary of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis performed for biological process, cellular component, and molecular function. Functional annotation of candidate selective genes that comprised significant SNPs (p < 1.00E−5). (Adamkewicz and Castagna, 1988;Qin et al., 2007;Teng et al., 2018), the precise loci/genes for specific color difference and the mechanism underneath shell color formation in bay scallop are elusive. With rapid development in biotechnology, constantly upgraded methods for studying shell color polymorphisms in bay scallop would dredge out specific genes behind target trait. GWAS has been proven to be an efficient approach to elucidate the genetic variations related to complex phenotypes and to identify the pivotal candidate genes of target economic traits, thereby providing referential makers for the selective breeding Frontiers in Marine Science | www.frontiersin.org FIGURE 6 | The relative expression levels of GRL101, PKS1, and PLC1 in the mantle of PP and BP bay scallops. Three replicates were performed for each adult mantle, and three technical replicates were conducted for each PCR. The comparison of the expression levels of GRL101, PKS1, and PLC1 in the mantle of PP and BP bay scallops was performed using an independent sample t-test. * and ** indicate statistically significant differences (p < 0.05 and p < 0.01, respectively).
projects (Luo et al., 2012). By this method, 469 SNPs that displayed statistically significant differences between PP and BP were identified in our research. As the Manhattan plot displayed, the top value of -log 10 (p-value) reached nearly 26, which is much more stringent than 5.8 (p < 1.72E−06) that referred to Bonferroni correction, demonstrating that our result is reliable (false positives are rare). Intriguingly, most of the significant SNPs associated with shell color were primarily allocated in chromosome 7 of bay scallop. With the similarity to that in Yesso scallop, the major SNPs associated with shell color were located in chromosome 11 , and almost all the significant SNPs associated with carotenoid coloration in scallop adductor muscle were distributed in chromosome 8 . Based on the fact that most qualitative traits are typically determined by one or several pivotal loci/genes, we speculated that color variants in scallop were potentially regulated by limited candidate loci/genes gathered closely in association with the genomic region(s) of the specific linkage group. In addition, the tandem or proximal (in a nearby chromosomal region but not adjacent) locations may not only facilitate coordinate expressions of candidate genes but also guarantee efficient expression regulations by pivotal loci. The top three SNPs with a p-value <1.00E−17 were successfully verified in another separate population via Sanger sequencing. The extremely significant (p < 0.01) or significant (p < 0.05) differences of the three selected loci between PP and BP suggested that these SNPs might be responsible for shell color formation and potentially used as selective breeding markers in shell color-specific projects of bay scallop. In addition, the top eight SNPs (p < 1.00E−15) were principally and separately located in genes (exon or intron) except on locus (SNP chr13-253787) located in intergenic regions. A similar phenomenon has also been concluded in the shell color of the Yesso scallop . The above results revealed that color traits in scallops may be principally controlled by the coding regions (where SNPs are located) of most critical candidate genes (Cargill et al., 1999) and meanwhile regulated through crucial SNPs within introns or intergenic regions via transcriptional regulation (Sturm et al., 2008;Visser et al., 2014).
Bay scallop shell color significance-associated candidate genes, namely, GRL101, PKS1, and PLC1, where the top three SNPs were located, were further selected for comparing their expression differences between PP and BP individuals. The scallop mantle constantly contacts with the shell; therefore, the mantle is considered to manufacture and transport pigments to the shell along the growing edge during shell formation (Fowler et al., 1992). Upon this basis, further expression patterns of the top three candidate genes were characterized in the mantle of PP and BP individuals, respectively. In this study, the mRNA expression level of GRL101 in PP was over six times that in BP (p < 0.05), suggesting GRL101 might function in purple-red shell coloration. As reported, GRL101 belongs to the G-protein-coupled receptors (GPCRs) superfamily (Karnik et al., 2003) that exerts multifunctions, such as a variety of hormones and neurotransmitters regulation (Nilaweera et al., 2006(Nilaweera et al., , 2008Rosenbaum et al., 2009). Specifically, a set number of GPCRs were proven to be involved in pigments formation, known as responding to red-pigment concentrating hormone (RPCH) in Carcinus maenas (Alexander et al., 2018) and interceding pigment granule dispersion in Xenopus laevis melanophores (Teh and Sugden, 2001). Furthermore, GPCRs activated by chromatophores can affect pigment redistribution to alter the color appearance of the animal (Lerner, 1994). In particular, visual pigment rhodopsin (the first order of GPCRs) performs functions on the absorption and conversion of light into chemical signals (Smith, 2010). Based on this, our results implied that GRL101 involved and functioned in the purple-red color shell formation of bay scallop, whereas GRL101 protein was regarded as an orphan receptor performing an uncertain function that was extensively explored in humans (Gene ID: 83550), mice (Gene ID: 245424), sea urchins (Gene ID: 115922752), king scallops (Gene ID: 117318188), and other species. Whether the GRL101 in the bay scallop is also considered an orphan receptor still needs to be further investigated. Coincidently, the majority of GPCRs were coupled to the phosphoinositide signaling pathway where the third candidate gene PLC1 played a pivotal role (Tobin, 1997). Previous researches reported that PLC1 has been classified into an important class of enzymes involved in lipids-related signaling in plants (Rupwate and Rajasekharan, 2012) and acted as signal transducers to generate messengers participating in lipid metabolism in mammals (Essen et al., 1996). Furthermore, the lipid-related genes were proved to be closely associated with carotenoids (Zhang Y. et al., 2014) that contributed to the formation of reddish-orange shells in Yesso scallop . Hence, it might then make sense that much higher expression of PLC1 detected in PP individuals leads to the accumulation of carotenoid which in turn resulted in a red-colored shell. Therefore, our results suggest that PLC1 might be related to lipid-associated carotenoid absorption and/or transportation. Collectively, we speculated that GRL101 might perform a common feature of coloration-related GPCRs and activate PLC1 in the phosphoinositide signaling pathway to affect lipid metabolism, including lipid-associated carotenoid regulation, resulting in the purple-red visualization of PP.
Another candidate gene PKS1 belongs to polyketide synthase-related genes (Moriwaki et al., 2004), and the ancient and diversified family members have been well studied in animal fatty acid synthase (Kroken et al., 2003). As reported, plenty of microorganisms employed PKS to produce pigments and other products of intermediate metabolism (Hutchinson, 2003). For example, the biosynthesis of black pigment melanin derived by PKS1 was reported in Colletotrichum lagenarium (Takano et al., 1995) and Aspergillus fumigatus (Heinekamp et al., 2013), respectively. Unlike PKS1 as the pivotal component in melanin biosynthesis in microorganisms, melanin biosynthesis in various tissues of vertebrates was prevailingly affected by tyrosinase genes (Sato et al., 2001). Similarly, melanin biosynthesis was more inclined to be regulated by tyrosinase that contributes to black-brown shell pigmentation in most mollusks (Bai et al., 2013;Sun et al., 2015). However, our study showed that a higher expression level (12.02-fold, p < 0.05) of PKS1 in the black-brown shell individuals was observed than that in purple-red ones; therefore, we speculated that black pigmentation in BP might be regulated by PKS1 as it functioned in microorganisms. Although melanin synthesis in scallop has not been reported yet, the typical biosynthesis pathway (in which tyrosinase is involved) as an alternative method seems worth further study.
Generally, causal SNPs were identified to be strongly associated with target traits, especially the SNPs that were located in the coding region would affect traits directly by changing their encoding amino acid and therefore resulting in the alteration of protein (Seong et al., 2011). However, in this study, the two most significant SNPs, namely, chr7-12764003, and chr7-13213864, were separately located in the CDS of AiGRL101 and AiPSK1, caused synonymous mutations in target amino acids. We speculated that SNPs located in the CDS of AiGRL101 and AiPSK1 might not control the scallop shell coloration directly via mutation of amino acids, but mediate regulating candidate genes by affecting codon bias (Plotkin and Kudla, 2011), translation efficiency (Miyuki and Kentaro, 2013), mRNA secondary structure stability (Chamary and Hurst, 2005), and splicing . Furthermore, linkage disequilibrium (LD) analysis of the three closely arranged SNPs (Figure 2) via MEGA (Version 10.1.8) (Sudhir et al., 2018) showed that LD coefficient r 2 between each pair of SNPs was higher than 0.7 (ranged from 0.72 to 0.78), demonstrating that shell color in bay scallop may be regulated by closely linked major SNPs and/or specific interactions among them which potentially affect the expression of candidate genes resulting in shell color formation and differentiation. Further research is needed for investigating the potential molecular mechanism of synonymous mutation of CDS in candidate genes and interactions of target genes in regulating shell coloration of bay scallops.
Furthermore, 713 SNPs (p < 1.00E−5) significantly associated with bay scallop shell color were mapped to 374 unigenes, followed by GO function annotation and KEGG pathway enrichment for elucidating their important molecular mechanism . GO function annotation in this study primarily participated in cellular process, membrane, and binding function, which is similar to previous research that cellular process and binding are the largest group of biological process and molecular function in the transcriptomic analysis of red/white-colored valves of moon scallop (Huang et al., 2015), as well as different shell color lines in Yesso scallop (Ding et al., 2015). Furthermore, transcriptome sequencing of orange shell color bay scallops revealed that differentially expressed genes were mainly annotated in signal transduction, metabolism, human diseases, and organismal systems (Teng et al., 2018) and that is in accordance with our result of KEGG pathway analysis. Taken together, the above results suggest that the shell color of the bay scallop is not only dominated by several loci/genes but also modified by plenty of related minor loci/genes.

CONCLUSION
We identified SNPs associated with the purple-red and black-brown shell colors bay scallop through GWAS. The genomic regions associated with shell color were mainly located in chromosome 7. The top three SNPs were located in the genic region of GRL101, PKS1, and PLC1 and successfully verified in another separate population. Expression patterns analyses suggested that GRL101 and PLC1 potentially functioned in purple-red shell coloration, and black pigmentation in BP might be principally regulated by PKS1. GO function annotation and KEGG pathway enrichment analyses revealed the possible molecular mechanism of shell color formation. Taken together, this evidence collectively provides valuable information for deciphering the phenotype differences of shell color and helps in the molecular marker-assisted breeding in bay scallops.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GenBank (accession number: PRJNA689862).

ETHICS STATEMENT
The animal study was reviewed and approved by Ocean University of China.

AUTHOR CONTRIBUTIONS
QX, JH, and ZB conceived and designed the experiments. XZ, JZ, and XHo collected the samples. QX, XZ, PL, and JL performed the experiments. XZ, QX, JZ, and XHu analyzed the data. QX and XZ wrote the manuscript. All authors have read and approved the final manuscript. All authors contributed to the article and approved the submitted version.