Mantle transcriptome sequencing of Mytilus spp. and identification of putative biomineralization genes

RNA extraction and pyrosequencing

Five populations of mussels (Mytilus spp.) representing four taxa were collected from several localities: M. edulis (North Sea, EduN), M. trossulus (Vancouver, TroV), M. galloprovincialis (Mediterranean Sea, GalM), M. galloprovincialis (Tasmania, GalT), M. chilensis (Chile, Chil) (Table 1). All taxa had been identified by single nucleotide polymorphism (SNP) genotyping (Zbawicka et al., 2012; Wenne et al., 2016; Larraín et al., 2017). Specimens chosen for RNA extraction were adults of indeterminate sex, spent. Tissue sampling was standardized across individuals: extraction of total RNA from the central and outer mantle tissue from right and left mussel side was performed using a GenElute Mammalian Total RNA Miniprep Kit (Sigma-Aldrich, St. Louis, MO, USA) according to the manufacturer’s protocol, with minor modifications. Homogenization was carried out in a custom buffer composed of 100 mM Tris, 1.4M NaCl, 20 mM EDTA, 2% CTAB, proteinase K (20 mg/ml) and 0.03 mM 2-mercaptoethanol (Malachowicz, Wenne & Burzynski, 2017). Samples were stored at −20 °C. RNA samples from two to four individuals were pooled for each group in equal proportions to make a total of five libraries. The quantity and quality of total RNA was determined at 260 nm on microplate using an Epoch Microplate Spectrophotometer (BioTek Instruments, Incorporated, Winooski, VT, USA) and gel electrophoresis. RNA integrity was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). cDNA libraries were prepared using a cDNA Rapid Library Preparation Kit (Roche Applied Science, Basel, Switzerland) and pyrosequencing was performed by Macrogen Incorporation (Seoul, Korea), using Roche GS-FLX sequencing system (Roche, Basel, Switzerland) according to manufacturer’s instructions.

De novo assembly, annotation, SSR and SNP discovery

All reads from the mantle tissue of five Mytilus groups were processed separately, through the same pipeline. Raw reads were pre-processed by removing adapters and low quality sequences (quality score = 0.05, discard reads <50 bp) with CLC Genomics Workbench software (v.7.5.5, CLC Bio, Qiagen, Aarhus, Denmark). Quality reads were assembled into contigs using de-Bruijn graphs in CLC Genomics Workbench with default parameters (Malachowicz, Kijewska & Wenne, 2015; Malachowicz, Wenne & Burzynski, 2017). Contigs were first searched against proteins from the NCBI non-redundant (NR) database using a BLASTX tool implemented in BLAST+ (v.2.2.29) (Altschul et al., 1990), with an E-value threshold of 10⁻⁵ and other default parameters. A list of the top BLAST hits was achieved by using bit score (12th column) and E-value values (11th column) with the code: sort -k1,1 -k12,12 nr -k11,11 g -k3,3 gr blastout.txt sort -u -k1,1–merge > besthit. For functional annotation, gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were assigned to the transcripts using Blast2GO software (Conesa et al., 2005) with E-value Hit-Filter = 1.0 E⁻¹⁰ (other default parameters) and the single directional best-hit method in the KEGG Automatic Annotation Server (Moriya et al., 2007), respectively.

The software, MSDB (Du et al., 2013) was used to discover potential microsatellite markers. Di- to hexa- nucleotide motifs, with minimum thresholds of 6-5-4-4-4 repeats, respectively and 50 bp of flanking sequence length were identified. SNP analyses were carried out using CLC Genomics Workbench. For each transcriptome, clean reads were mapped back to corresponding contigs with an overlap criterion of 70% and a similarity of 90%. Candidate variants were called with the following minimum value parameters: neighborhood quality = 15; quality of central bases = 20; coverage = 20; count = 8; variant frequency = 35.0.

Identification of potential biomineralization and melanogenesis-related contigs

Candidate contigs related to biomineralization and melanogenesis were identified by searching top BLASTX hits with key words from a reference list of previously reported candidate biomineralization genes (Table S1). Melanogenesis-related transcripts were identified based on KEGG pathways analysis.

Several biomineralization proteins were selected for additional analysis: carbonic anhydrase (CA), chitinase (CHIT), tyrosinase (TYR), dermatopontin (DPT), arginine kinase (AK) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH). Proteins of selected sequences were predicted using OrfPredictor (Min et al., 2005). Functional domains and motifs were identified using Simple Modular Architecture Research Tool (SMART) software, allowing identification and annotation of domain architectures (http://smart.embl-heidelberg.de) (Letunic, Doerks & Bork, 2015). To further scrutinize and compare selected contigs, representative reference bivalve databases were created. Protein sequences described as CA, CHIT, DPT, AK and GAPDH, available at the GenBank database (www.ncbi.nlm.nih.gov) were downloaded (217, 155, 8, 43, 24 sequences, respectively). Only sequences with annotated domains, were selected (199, 123, 8, 43, 24 sequences, respectively). Following this, mussel contigs from the present study were searched against the created databases with an E-value threshold of 10⁻¹⁵. For phylogenetic reconstruction of Mytilus spp. TYR proteins, contigs from this study, amino acid sequences from GenBank (Acc. No: AKS48166, ALA16023, OPL3388), and sequences from previous studies (Aguilera, McDougall & Degnan, 2014; Moreira et al., 2015; Hüning et al., 2016; Sleight et al., 2016) were used. Homologous sequences were aligned using multiple alignment program MUSCLE (Edgar, 2004) with default parameters, and trimmed using TrimAl (Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009) with an 80% threshold. Best substitution models were inferred using ModelFinder (Kalyaanamoorthy et al., 2017), as follow: CA—the Le Gascuel (LG) model (Le & Gascuel, 2008), gamma distribution, four discrete rate categories (+G4); GH18—the Whelan and Goldman (WAG) model (Whelan & Goldman, 2001) + G4 with invariant sites (+I); AK and GAPDH—WAG+ G4; DPT—the Block substitution matrices (Blosum62) model (Henikoff & Henikoff, 1992) + G4, with frequencies (+F); TYR—Blosum62+F+I+G4 model. Three phylogenetic models were run for each alignment. Neighbor-joining (NJ) trees were performed using MEGA X (Kumar et al., 2018) and the Jones–Taylor–Thornton (JTT) substitution model (Jones, Taylor & Thornton, 1992). Maximum likelihood reconstructions were performed using IQ-TREE (Nguyen et al., 2015) with 1,000 bootstrap replicates. Bayesian inferences were constructed using MrBayes v.3.2 (Ronquist et al., 2012) with parameters: ngen = 1,500,000; samplefreq = 100; printfreq = 100; nchain = 4; starttree = random. Trees were visualized and edited using MEGA X and manually combined to produce a consensus tree for each protein alignment.

Comparative transcriptomic analysis

Homologous transcripts present in all five groups were identified by using bi-directional tBLASTx analyses conducted between each transcriptome with an E-value cut-off of 10⁻¹⁵. A list of top BLAST hits was achieved by using bit score (12th column) and E-value values (11th column) with the code: sort -k1, 1 -k12, 12 nr -k11, 11 g -k3,3 gr blastout.txt | sort -u -k1,1–merge>besthit. Further, results from this study were compared with available Mytilus spp. mantle sequences. The accession numbers are: for M. edulis: PRJEB4516 (Freer et al., 2014), PRJNA30561 (Yarra et al., 2016); for M. galloprovincialis: PRJNA230138 (Moreira et al., 2015), PRJNA295512 (Bjärnmark et al., 2016). For this purpose, a de novo transcriptome assembly was performed using the raw reads deposited in the NCBI Sequence Read Archive (SRA). The Illumina HiSeq raw reads were filtered using Trimmomatic v.0.32 (Bolger, Lohse & Usadel, 2014) and retained reads were assembled using Trinity program with default parameters (Grabherr et al., 2011).

To find enriched GO terms between groups, Fisher’s exact test was performed with Blast2GO software, applying a one-tailed test, removing double IDs, with a false discovery rate (FDR) cut-off of 0.05. As a test-set, list of contigs from each group was used. As a background (reference-set), a fusion of the five transcriptomes was used.

To identify putative orthologs among the genus of Mytilus, transcripts from this study and from M. coruscus (Xu et al., 2016), M. californianus (Romiguier et al., 2014) were analyzed using OrthoMCL software with default settings (Li, Stoeckert & Roos, 2003). The putative orthologous clusters were aligned with a MUSCLE tool and then concatenated in CLC Genomics Workbench. A splits network was created with SplitsTree 4.13.1 (Huson & Bryant, 2006), using the neighbor-net algorithm under pdistance model. A phylogenetic tree was inferred using a NJ model in MEGA X, with the JTT + G4 substitution model and 1,000 bootstrap replicates.

Assembly, annotation, SSR and SNP detection

A total of 753,135 raw reads with an average read length of 342 bp were generated from sequencing of the Mytilus spp. mantle transcriptomes. After filtering, 743,773 (98.75%) high-quality reads were used for de novo assembly. In total, 20,982 contigs with an average length of 529 bp were obtained (Table S2). Mussel raw sequences from this study were deposited in NCBI SRA under accession number PRJNA419475.

BLAST analysis against NCBI NR database showed that between 43.19% (TroV) and 52.18% (GalM) of transcripts had a significant hit (Tables S2 and S3). As expected, most of the annotated transcripts were homologous to proteins from the Mollusca phylum: Pacific oyster (C. gigas, av. 50%), Mytilus spp. (av. 10%) and owl limpet (L. gigantea, av. 5%) (Fig. S1). KEGG analysis showed that a total of 3,679 (17.53%) transcripts were assigned to 271 biochemical pathways. The number of sequences assigned to each pathway was different across groups (Table S2). However, the most represented pathways in all groups included focal adhesion, phagosome, ribosome and apoptosis (Table S4). The signal transduction, endocrine system, immune system and cellular community were the most represented KO (KEGG orthology) subcategories in all groups (Fig. S2). In total, 35,383 GO terms were assigned to 7,824 (37.28%) contigs (Table S2). The main GO category distribution was similar in all Mytilus spp. groups, the dominant was the biological process category (BP), followed by cellular component and molecular function (MF) (Fig. S3). In this study, 483 potential simple sequence repeats (SSRs) were identified in 645 contigs (3.07% of the total number of contigs) (Table S5A). Among the identified putative markers, the most highly represented in the EduN, TroV, GalM and GalT groups were di-nucleotide repeats, and tri-nucleotide in Chil. AT was the dominating motif in all groups (Table S5B). In addition, a total of 1,497 putative SNPs were detected in all groups. The A/G variation type was dominant in EduN and TroV groups, whereas in other groups the C/T was dominant (Fig. S4).

Contigs potentially involved in biomineralization and melanogenesis

BLASTx searching revealed 1,292 contigs with significant sequence similarity to candidate biomineralization proteins (Table S6). SMPs such as, nacrein, calmodulin, collagens, perlucin, perlustrin, PIF were found among identified contigs (see Table S6). Transcripts involved in ion and acid-base regulation, protein homeostasis and anti-oxidation were also found including V-type proton ATPases, sodium/potassium-transporting ATPases, heat shock 70 kDa protein (HSP70) and glutathione peroxidase. Other annotations included C1q-domain-containing protein, cAMP-response element-binding protein and cold-shock protein (CSP) (Table S6). Of annotated biomineralization contigs, 68.65% had domains, including chitin binding, EF hand, sushi, Von Willebrand factor A (Table S6).

Several biomineralization proteins were selected for more comprehensive analyses. In M. galloprovincialis, M. edulis and M. trossulus groups, five contigs (GalT_1066, GalT_1436, EduN_2761, GalM_2384, TroV_1436) were identified as DPT with DERM domain. A phylogenetic analysis of DPT proteins in Bivalvia generated a separate, well-supported clade of Mytiloida members, suggesting that these contigs are Mytilus spp. orthologues of DPT (Fig. 1A). Contigs with 81–95% sequence identity to AK were identified in all groups. The phylogenetic analysis showed two separate clades of Pterimorphia and Heteroconchia subclasses (Fig. 1B). Contigs Chil_958, EduN_492, GalM_470 and TroV_123 encoded proteins with similarity to GAPDH and possessed the Gp_dh_C domain. Alignment with other known GAPDHs from Ostreoida, Pectinoida and Pterioida members revealed many conserved regions among these species (Fig. S5). However, a phylogenetic analysis showed two well-supported separate clades of Mytiloida and Ostreoida with Pectinoida orders (Fig. 1C). TYR amino acid sequences derived were aligned with available molluscan TYR sequences. A phylogenetic analysis revealed that mussel TYR, in this study, created two separate clades: clade A and clade B (Fig. 2), described earlier as ancestral and bivalve-specific (Aguilera, McDougall & Degnan, 2014). It is possible that there were two small expansions within Mytilus, which is a common feature in bivalves. In M. galloprovincialis there were TyrA2, TyrA3 orthologues shared with the genera Crassostrea and Pinctada, indicating gene duplication before the divergence of lineages A and B. Searches in mussel transcriptomes revealed five contigs (in EduN, TroV and GalT groups) that pertain to the CA gene family and 10 contigs (in all transcriptomes) that showed similarity to the glycoside hydrolase family 18 (GH18) CHIT members (Table 2). A SMART analysis revealed CA domains in annotated CA sequences, and GH18_CHIT-like superfamily domain in GH18 sequences. The chitin binding peritrophin-A domain (CBM_14 superfamily) was found in transcripts recognized as acidic mammalian CHIT (GalT_650, GalT_649 and TroV_158). The phylogenetic analysis of CA in the Bivalvia class revealed three major clusters: α-carbonic anhydrases related proteins (α-CARPs), cytosolic/extracellular α-carbonic anhydrases (α-CAs) and membrane-bound α-CAs. There were two visible sub-clusters: CA2 and CA14 (Fig. 3A). Two major clades, A and B were resolved in the GH18 CHIT family phylogenetic tree. The A cluster consisted of CHIT domain-containing proteins (CHIDs, with GH18_SI-CLP domain), and Di-N-acetylchitobiases. In contrast to A, cluster B contained genes involved in chitin degradation (CHIT and CHIT-like proteins), creating three sub-clades (B1, B2, B3) (Fig. 3B). Presented phylogenetic reconstructions were supported by the high values of the nodes.

A total of 42, 99, 19 and 62 contigs from the mussel mantle transcriptomes were assigned to 13, 40, 8 and 15 KO terms in melanogenesis, endocytosis, TYR and calcium signaling KEGG pathways, respectively (Table S7).

Comparative transcriptomics

Bi-directional Best Hit (BBH) method was used to identify 6,130 homologous pairs between transcriptomes of which 552 were shared between all five groups (Fig. S6). The highest number of homologs were found between pairs: EduN_TroV, EduN_GalT and EduN_GalM (1,029, 946, 818, respectively). Putative homologous pairs were ordered according to sequence similarity of BLAST alignments (Fig. S6). In general most homologous pairs exhibited a sequence similarity of 90–100%. The highest similarity per no. of homologs was found between EduN_GalM (85.45%), followed by EduN_Chil (84.89%), GalM_Chil (84.55%), GalM_GalT (84.12%). The lowest similarity to other groups was shown by TroV. Transcriptomes obtained from this study were compared to available mantle sequences for M. edulis and M. galloprovincialis using tBLASTx. In total, 18,456 (87.96%) contigs from this study were homologous with one of the four created mantle transcriptomes, but 12.04% might represent new sequence information. Additionally, between 47.8% and 81.5% of homologous contigs had high similarity (>90%; Fig. S7).

For more accurate analysis, Fisher’s exact test was applied to discover significantly over-represented GO terms in each group. The most interesting GO enrichments associated with BP and MF were identified only in EduN and GalM groups. The results were reduced to the most specific terms (FDR = 0.01) (Figs. S8 and S9). A total of 25 GO terms were enriched in the M. edulis from the North Sea including cellular response to chemical stimulus (GO:0070887), positive regulation of cellular process (GO:0048522) in the BP category, ATP binding (GO:0005524), calcium ion binding (GO:0005509) in the MF category (Fig. S8). A total of 46 terms were enriched in the M. galloprovincialis from the Mediterranean Sea such as organonitrogen compound biosynthetic process (GO:1901566), regulation of biological quality (GO:0065008) in the BP category, metal ion binding (GO:0046872), purine ribonucleotide binding (GO:0032555) in the MF category (Fig. S9).

Both the neighbor-net splits network (Fig. 4A) and NJ phylogenetic tree (Fig. 4B) showed that M. chilensis was the most similar to the M. edulis taxon, followed by M. galloprovincialis and M. trossulus. A separate clade was identified containing M. californianus Conrad, 1837 and M. coruscus Gould, 1861 (presently accepted as M. unguiculatus Valenciennes, 1858) mussels (Fig. 4).