Unbiased high-throughput characterization of mussel transcriptomic responses to sublethal concentrations of the biotoxin okadaic acid

Specimen collection and experimental simulation of DSP HABs

Mussel specimens (Mytilus galloprovincialis (Lam.)) were collected in Valcobo beach, Galicia, NW Spain (43°19′02.71″N 8°21′56.35″W) in an area free of OA pollution during the resting period of the reproductive cycle (Banni et al., 2011). Sampled individuals (adults between 10 and 15 cm) were randomly divided into two experimental groups; exposed (exposed group) and non-exposed (control group) to the OA-producing dinoflagellate Prorocentrum lima. Both groups were kept in aerated seawater tanks and fed continuously with a suspension of the microalgae Tetraselmis suecica and Isochrysis galbana. After acclimation (one week), the exposed group was fed with 200 cells/L of a Prorocentrum lima culture (exponential phase) for 24 h. Specimens were dissected immediately after exposure, collecting samples from digestive gland and gill tissues. Each experimental sample consisted of tissue obtained from 5 individuals per group, dissected and pooled for RNA extraction.

RNA extraction and construction of cDNA libraries

The OA content in exposed and control samples was quantified using high-resolution mass spectrometry (Domenech et al., 2014). Total RNA was extracted from digestive gland and gill tissues using TRIzol^® (Thermo Scientific, Waltham, Massachusetts, USA) following the manufacturer’s instructions. RNA concentration and quality check was measured using a NanoDrop spectrophotometer (Thermo Scientific, Waltham, Massachusetts, USA) and a Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). cDNA library construction and pyrosequencing were performed using digestive gland samples, based on the larger absorption and accumulation of OA in this tissue. cDNA libraries were obtained from digestive gland tissue using the SMARTer™ PCR cDNA synthesis kit (Clontech, Mountain View, California, USA) and purified with GeneJET™ PCR Purification Kit (Thermo Scientific, Waltham, Massachusetts, USA) according to the manufacturer’s instructions.

The construction of normalized cDNA libraries (norm), for both exposed (mgt) and control (mgc) samples was carried out using the Trimer cDNA Normalization Kit (Evrogen, Moscow, Russia) following manufacturer’s protocol. This method enhances the detection of rare (lower concentration) transcripts by decreasing the prevalence of highly abundant transcripts (Bogdanova et al., 2011). The Suppression Subtractive Hybridization (SSH) libraries were constructed using the PCR-Select™ cDNA subtraction kit (Clontech, Mountain View, California, USA), following manufacturer’s instructions. Accordingly, two types of SSH libraries were produced: forward (fwd) and reverse (rev), representing upregulated and downregulated transcripts, respectively. This method was used to optimize the isolation of differentially expressed transcripts by removing commonly abundant cDNAs (Diatchenko et al., 1996).

Library sequencing and characterization

Normalized (exposed and control) and SSH (forward and reverse) libraries were sequenced by means of Roche-454 FLX + Titanium pyrosequencing (Roche Diagnostics, Indianapolis, Indiana, USA), with a sequencing depth of 40×. The obtained read datasets were preprocessed, assembled de novo and mapped to Gene Ontology (GO) and KEGG databases (Kanehisa, 2002). Additionally, low quality reads were discarded, and adaptors and low quality ends were trimmed before de novo assembly using MIRA v.3.9.16 (Chevreux, Wetter & Suhai Suhai, 1999). Both normalized and SSH read datasets are available at NCBI’s Bioproject database under the accession number PRJNA167773.

The generated contigs were annotated using BLAST (blastx) against the non-redundant protein sequence database (nr), setting a threshold e-value of 1e⁻⁶ (Altschul et al., 1997). Contigs were subsequently annotated with GO terms using the Blast2GO suite (Conesa et al., 2005; Gotz et al., 2008), including those terms obtained from InterPro and Annex analyses (Apweiler et al., 2001; Myhre et al., 2006).

Custom-made microarray construction and differential expression analysis

The sequencing and assembly of normalized and SSH libraries allowed to design specific probes targeting many of the transcripts identified. Accordingly, an Agilent oligonucleotide microarray encompassing 51,300 probes was constructed using the eArray™ design tool (Agilent Technologies, Santa Clara, California, USA) following a two-color Microarray-Based Gene Expression Analysis v.6.5 Agilent-specific protocol with dye swap. Two biological replicates per tissue sample were analyzed in microarray experiments. Expression analyses were conducted using the R package limma from the Bioconductor repository (Ritchie et al., 2015). Results are organized based on the magnitude of the observed change in expression or Fold Change in a logarithmic scale (logFC) and the statistical significance of the observed change in expression represented by an adjusted p-value or False Discovery Rate by the Benjamini–Hochberg procedure (FDR). Probes showing an FDR < 0.05 were considered as differentially expressed. The correlation between logFC values of differentially expressed transcripts commonly observed in both digestive gland and gill tissues was analyzed using a linear regression based on Pearson’s coefficient of determination. The GO terms for the most representative biological processes in both upregulated and downregulated groups of transcripts were determined using topGO with statistical significance (p-values) calculated according to the weight algorithm (Alexa & Rahnenfuhrer, 2010). Lastly, contigs were also mapped to the KEGG database for pathway analysis (Kanehisa, 2002).

Characterization of OA-specific cDNA libraries in the mussel Mytilus

The analysis of OA in pooled digestive gland tissue of exposed individuals revealed a concentration of 18.27 ng of OA per gram of fresh tissue in exposed individuals (OA content in controls individuals is below detection limit), an order of magnitude below the legal OA limit established for safe consumption of shellfish in the European Union (Reguera et al., 2014). This result reinforces the focus of the present study on early stages of DSP HAB episodes, at a moment when mussels start accumulating OA in their tissues but their commercialization is still allowed by law. The construction of normalized (norm) cDNA libraries yielded 919,177 good quality reads, 514,276 for the exposed group (mgt) and 404,901 for the control group (mgc). After assembly, a total of 24,624 and 16,395 consensus sequences (contigs) were obtained, respectively. Complementary, the SSH libraries produced a set of 1,221,928 good quality reads (SSH) with 469,795 corresponding to the forward (fwd) library and 752,133 to the reverse (rev) library. Once assembled, a total of 21,591 contigs and 33,437 contigs were obtained, respectively. Overall, blastx searches against the nr database resulted in the identification of 17,952 contigs from normalized libraries and 25,001 contigs from SSH libraries (see details in Table 1).

Given the high level of redundancy among de novo assembled libraries (Fig. 1), contigs were combined into unigenes according to their annotation and were considered equivalent to the annotated transcripts (unigenes are therefore considered a set of uniquely identified transcripts). The normalized and SSH libraries constructed expand and complement partial sequence data previously released by our group in the Chromevaloa database (Suarez-Ulloa et al., 2013a). By combining both sets of sequences, the present work was able to produce a microarray tool increasing the coverage of OA-specific transcriptome in the mussel Mytilus, improving the unbiased analysis of the differences in gene expression.

Microarray-based analysis of transcriptomic responses to OA

The present work expanded previous analysis of the mussel’s response to OA exposure using an omic approach using an oligonucleotide microarray designed from the sequences identified in pyrosequencing libraries. Accordingly, a medium-high coverage Agilent microarray (51,300 probes) was designed and developed using the sequences (contigs) obtained from the cDNA libraries constructed in this work. The hybridization of the microarray with RNA samples from exposed and control groups revealed a total number of 14,160 probes (digestive gland) and 6,913 probes (gill) differentially expressed (Fig. 2). The consistency between expression profiles in digestive gland and gill was assessed performing a linear regression of the logFC values of differentially expressed transcripts common for both tissues (i.e., those showing FDR < 0.05 in both cases), showing a good correlation between both sets of transcripts (Fig. 3). The detailed description of the transcripts displaying the highests differences in expression levels in both tissues, along with the maximum observed logFC value in the microarray analysis, is indicated in Supplemental Information 1 and 2.

The microarray analysis identified a set of transcripts displaying sharp expression differences between exposed and control treatments Supplemental Information 1 and 2, expanding the list of transcripts potentially involved in the response to OA (Manfrin et al., 2010; Suarez-Ulloa et al., 2013a). This was primarily facilitated by a larger coverage in the transcriptomic assessment, but also by the increase in bivalve genomic information that has been incorporated to molecular databases in recent years (Suarez-Ulloa et al., 2013b; Gerdol et al., 2014). Differentially expressed transcripts identified in this study include heat shock 70 kda protein 12b, proteases like cathepsins b and d, polyubiquitin and and proteasomal subunit beta type-4, commonly associated with an accumulation of misfolded or oxidized proteins observed under different types of environmental stress (Gotze et al., 2014). A subset of the identified transcripts showing the highest fold-change classified according to their main functional role is presented in Table 2. Our results corroborate previous analyses describing the responses of Mytilus galloprovincialis to OA stress (Manfrin et al., 2010), particularly the strong upregulation of vdg3 and elongation factor 2. In the case of vdg3, this gene is associated with developmental changes during the benthic settlement stage (He et al., 2014) and it has only been identified in bivalves, being particularly abundant in the digestive gland. On the other hand, the elongation factor 2 (EEF-2) is widely ubiquitous across eukaryotic taxa, playing an essential regulatory role in protein synthesis as a housekeeping gene. Although the functional implications of vdg3 and EEF-2 in the context of this study are still unclear, the obtained results support previous reports discouraging the use of EEF-2 as an internal control for qPCR analyses on bivalves without further validation (Du et al., 2013).

Opposite to these findings, a number of transcripts hypothesized to be responsive to OA failed to show substantial expression modifications under the conditions of this study. Notably, the Serine/Threonine phosphatases PP1 and PP2A, specific targets in OA toxicity mechanisms, did not show significant expression changes between treatments. OA is a well known selective inhibitor of the enzymatic activity of PP1 and PP2A phosphatases with critical consequences for the cell’s fate (Shenolikar, 1994). However, our results suggest that the upregulation of the PP1 and PP2A genes is not a relevant strategy versus the antagonist effects of OA. Similarly, Multi-Xenobiotic Resistance proteins (MXRs), good candidates to explain the high tolerance of bivalves versus pollution (Contardo-Jara, Pflugmacher & Wiegand, 2008), failed to show significant changes in expression. It is possible that their attributed role in OA uptake could be supplied by other proteins (e.g., the highly upregulated nose resistant to fluoxetine protein 6, a transport mediator of xenobiotics accross tissues). Indeed, lysosomal uptake has been suggested as a possible explanation for the extraordinary tolerance of mussels to the effects of DSP pollution (Svensson, Sarngren & Forlin, 2003).

In addition to transcripts previoulsy linked to OA responses, our results found an upregulation of an antimicrobial peptide (mytimacin) as well as an antifungal peptide (mytimycin) specific from mussels (Sonthi et al., 2011; Gerdol et al., 2012). Interestingly, mytimacin-5 (partial) was identified as one of the most upregulated transcripts in both gill and digestive gland. This peptide is especially interesting among the mytimacin family due to two additional cysteines in conserved positions predicted to form an extra disulfide bridge with yet unknown functional implications (Gerdol et al., 2012). C1q domain-containing proteins 1q3 and 1q25 showed a strong upregulation in the digestive gland. C1q is involved in the mammalian classical component pathway, playing an important role in innate immunity. Although no clear homologues to the vertebrate C1q complex subunits have been found in invertebrates yet, a massive expansion in the C1q domain-containing protein family has been suggested in bivalves, including Mytilus (Gerdol, Venier & Pallavicini, 2015). C1q domain-containing proteins are very versatile and might display a wide range of ligand interactions and functions such as clearance of apoptotic cells through direct binding (Kishore et al., 2004). They have been found upregulated in molluscs challenged with different pathogens (Perrigault, Tanguy & Allam, 2009; Taris et al., 2009). Although their specific function remains unclear, the substantial upregulation found in the present work might be indicative of a relevant role during environmental stress responses.

Altogether, the obtained results provide valuable insights into the molecular effects of OA in the mussel Mytilus and will be improved by considering the following: (a) kinetic effects during transcription and its regulation could bias estimations of differential expression (Bai et al., 2015); (b) the strong upregulation observed in endo-beta xylanases and endo-beta glucanases (although coherent with energy production) might be emphasized by the composition of the cell wall from dinoflagellates; (c) the differential regulation of antimicrobial and antifungal peptides might be influenced by the presence of infiltrated hemocytes in digestive gland tissue.

Expression and function profiles of transcripts differentially expressed in response to OA

The GO term annotation of transcripts differentially expressed in response to OA allowed the analysis of the biological processes in which their enconding genes are involved. A comparison of the functional profile for the two tissues studied is shown in Fig. 4. These profiles are based on the levels of representation for the most general sub-categories in GO stemming from Biological Process (Ashburner et al., 2000). Although absolute differences in magnitude are evident between digestive gland and gill (Fig. 2), no major functional differences were found when comparing the profiles for both tissues (Figs. 3 and 4). Nonetheless, such comparison might be hampered by sample size differences (e.g., subtle tissue-specific differences could remain undetected) and the fact that the microarray could lack gill-specific transcripts. Indeed, recent reports suggest that OA might display tissue-specific effects. Accordingly, different cytotoxic effects of OA specific for different human cell types had been demonstrated in vitro (Rubiolo et al., 2011). Furthermore, it has been reported that mussel gills display higher sensitivity to OA than hemocytes after one hour exposure (Prego-Faraldo et al., 2015). Tissue specificity is further evidenced by comparisons among enriched GO terms determined for transcripts upregulated and downregulated in digestive gland and gill (Table 3).

GO terms related with transcription regulation and cell cycle are enriched in the set of transcripts downregulated in the digestive gland (e.g., transcription from RNA polymerase II promoter, histone acetylation, mitotic nuclear division, mitotic S phase). On the contrary, these terms are mostly represented in the upregulated set of transcripts in the gill (e.g., positive regulation of cell growth, cell cycle, positive regulation of transcription, DNA-templated). Although this might suggest a higher degree of stress in digestive gland, both tissues consistently show enrichment in GO terms connected to DNA repair and degradation of damaged proteins. Therefore, while the mechanisms involved in OA toxicity could be consistent across these tissues, different responses could be elicited depending on the level of the accumulation. Further research will be required to clarify the extent in which the effects of OA are determined by the nature of the tissue, the time/dose or a combination of both.

Our results show an overall larger number of upregulated transcripts compared with those downregulated, in agreement with previous reports although a strong dependence of the expresion profiles with time was demonstrated (Manfrin et al., 2010). Such findings are further supported by the analysis of the response of the Pacific oyster Crassostrea gigas to OA exposure using time-series (Romero-Geraldo, Garcia-Lagunas & Hernandez-Saavedra, 2014), showing a strong dependence on time and dose. Altogether, it seems that expression profiles can hardly be extrapolated to other conditions different to those being studied. Given the highly dynamic nature of the transcriptome, only consistent patterns in expression can be informative of environmental stress conditions (Aardema & MacGregor, 2002). This supports the use of expression signatures rather than individual biomarkers for biomonitoring purposes. Modeling systems of greater complexity including time and dose as variables would provide valuable information about the dynamics of the expression profiles.

The present work was completed by investigating the metabolic pathways associated with those enzymes identified as differentially expressed under OA exposure conditions (Supplemental Information 3). Most of these pathways are involved in energy production (e.g., glycolysis/gluconeogenesis pathway, the citrate cycle, the pentose phosphate pathway and the oxidative phosphorylation pathway) as well as the regulation of the cell cycle and metabolism of drugs and xenobiotics. The observed functional profiles are consistent between tissues and also with observations in other organisms and types of abiotic stress. Accordingly, the role of metabolic functions was observed at the proteomic level in the mussel Perna viridis exposed to OA pollution (Huang et al., 2015). The differential expression of enzymes involved in metabolic pathways such as Glycolisis, TCA and oxidative phosphorylation suggests that energy production becomes critical in situations of environmental stress. Such observations agree with the responses found in the Eastern oyster Crassostrea virginica exposed to different types of abiotic stress (Chapman et al., 2011). Furthermore, the role of the mTOR pathway as key regulator of the balance between energy consumption and cellular development was also evidenced in bivalves under environmental stress (Clark et al., 2013). An upregulation of enzymes PI3K, AMPK, LKB1 and ERK1/2 from this pathway (responsible for arresting the cell cycle when energy is required for resisting stress conditions) was found in the present work, suggesting that such mechanism is activated in the mussel as a response to OA toxicity. Lastly, the differential expression of enzymes involved in immunity-related pathways like biosynthesis of antibiotics further supports a link between environmental stress and changes in the immunity system (Malagoli et al., 2007).