141 ISJ 16: 141-151, 2019 ISSN 1824-307X RESEARCH REPORT Identification of BMP2 and BMP7 genes and association of their SNPs with growth traits in the hard clam (Meretrix meretrix) G Pan1, 2, X Gao1, Z Lin1, H Yao1, Y Dong1* 1Key Laboratory of Aquatic Germplasm Resources of Zhejiang, Zhejiang Wanli University, Ningbo, 315100 2National Demonstration Center for Experimental Fisheries Science Education, Shanghai Ocean University, Shanghai, 201306 Accepted July 18, 2019 Abstract The hard clam (Meretrix meretrix) is a commercially and ecologically important benthic mollusk in South and Southeast Asia. Growth is one of the most important and complex traits affecting clam breeding. Therefore, identification of key growth factor genes and their association with growth traits are the primary prerequisites for in-depth study of the mechanisms of growth regulation and molecular marker-assisted selection of the hard clam. In the present study, we generated the full-length cDNA sequences of the growth-related genes BMP2 and BMP7 in M. meretrix (MmBMP2 and MmBMP7) and then investigated their gene expression patterns and associations with growth traits. The full-length cDNA of MmBMP2 was 2641 base pairs (bp), with a 1395 bp open reading frame (ORF) encoding 464 amino acids. The full-length MmBMP7 cDNA was 3161 bp, with a 1257 bp ORF encoding 418 amino acids. In adult clams, expression of MmBMP2 and MmBMP7 were significantly higher in the mantle than in the other five tested tissues (p < 0.05). Among the seven developmental stages of M. meretrix, expression of MmBMP2 was the highest in juveniles (p < 0.05), and expression of MmBMP7 was very low in 4-cell embryo and eyespot larva stages (p < 0.05). Two hundred and sixty-three individuals from two populations were screened for MmBMP2 and MmBMP7 polymorphisms, and single nucleotide polymorphisms (SNPs) associated with growth traits in MmBMP2 was found in both populations. Subsequent analysis revealed that two SNPs in MmBMP2 were shared between the populations. One was A726G, for which individuals with genotypes AG had significantly higher growth traits than those with genotype AA (p < 0.05). The other was A732T, for which individuals with genotype AT had significantly higher growth traits than those with genotype AA (p < 0.05). These results suggest that BMP2 markers may have potential value in molecular selective breeding to improve clam growth. Key Words: Meretrix meretrix; BMP2; BMP7; SNP; growth traits Introduction The transforming growth factor-beta (TGF-β) supergene family plays various roles in the regulation of cell growth, cell specialization, and morphogenesis (Massague, 1998; Derynck and Akhurst, 2007). More than 40 members of TGF-β have been identified in mammals (Li et al., 2011), including bone morphogenetic proteins (BMPs), Smads, and activin. BMPs are multi-functional growth factors, and their major function is to induce the formation of both bone and cartilage in vertebrates (i.e. the process of biomineralization) ___________________________________________________________________________ Corresponding author: Yinghui Dong College of Biological & Environmental Sciences Zhejiang W anli University 8 South Qianhu Rd, Ningbo, Zhejiang 315100, P. R. China E-mail: dongyinghui118@126.com (Chen et al., 2004; Xiao et al., 2007). Thus, most bone-related molecular studies focus on BMP2, BMP4, and BMP7 (Salazar et al., 2016). BMPs also play a role in cell growth and differentiation, embryogenesis, hematopoiesis, and neurogenesis (Chen et al., 2004; Bragdon et al., 2011). BMPs have been found in invertebrates, including echinoderms (Hwang et al., 1994; Shih et al., 2002), arthropods (Wharton et al., 1991; Fritsch et al., 2010), mollusks (Nederbragt et al., 2002; Iijima et al., 2008; Miyashita et al., 2008), platyhelminths (Freitas et al., 2009), cnidarians (Reinhardt et al., 2004; Zoccola et al., 2009), and poriferans (Müller et al., 2011). Among the members of the BMPs family, BMP2 is the factor that can induce bone formation alone, inhibit myogenic differentiation, and promote skeletal cell differentiation (Musgrave et al., 2001). During mailto:dongyinghui118@126.com 142 Table 1 Primers used for MmBMP2 and MmBMP7 genes in the experiments embryonic development, BMP2 not only promotes embryonic muscle growth by inducing the expression of myogenic genes, but also inhibits muscle growth by inducing apoptosis (Amthor et al., 1998). BMP7 is a key factor in bone homeostasis and spinal development (Cook and Rueger, 1996). Previous studies indicated that BMP7 is abundantly expressed in articular cartilage and synovial membrane obtained from human patients undergoing autologous chondrocyte implantation (Schmal et al., 2012). Furthermore, researchers discovered that BMP7 in mouse was expressed in a wide variety of tissues during development (Lyons et al., 1995) and was strongly expressed in the developing myocardium of embryos (Dudley and Robertson, 2010). BMP7-deficient mice were found to exhibit defects in eye development, renal development, and skeletal patterning. The deletion of BMP7 gene results in lens development defects, affects the development of ribs and skulls, and causes polydactyly of hind limbs. (Dudley et al., 1995; Luo et al., 1995). Because of their important roles in growth and development, BMP2 and BMP7 have become widely studied genes in mollusks. Several researchers reported that BMPs and their homologs are involved in the regulation of morphogenesis in the shells of mollusks (Herpin et al., 2010; Lelong et al., 2010). Furthermore, dpp-BMP2/4 is expressed around the shell gland of the embryo in gastropods such as Lymnaea stagnalis and Patella vulgata, which suggests that it may be involved in determination of the boundary of the embryonic shell region and in shell formation (Nederbragt et al., 2002). The high-level expression of Pinctada fucata BMP2 in the inner layer of mantle tissue suggests that it plays an important role in the formation of nacre (Miyashita et al., 2008). In Pinctada martensii, the microstructure observed using a scanning electron microscope indicated a disordered growth status in the nacre and obvious holes in the prismatic layer in the dsRNA-Pm-BMP7-injected group, which suggests that BMP7 plays a crucial role in the nacre and prismatic layer formation process of the shell (Yan et al., 2014). The hard clam (Meretrix meretrix) is a commercially important bivalve species that is distributed in the intertidal zone of South and Primer name Sequence from 5’to 3’ Application MmBMP2-F1 ATGGCGTGTCCTCTCAGTTTAC 3’ RACE cloning MmBMP2-R1 ATCACTCCCCAGCCTGTCCTTCAATG 5’ RACE cloning MmBMP2-F2 GTGGGCAGGATTTCACTGTAGG Verifying the sequence of cDNA MmBMP2-R2 GTAGCATTGAAGGACAGGCTGGGGAGT Verifying the sequence of cDNA MmBMP2-F3 CTAACAGCAGCGACAGTGTGT SNP MmBMP2-R3 CTGGACTTCATCTTCGTGGAA SNP MmBMP2-F4 CCAGTCTACCTACAGAACAGCC SNP MmBMP2-R4 TCCGTTTAGTCCGTGATAAGC SNP MmBMP2-F5 TATCACGGACTAAACGGAGCA SNP MmBMP2-R5 CGTAAACTACCTTGACTGGAAGCA SNP Real-2-F TCAGTCGGGTGTATTTTAGCAGA Real time PCR Real-2-R CATTTTTTCCAGCAGGTCAAGTG Real time PCR MmBMP7-F1 GTCAGGTGAGCGATGTAGGGTCGTTTC 3’ RACE cloning MmBMP7-R1 CAAGTGTCTGCTGGTTCCCGATTCC 5’ RACE cloning MmBMP7-F2 AAGTTGTGATAAAACTGCGCG Verifying the sequence of cDNA MmBMP7-R2 AACTATATATAGCTTGCATC Verifying the sequence of cDNA MmBMP7-F3 TTTAGGCAATTTCCCGAGAC SNP MmBMP7-R3 CCTTGTCGTGCCTAAGATGC SNP MmBMP7-F4 TTCCGCATCTTAGGCACGAC SNP MmBMP7-R4 GGGTAACCTTCAGGGGCGAT SNP MmBMP7-F5 ATTATCGCCCCTGAAGGTTAC SNP MmBMP7-R5 AAACGACCCTACATCGCTCAC SNP Real-7-F GTTTTCGCCTGTTGTGGTG Real time PCR Real-7-R GAGTGTGCGGTGTTCGTGTA Real time PCR 18S-F CTTTCAAATGTCTGCCCTATCAACT Real time PCR 18S-R TCCCGTATTGTTATTTTTCGTCACT Real time PCR 143 Southeast Asia (Liu et al., 2006; Tang et al., 2006). With the success of artificial breeding, clam culture has become an economically important aquaculture industry in China. Although a selective breeding program for genetic improvement of growth traits of M. meretrix has begun (Wang et al., 2011), association analysis between phenotype and genotype has still been relatively rare to date. However, molecular marker-assisted selection of growth traits would probably further increase selection efficiency and lead to rapid progress in improving the process of clam breeding (Lande and Thompson, 1990; Hospital et al., 1997; Moreau et al., 1998; Xie and Xu, 2010). Therefore, identification of growth-related markers for further development of phenotype-selective breeding is very necessary. Studies of growth-related functional genes of M. meretrix, including growth factor receptor-bound protein 2 (Gao et al., 2015), cathepsin B (Yao et al., 2011), sulfotransferase (Wang et al., 2012), and fructose-1, 6-diphosphate aldolase 9 (Wang et al., 2012), have been conducted. However, the M. meretrix BMP2 and BMP7 genes (MmBMP2 and MmBMP7) have not been studied to date. In the current study, we investigated MmBMP2 and MmBMP7 by cloning their full-length cDNAs and examining their expression levels during embryonic development and in adult tissues. We also screened and analyzed associations between single nucleotide polymorphisms (SNPs) and growth traits. The results of this study will improve our understanding of the functions of MmBMP2 and MmBMP7 and provide a basis for molecular breeding programs for M. meretrix. Materials and methods Experimental materials and sample collection For cloning and gene expression analysis of MmBMP2 and MmBMP7, adult clams were obtained from the Ningbo Danyan Aquaculture Company, Zhejiang, China. Six tissues (digestive gland, mantle, adductor muscle, foot, siphon, and gills) were dissected, immediately frozen in liquid nitrogen, and stored at -80 °C. Samples from seven developmental stages (unfertilized mature eggs, 4-cell embryos, blastulae, gastrulae, trochophores, eyespot larvae, and juvenile clams) were collected and preserved at -80 °C. A total of 263 individuals were randomly sampled at the same time from two populations: 136 individuals from the Jiangsu population (JS) and 127 individuals from Wanli No. 2 strain (WL2). JS was collected from the wild population located in Jiangsu province, China. WL2, a new fast-growing variety with a dark gray shell color and zigzag patterns, had been selected for five generations by our research team. The following six growth index parameters, total weight, soft weight, shell weight, shell width, shell length, and shell height, were measured for each clam. The mantles of each clam were collected and stored at -80 °C. Cloning the full-length cDNA of MmBMP2 and MmBMP7 Total RNA was extracted from the mantle using TRIzol reagent (Sangon, Shanghai, China). RNA integrity was examined by electrophoresis on a formaldehyde-denatured 1.2% agarose gel and staining with ethidium bromide, while the quality and quantity were assessed by ultraviolet spectrophotometry. First-strand cDNA was synthesized using the simple modular architecture research tool (SMART) rapid amplification of cDNA ends (RACE) reagents (Clontech, USA) according to the manufacturer’s instructions. Based on the 454 cDNA library of M. meretrix (GenBank accession no. SRA021052), the expressed sequence tags of BMP2 and BMP7 genes were retrieved. The primers for 3’-RACE (MmBMP2-F1, MmBMP7-F1) and 5’-RACE (MmBMP2-R1, MmBMP7-R1) were designed (Table 1). The amplification was conducted as follows: 5 cycles of 94 °C for 30 s, 72 °C for 3 min; 5 cycles of 94 °C for 30 s, 70 °C for 30 s, and 72 °C for 3 min; 25 cycles of 94 °C for 30 s, 68 °C for 30 s, and a final extension at 72 °C for 3 min. The PCR products were purified using a gel extraction kit (Tiangen, Sichuan, China). The purified PCR products were cloned into the pMD18-T vector (Takara), which was transformed into DH5α cells (Tiangen) according to the manufacturer’s protocols, and the positive clones were sequenced. The full-length cDNA sequence was determined by piecing together the sequences of the 3’ and 5’ RACE products. To confirm the accuracy of cloning and sequencing, the full-length cDNA was re-amplified with high fidelity polymerase (Takara) using gene-specific primers (MmBMP2-F2 and MmBMP2-R2, MmBMP7-F2 and MmBMP7-R2) (Table 1) that were designed based on the above-mentioned cDNA sequences. PCR products were cloned and sequenced following the procedures described above. Sequence analysis The sequences were spliced using the BLAST algorithm in the National Center for Biotechnology Information database (http://www.ncbi.nlm.nih.gov/blast/). The deduced amino acid sequences were analyzed using the SMART tool (http://smart.embl-heidelberg.de) to predict conserved domains. The presence and locations of the signal peptide and cleavage sites in the amino-acid sequence were predicted using the Signal P 4.0 server (http://www.cbs.dtu.dk/services/SignalP/). Protein transmembrane were analyzed using the TMHMM server (http://www.cbs.dtu.dk/services/TMHMM). Multiple alignments of BMP2 and BMP7 domains structures from M. meretrix and other species were performed using the ClustalW2 multiple alignment program (http://www.ebi.ac.uk/Tools/msa/clustalw2/). A phylogenetic tree was constructed using the neighbor-joining method with MEGA 6.0. Quantitative expression analysis of MmBMP2 and MmBMP7 The expression levels of MmBMP2 and MmBMP7 at seven developmental stages (n > 500, three sets of samples for each stage) and in six tissues from adults (four sets of samples for each tissue) were analyzed using real-time quantitative reverse transcription PCR (qRT-PCR), with three 144 technical repeats for each PCR reaction. Total RNA was extracted from the samples using the method described above. The primers for MmBMP2 and MmBMP7, Real-2-F and Real-2-R, and Real-7-F and Real-7-R (Table 1) were designed. Primers 18S-F and 18S-R (Table 1) were used as the internal reference for qRT-PCR. PCR amplification was performed in a 20 μl reaction volume containing 10 μl of iTaq Universal SYBR Green Supermix (Bio-Rad, CA, USA), 7.2 μl of deionized water, 0.8 μl of the first strand cDNA, and 1 μl of forward and reverse primers. Amplification was performed using the following thermal cycling conditions: incubation at 94 °C for 20 s, 40 cycles of 94 °C for 3 s, 60 °C for 15 s, and 72 °C for 10 s. All qRT-PCR analyses were performed using three biological replicates, and negative controls were run with each primer pair. The results of qRT-PCR analysis were based on the CT values of the PCR products, and the expression levels of MmBMP2 and MmBMP7 were analyzed using the comparative CT method. Statistical analysis was performed using SPSS 20.0 (Chicago, IL, USA). Differences in relative MmBMP2 and MmBMP7 mRNA expression levels among different developmental stages and different adult tissues were compared by one-way analysis of variance (ANOVA). Multiple comparisons were conducted using the Student-Newman-Keuls test. Differences were considered significant at P < 0.05 and extremely significant at P < 0.01. SNP identification and association analysis RNA samples were extracted from all 263 clams. Based on the cDNA sequences of MmBMP2 and MmBMP7, six pairs of primers (2-F3, 2-R3, 2-F4, 2-R4, 2-F5, 2-R5, 7-F3, 7-R3, 7-F4, 7-R4, 7-F5 and 7-R5) were selected (Table 1). PCR amplification was performed in 50 µl reaction volumes containing 2 µl of template cDNA, 19 µl of dH2O, 25 µl of PCR Mix, and 2 µl of each primer. The PCR conditions were as follows: 94 °C for 3 min; 35 cycles at 94 °C for 30 s, 58 °C for 1 min, and 72 °C for 2 min; and a final extension at 72 °C for 5 min. Each amplification product was verified by electrophoresis on a 1.0% agarose gel and staining with ethidium bromide. The amplicons representing unique banding patterns were sent to Sangon Biotech (Shanghai, China) for sequencing in both directions. The sequences were aligned using MEGA6.0 software. Mutation sites were named according to the position of the initiation codon. Total weight, soft weight, shell width, shell length, shell height, and shell weight of the clams were analyzed by one-way ANOVA using SPSS 20.0. The association of SNPs with growth traits in the two populations was analyzed using the above method. SNP markers with genotypes that showed a significant correlation with growth traits were studied by post hoc multiple comparison (Duncan method). Results cDNA sequence analysis The full-length MmBMP2 cDNA was 2641 bp (GenBank accession no. KP250876) and contained a 1395-bp open reading frame (ORF) encoding 464 amino acids. The calculated molecular mass of the deduced mature MmBMP2 was 52.906 kDa, and the theoretical isoelectric point was 9.33. Analysis with SignalP revealed that the deduced amino acid sequence contained a signal peptide of 25 amino acids (MTPSYRVSILVLAVLLSELITSTYT). The full-length MmBMP7 cDNA was 3161 bp (GenBank accession no. KP250875) and contained a 1254-bp ORF encoding 418 amino acids. The calculated molecular mass of the deduced mature MmBMP7 was 47.54 kDa, and the theoretical isoelectric point was 7.81. Analysis with TMHMM showed that the MmBMP7 protein contained a potential transmembrane domain that was located between amino acid residues 12 and 31. SMART analysis showed that the deduced amino acid sequence of them both contained a TGF-β superfamily domain (Table 2). A phylogenetic tree was constructed using the neighbor-joining method and the deduced amino acid sequences of MmBMP2 and MmBMP7 and BMPs from some selected animals available in the NCBI database (Fig. 1). The tree obtained showed that BMP2 and BMP7 were divided into two major groups. One group contained all mollusks, whereas the other group contained vertebrates. In the BMP2 mollusk group, Azumapecten farreri clustered with Table 2 Comparison of cDNA sequence between MmBMP2 and MmBMP7 cDNA sequence information MmBMP2 MmBMP7 full-length 2641bp 3161bp ORF 1395bp 1257bp number of amino acids 464aa 418aa 5’UTR 546bp 431bp 3’UTR 703bp 1476bp protein molecular mass 52.906kDa 47.54kDa TGF-β superfamily domain 363Cys-464Arg 317Cys-418His 145 Fig. 1 Neighbor-joining phylogenetic tree of BMP2 and BMP7 between M. meretrix and other selected species using MEGA6.0 software. The abbreviation and the GenBank accession numbers used to construct phylogenetic tree are given in Table 3 Table 3 Species and GenBank accession numbers of BMPs used for multiple alignment and phylogenetic analysis Species Abbreviation GenBank no. Species Abbreviation GenBank no. Meretrix meretrix MmBMP2 ALG64479 Meretrix meretrix MmBMP7 ALG64478 Azumapecten farreri AfBMP2 AGF68558.1 Tegillarca granosa TgBMP7 AFP57673.1 Crassostrea gigas CgBMP2 EKC18750.1 Crassostrea gigas CgBMP7 EKC34211.1 Patella vulgata PvBMP2 AAM33143.1 Pinctada martensii PmBMP7 AGS32053.1 Carassius auratus CaBMP2 BAN17326.1 Danio rerio DrBMP7 AAF17558.1 Danio rerio DrBMP2 AAC25595.1 Ictalurus punctatus IpBMP7 AHH42433.1 Paralichthys olivaceus PoBMP2 BAD16743.1 Xenopus laevis XlBMP7 AAI08478.1 Xenopus tropicalis XtBMP2 AAI70888.1 Rattus norvegicus RnBMP7 EDL85136.1 Bos taurus BtBMP2 AAI42130.1 Equus caballus EcBMP7 ADK93983.1 Mus musculus MuBMP2 EDL28371.1 Mus musculus MuBMP7 EDL06610.1 Pan troglodytes PtBMP2 JAA30020.1 Homo sapiens HsBMP7 AAH08584.1 Homo sapiens HsBMP2 ACV32596.1 146 Fig. 2 (A) Multiple alignment of the MmBMP2 TGF-β domain structures peptide sequence with other mollusc species, (B) Multiple alignment of the MmBMP7 TGF-β domain structures peptide sequence with other mollusc species. Identities are shaded very light red and similarities are shaded lightest red. The * indicate the shared cysteine residues in compared species, the letter box is additional cysteine residue in MmBMP7. (C) Prediction of BMP2 TGF-β domains in mollusc species, (D) Prediction of BMP7 TGF-β domains in mollusc species. Blue boxes indicate transmembrane region, green box indicates coiled coil region, and pink and red boxes mean low complexity region Crassostrea gigas and then with M. meretrix and P. vulgata. In the BMP7 mollusk group, P. martensii first clustered with C. gigas and then with Tegillarca granosa and M. meretrix. MmBMP2 shared the highest identity (53.1%) with A. farreri, and MmBMP7 shared the highest identity (58.7%) with P. martensii. For the mature BMP2 protein of the compared mollusks, the TGF-β domain structure peptide sequences shared seven cysteine residues and the spacing between cysteine residues was well conserved (Fig. 2). The BMP7 protein had seven cysteine residues in M. meretrix but six common cysteines in T. granosa, P. martensii and C. gigas (Fig. 2). Quantitative expression analysis The expression patterns of MmBMP2 and MmBMP7 in the embryos/larvae and adult tissues were analyzed using qRT-PCR. In adult clams, the expression of MmBMP2 in the mantle was extremely significantly higher than that in the other five tissues (p < 0.01) (Fig. 3), suggesting that it might be involved in shell growth. Among the seven developmental stages, MmBMP2 transcript levels were the highest in juvenile clams (p < 0.05), followed by unfertilized mature eggs, blastulae, gastrulae, and eyespot larvae, with the lowest expression in 4-cell embryos and trochophores (Fig. 3). This result may be related to shell formation in the juvenile clams. In adults, MmBMP7 expression was also significantly higher in the mantle than in the other five tissues (p < 0.05) (Fig. 3). The mantle is a shell-forming tissue, and its high expression of BMPs is in agreement with its role in bone formation. Among the developmental stages, MmBMP7 expression was significantly higher in unfertilized mature eggs, blastulae, gastrulae, trochophores, and juvenile clams than in the other two stages (p < 0.05) (Fig. 3), which suggests that it may be involved in morphogenesis and formation of the hard shell. MmBMP2 and MmBMP7 were expressed to different degrees during development, which indicates that BMPs might play a role in ontogeny. Association of SNPs in MmBMP2 and MmBMP7 with growth traits The SNP analysis revealed that eight MmBMP2 SNPs (T134C, A476G, C570T, T585C, A715T, A726G, A732T, and T786C) were present in the JS population, of which three (T585C, A726G, and A732T) were associated with growth traits. Five SNPs (G390A, C570T, A726G, A732T, and T786C) in exons of MmBMP2 were identified in the WL2 population, including two (A726G and A732T) that were potentially associated with growth traits. JS individuals with genotype TC at position T585C grew faster than those with the TT genotype in terms of shell length, shell width, and shell height (p < 0.05). 147 Fig. 3 Analysis of expression differences of MmBMP2 and MmBMP7. (A) Analysis of expression difference of MmBMP2 in six tissues, (B) Analysis of expression difference of MmBMP7 in six tissues, (C) Analysis of expression difference of MmBMP2 in seven developmental stages, (D) Analysis of expression difference of MmBMP7 in seven developmental stages. The same letters indicate no difference in the level of expression, whereas different letters indicate significant differences in expression levels among various developmental stages (p < 0.05). **p < 0.01, *p < 0.05 Clams with genotype AG at A726G exhibited faster growth in all growth traits except shell weight compared with AA individuals (p < 0.05). The SNP A732T AT genotype was also associated with significantly higher growth traits than the AA genotype (p < 0.05) (Table 4). WL2 individuals with genotype AG at position A726G had significantly better growth traits compared with genotype AA (p < 0.05). Multiple comparison analysis showed that individuals with the A732T AT genotype grew faster than those with the AA genotype in all growth traits measured (p < 0.05) (Table 4). Comparison of the loci between the JS and WL2 groups revealed that the populations shared two SNPs at A726G and A732T in the MmBMP2 coding region. Association analysis showed that individuals with the AG genotype at locus 726 had significantly higher growth traits than those with the AA genotype in both groups (p < 0.05). Additionally, individuals with the AT genotype at locus 732 had significantly higher growth traits than those with the AA genotype in both groups (p < 0.05). SNP A732T was significantly associated with better outcomes for all growth traits tested in both the JS and WL2 populations (p < 0.05). The SNP analysis revealed the presence of three MmBMP7 SNPs (C32T, C126G, and G144A) in the JS population, including two (C32T and G144A) associated with growth traits. JS individuals with genotype CT at position C32T grew faster than those with the CC genotype in terms of shell height, shell weight, soft weight, and shell weight (p < 0.05). Clams with genotype GA at G144A exhibited faster growth in terms of shell width, 148 Table 4 Correlation analysis of SNPs in the MmBMP2 with growth traits in two groups Note: bold parts are the SNPs associated with growth traits, p < 0.05; underline parts are the SNPs potentially associated with growth traits, p < 0.01 total weight, soft weight, and shell weight compared with GG individuals (p < 0.05) (Table 5). Three SNPs (C32T, C126G, and G306A) were found in exons of MmBMP7 from the WL2 population, but none was associated with growth traits. Discussion In vertebrates, every mature BMP shares seven conserved cysteines (Bragdon et al., 2011), six of which build a cystine knot by forming three intrachain disulfide bonds to resist heat, denaturants, and extreme pH levels. The remaining cysteine residue forms an interchain disulfide bond that links two monomers into active hetero- or homodimers (Kingsley, 1994). This unique cysteine is absent in some mollusks (T. granosa, C. gigas, and P. martensii) (Yan et al., 2014). This absence raises the possibility that BMP7 in some mollusks might not form dimers or might form noncovalently linked dimers between the subunits, which could be dynamic and subject to regulation. However, both BMP2 and BMP7 in M. meretrix contained all seven cysteine residues. Gene expression analysis showed that expression of MmBMP2 and MmBMP7 was the highest in the mantle, which suggests that they might be involved in shell growth in M. meretrix. In P. martensii, expression of BMP2 and BMP7 was also the highest in the mantle (Miyashita et al., 2008; Yan et al., 2014). Previous studies showed that BMPs of mollusks may play an important role in formation of the shell boundary. Using northern blot and in situ hybridization assays, Miyashita et al. (2008) found that expression of BMP2 occurred predominantly in the mantle, suggesting that BMP2 played a key role in nacre formation. Zhou et al. (2010) reported that P. martensii BMP2 was involved in repair after shell incision. Furthermore, the expression of BMP7b increased significantly 24 h after shell incision in P. martensii. In the current study, both MmBMP7 and MmBMP2 were expressed in all six tissues tested, which suggests that they may be involved in a wide range of biological processes. For example, mGDF, which has high similarity to BMP2, in C. gigas, along with other growth factors, is involved in differentiation of the digestive gland from stem cells into functional cells (Lelong et al., 2010). Expression of MmBMP2 and MmBMP7 was detected at every embryonic stage of M. meretrix during development, suggesting that they may be involved in differentiation of early tissues and organs. Bivalves have a common process of shell formation that occurs at the very early stage of larval development. During the gastrula stage, the shell formative region sinks, evaginates and enlarges until Group Locus Genotype Sample Percent (%) Shell length (cm) Shell width (cm) Shell height (cm) Total weight (g) Soft weight (g) Shell weight (g) Jiangsu 585 T>C TT 121 88.97 4.11±0.61a 2.05±0.32a 3.32±0.48a 17.28±6.68 8.38±3.37 8.89±3.42 TC 15 11.03 4.52±0.57b 2.27±0.33b 3.65±0.47b 20.51±5.86 9.83±2.70 10.68±3.39 726 A>G AA 63 46.32 3.97±0.62a 1.97±0.33a 3.20±0.48a 15.50±6.32a 7.49±3.25a 8.01±3.16a AG 50 36.76 4.26±0.60b 2.14±0.32b 3.42±0.48b 18.35±6.21b 8.92±3.05b 9.43±3.29b GG 23 16.91 4.47±0.45b 2.22±0.24b 3.64±0.40b 21.92±6.34c 10.61±3.08 11.31±3.45c 732 A>T AA 70 51.47 4.01±0.63a 2.00±0.34a 3.24±0.49a 15.91±6.39a 7.68±3.25a 8.23±3.22a TA 50 36.76 4.29±0.59b 2.15±0.32b 3.44±0.47b 18.95±6.47b 9.27±3.19b 9.68±3.43b TT 16 11.76 4.40±0.47b 2.17±0.24 3.60±0.42b 21.04±6.52b 10.05±3.20 10.99±3.50b Wanli 2 726 A>G AA 41 32.28 35.95±2.56a 17.28±1.34a 30.04±2.23a 11.35±2.70a 5.57±1.31a 5.78±1.44 AG 59 46.46 37.36±3.14b 17.96±1.68b 31.07±2.51b 12.71±3.43b 6.30±1.60b 6.42±1.88 GG 27 21.26 36.18±3.27 17.34±1.59 30.29±2.77 11.64±3.60 5.70±1.67 5.94±1.96 732 A>T AA 51 40.16 36.25±2.63a 17.40±1.38a 30.28±2.32 11.56±2.83a 5.70±1.33a 5.86±1.55a AT 52 40.94 37.51±3.22b 18.09±1.70b 31.17±2.52a 12.95±3.51b 6.38±1.64b 6.57±1.91b TT 24 18.90 35.70±3.11a 17.05±1.46a 29.94±2.69b 11.15±3.36a 5.50±1.62a 5.65±1.77a 149 Table 5 Correlation analysis of SNPs in the MmBMP7 with growth traits in two groups Note: bold parts are the SNPs associated with growth traits, p < 0.05; underline parts are the SNPs potentially associated with growth traits, p < 0.01 the whole embryo is encapsulated (Kin et al., 2009). Afterward, the conchiolin is secreted and calcified, and simultaneously the shell forms a hinged part in the back and divides into two pieces until bivalve stage I (prodissoconch I) is reached. During formation of the velum, the phase II embryo shell (prodissoconch II) continues to form (Kakoi et al., 2008). Miyazaki et al. (2010) reported that the shell of P. martensii begins to form prodissoconch I at the trochophore stage and then forms prodissoconch II at the D-shaped larvae stage. In the current study, MmBMP2 and MmBMP7 were highly expressed at the trochophore stage, which suggests that they might play an important role in the formation of the protegulum. Moreover, expression of MmBMP2 was the highest in the juvenile clam stage, which indicates that organs and shells grew fast and developed rapidly at this stage. Molecular marker-assisted breeding can be used to select new varieties with certain traits at the molecular level and thereby greatly improve breeding efficiency. SNPs, as the third generation of molecular markers, are a powerful tool for assessing genetic diversity and conducting association analysis. In marine mollusks, association analysis between SNPs of candidate genes and economic traits such as growth and stress tolerance has been conducted. For example, Guo et al. (2012) identified a growth-related SNP locus in the 3'UTR region of TGF-β receptor type 1 in Chlamys farreri. In the exon region of myostatin in Argopecten irradians, two SNP loci associated with growth were identified (Guo et al., 2011). Bao et al. (2013) characterized three immune-related SNP loci in the hemoglobin of T. granosa. In the current study, we identified three SNPs potentially associated with growth traits in MmBMP2 of the JS population and two in the WL2 population, and two of them were common to the two populations. For MmBMP7, we found two SNPs potentially associated with growth traits in the JS population and none in the WL2 population. The individuals with the AG genotype of SNP A726G grew significantly faster than those with the AA genotype, and individuals with the AT genotype of SNP A732T grew significantly faster than those with the AA genotype in both populations. Furthermore, these two SNPs were located within the coding region of MmBMP2 without amino acid changes (A726G encoding threonine, A732T encoding valine). Many studies have shown that synonymous mutations regulate gene transcription and translation. Greenwood and Kelsoe (2003) found that the synonymous mutation of many genes in promoters and introns could affect transcriptional efficiency. Kimchi-Sarfaty et al. (2007) reported that the synonymous mutation of the MDR1 gene could change the spatial structure of the protein. Capon et al. (2004) demonstrated that synonymous mutations of the corneodesmosin gene could improve the stability of mRNA molecules. Therefore, we hypothesize that clams with the A726G AG genotype or the A732T AT genotype are favorable for breeding. Furthermore, the MmBMP2 SNP could influence growth performance and may be a suitable marker for marker-assisted selection in these M. meretrix populations. Acknowledgments This work was financially supported by National Natural Science Foundation of China (31772846), Modern Agro-industry Technology Research System (CARS-49), National Infrastructure of Fishery Germplasm Resources Programme (2018DKA30470) and Zhejiang Provincial First-Class Discipline of Bioengineering -A （ZS2019001）. References Amthor H, Christ B, Weil M, Patel K. The importance of timing differentiation during limb muscle development. Curr. Biol. 8(11): 642-652, 1998. Bao Y, Li P, Dong Y, Xiang R, Gu L, Yao H, et al. Polymorphism of the multiple hemoglobins in blood clam Tegillarca granosa and its association with disease resistance to Vibrio parahaemolyticus. Fish Shellfish Immun. 34(5): 1320-1324, 2013. Group Locus Genotype Sample Percent (%) Shell length (cm) Shell width (cm) Shell height (cm) Total weight (g) Soft weight (g) Shell weight (g) Jiangsu 32 C>T CC 90 66.67 4.12±0.60 2.06±0.32 3.33±0.47a 17.03±6.36a 8.19±3.14a 8.83±3.34a CT 27 20.00 4.42±0.53 2.21±0.28 3.58±0.46b 20.98±6.08b 10.43±3.0b 10.55±3.20b TT 18 13.33 4.08±0.69 2.05±0.37 3.37±0.49 16.67±7.34a 7.97±3.69a 8.86±3.68a 144 G>A GG 78 57.78 4.08±0.60a 2.03±0.31a 3.29±0.47a 16.56±6.25a 8.03±3.17a 8.52±3.21a GA 38 28.15 4.33±0.59b 2.15±0.31b 3.49±0.48b 19.91±6.91b 9.63±3.41b 10.28±3.65b AA 19 14.07 4.31±0.61 2.20±0.37b 3.49±0.52 18.48±6.42 8.93±3.23 8.82±3.34 150 Bragdon B, Moseychuk O, Saldanha S, King D, Julian J, Nohe A. Bone morphogenetic proteins: a critical review. Cell. Signal. 23(4): 609-620, 2011. Capon F, Allen MH, Ameen M, Burden AD, Tillman D, Barker JN, et al. A synonymous SNP of the corneodesmosin gene leads to increased mRNA stability and demonstrates association with psoriasis across diverse ethnic groups. Hum. Mol. Genet. 13(20): 2361-2368, 2004. Chen D, Zhao M, Mundy GR. Bone morphogenetic proteins. Growth Fact. 22(4): 233-241, 2003. Cook SD, Rueger DC. Osteogenic protein-1: biology and applications. Clin. Orthop. Relat. R. 324: 29, 1996. Derynck R, Akhurst RJ. Differentiation plasticity regulated by TGF-β family proteins in development and disease. Nat. Cell Biol. 9(9): 1000-1004, 2007. Dudley AT, Lyons KM, Robertson EJ. A requirement for bone morphogenetic protein-7 during development of the mammalian kidney and eye. Gene. Dev. 9(22): 2795-2807, 1995. Dudley AT, Robertson EJ. Overlapping expression domains of bone morphogenetic protein family members potentially account for limited tissue defects in BMP7 deficient embryos. Dev. Dyn. 208(3): 349-362, 1997. Freitas TC, Jung E, Pearce EJ. A bone morphogenetic protein homologue in the parasitic flatworm, Schistosoma mansoni. Int. J. Parasitol. 39(3): 281-287, 2009. Fritsch C, Lanfear R, Ray RP. Rapid evolution of a novel signalling mechanism by concerted duplication and divergence of a BMP ligand and its extracellular modulators. Dev. Genes Evol. 220(9-10): 235-250, 2010. Gao XY, Dong YH, Shi SJ, Yao HH, Ruan WB, Zhao JX, et al. Cloning, spatiotemporal expression and SNPs identification of GRB2 gene in hard clam Meretrix meretrix. J. Fish. China. 39 (9): 1324-1332, 2015. Greenwood TA, Kelsoe JR. Promoter and intronic variants affect the transcriptional regulation of the human dopamine transporter gene. Genomics. 82(5): 511-520, 2003. Guo H, Bao Z, Li J, Lian S, Wang S, He Y, et al. Molecular characterization of TGF-β type I receptor gene (Tgfbr1) in Chlamys farreri, and the association of allelic variants with growth traits. Plos One. 7(11): e51005, 2012. Guo L, Li L, Zhang S, Guo X, Zhang G. Novel polymorphisms in the myostatin, gene and their association with growth traits in a variety of bay scallop, Argopecten irradians. Anim. Genet. 42(3): 339-340, 2011. Herpin A, Lelong C, Becker T, Rosa F, Favrel P, Cunningham C. Structural and functional evidence for a singular repertoire of BMP receptor signal transducing proteins in the lophotrochozoan Crassostrea gigas suggests a shared ancestral BMP/activin pathway. Febs J. 272(13): 3424-3440, 2010. Hospital F, Moreau L, Lacoudre F, Charcosset A, Gallais A. More on the efficiency of marker-assisted selection. Theor. Appl. Genet. 95(8): 1181-1189, 1997. Hwang SP, Partin JS, Lennarz WJ. Characterization of a homolog of human bone morphogenetic protein 1 in the embryo of the sea urchin, Strongylocentrotus purpuratus. Development. 120(3): 559-568, 1994. Iijima M, Takeuchi T, Sarashina I, Endo K. Expression patterns of engrailed, and dpp, in the gastropod Lymnaea stagnalis. Dev. Genes Evol. 218(5): 237-251, 2008. Kakoi S, Kin K, Miyazaki K, Wada H. Early development of the Japanese spiny oyster (Saccostrea kegaki): characterization of some genetic markers. Zool. Sci. 25(5): 455-464, 2008. Kimchi-Sarfaty C, Oh JM, Kim IW, Sauna ZE, Calcagno AM, Ambudkar SV, et al. A "silent" polymorphism in the MDR1 gene changes substrate specificity. Science. 315(5811): 525-528, 2007. Kin K, Kakoi S, Wada H. A novel role for dpp, in the shaping of bivalve shells revealed in a conserved molluscan developmental program. Dev. Biol. 329(1): 152-166, 2009. Kingsley DM. The TGF-beta superfamily: new members, new receptors, and new genetic tests of function in different organisms. Gene. Dev. 8(2): 133-146, 1994. Lande R, Thompson R. Efficiency of marker-assisted selection in the improvement of quantitative traits. Genetics. 124(3): 743-756, 1990. Lelong C, Mathieu M, Favrel P. Structure and expression of mGDF, a new member of the transforming growth factor-beta superfamily in the bivalve mollusc Crassostrea gigas. Eur. J. Biochem. 267(13): 3986-3993, 2010. Li Q, Agno JE, Edson MA, Nagaraja AK, Nagashima T, Matzuk MM. Transforming growth factor β receptor type 1 is essential for female reproductive tract integrity and function. Plos Genet. 7(10): e1002320, 2011. Liu B, Dong B, Tang B, Zhang T, Xiang J. Effect of stocking density on growth, settlement and survival of clam larvae, Meretrix meretrix. Aquaculture. 258(1): 344-349, 2006. Luo G, Hofmann C, Bronckers AL, Sohocki M, Bradley A, Karsenty G. BMP-7 is an inducer of nephrogenesis, and is also required for eye development and skeletal patterning. Gene. Dev. 9(22): 2808-2820, 1995. Lyons KM, Hogan BL, Robertson EJ. Colocalization of BMP7 and BMP2 RNAs suggests that these factors cooperatively mediate tissue interactions during murine development. Mech. Dev. 50(1): 71-83, 1995. Massague J. TGF-β signal transduction. Adv. Cem. Based Mater. 2(95): 30-38, 1998. Miyashita T, Hanashita T, Toriyama M, Takagi R, Akashika T, Higashikubo N. Gene cloning and biochemical characterization of the BMP-2 of Pinctada fucata. Biosci. Biotech. Bioch. 72(1): 37-47, 2008. Miyazaki Y, Nishida T, Aoki H, Samata T. Expression of genes responsible for biomineralization of Pinctada fucata during development. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 155(3): 241-248, 2010. 151 Moreau L, Charcosset A, Hospital F, Gallais A. Marker-assisted selection efficiency in populations of finite size. Genetics. 148(3): 1353-1365, 1998. Musgrave DS, Pruchnic R, Wright V, Bosch P, Ghivizzani SC, Robbins PD, et al. The effect of bone morphogenetic protein-2 expression on the early fate of skeletal muscle-derived cells. Bone. 28(5): 499-506, 2001. Nederbragt AJ, van Loon AE, Dictus WJ. Expression of Patella vulgata orthologs of engrailed and dpp-BMP2/4 in adjacent domains during molluscan shell development suggests a conserved compartment boundary mechanism. Dev. Biol. 246(2): 341-355, 2002. Reinhardt B, Broun M, Blitz IL, Bode HR. HyBMP5-8b, a BMP5-8 orthologue, acts during axial patterning and tentacle formation in hydra. Dev. Biol. 267(1): 43-59, 2004. Salazar VS, Gamer LW, Rosen V. BMP signalling in skeletal development, disease and repair. Nat. Rev. Endocrinol. 12(4): 203-221, 2016. Schmal H, Mehlhorn AT, Pilz IH, Doviakue D, Kirchhoff C, Sdkamp NP, et al. Immunohistological localization of BMP-2, BMP-7, and their receptors in knee joints with focal cartilage lesions. The Scientific World J. 2012-1-3(6): 467892, 2012. Shih LJ, Chen C, Chen CP, Hwang SP. Identification and characterization of bone morphogenetic protein 2/4 gene from the starfish Archaster typicus. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 131(2): 143-151, 2002. Tang B, Liu B, Wang G, Zhang T, Xiang J. Effects of various algal diets and starvation on larval growth and survival of Meretrix meretrix. Aquaculture. 254(1): 526-533, 2006. Wang C, Wang H, Li Y, Liu B. Identification of a fructose-1,6-bisphosphate aldolase gene and association of the single nucleotide polymorphisms with growth traits in the clam Meretrix meretrix. Mol. Biol. Rep. 39(4): 5017-5024, 2012. Wang C, You Y, Wang H, Liu B. Genetic diversity of the sulfotransferase-like gene and one nonsynonymous SNP associated with growth traits of clam, Meretrix meretrix. Mol. Biol. Rep. 39(2): 1323-1331, 2012. Wang H, Chai X, Liu B. Estimation of genetic parameters for growth traits in cultured clam Meretrix meretrix (Bivalvia: Veneridae) using the Bayesian method based on Gibbs sampling. Aquac. Res. 42(2): 240-247, 2011. Werner E.G. Müller, Binder M, Lintig JV, Guo YW, Wang X, Kaandorp JA, et al. Interaction of the retinoic acid signaling pathway with spicule formation in the marine sponge Suberites domuncula through activation of bone morphogenetic protein-1. BBA-Gen. Subjects. 1810(12): 1178-1194, 2011. Wharton KA, Thomsen GH, Gelbart WM. Drosophila 60A gene, another transforming growth factor β family member, is closely related to human bone morphogenetic proteins. P. Natl. Acad. Sci. USA. 88(20): 9214-9218, 1991. Xiao YT, Xiang LX, Shao JZ. Bone morphogenetic protein. Biochem. Biophys. Res. Commun. 362(3): 550-553, 2007. Xie C, Xu SZ. Efficiency of multistage marker-assisted selection in the improvement of multiple quantitative traits. Heredity. 80(4): 489-498, 2010. Yan F, Luo S, Jiao Y, Deng Y, Du X, Huang R, et al. Molecular characterization of the BMP7 gene and its potential role in shell formation in Pinctada martensii. Int. J. Mol. Sci. 15(11): 21215-21228, 2014. Yao X, Zhang J, Sun J, Liu B. Recombinant expression, characterization and expressional analysis of clam Meretrix meretrix, cathepsin B, an enzyme involved in nutrient digestion. Mol. Biol. Rep. 38(3): 1861-1868, 2011. Zhou YJ, He ZX, Li CZ, Xiang L, Zhu FJ, Zhang YG, et al. Correlations among mRNA expression levels of engrailed, BMP2 and Smad3 in mantle cells of pearl oyster Pinctada fucata. Prog. Biochem. Biophys. 37(7): 737-746, 2010. Zoccola D, Moya A, Béranger GE, Tambutté E, Allemand D, Carle GF, et al. Specific expression of BMP2/4 ortholog in biomineralizing tissues of corals and action on mouse BMP receptor. Mar.Biotechnol. 11(2): 260-269, 2009.