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

Background. Harmful Algal Blooms (HABs) responsible for Diarrhetic Shellfish Poisoning (DSP) represent a major threat for human consumers of shellfish. The biotoxin Okadaic Acid (OA), a well-known phosphatase inhibitor and tumor promoter, is the primary cause of acute DSP intoxications. Although several studies have described the molecular effects of high OA concentrations on sentinel organisms (e.g., bivalve molluscs), the effect of prolonged exposures to low (sublethal) OA concentrations is still unknown. In order to fill this gap, this work combines Next-Generation sequencing and custom-made microarray technologies to develop an unbiased characterization of the transcriptomic response of mussels during early stages of a DSP bloom. Methods. Mussel specimens were exposed to a HAB episode simulating an early stage DSP bloom (200 cells/L of the dinoflagellate Prorocentrum lima for 24 h). The unbiased characterization of the transcriptomic responses triggered by OA was carried out using two complementary methods of cDNA library preparation: normalized and Suppression Subtractive Hybridization (SSH). Libraries were sequenced and read datasets were mapped to Gene Ontology and KEGG databases. A custom-made oligonucleotide microarray was developed based on these data, completing the expression analysis of digestive gland and gill tissues. Results. Our findings show that exposure to sublethal concentrations of OA is enough to induce gene expression modifications in the mussel Mytilus. Transcriptomic analyses revealed an increase in proteasomal activity, molecular transport, cell cycle regulation, energy production and immune activity in mussels. Oppositely, a number of transcripts hypothesized to be responsive to OA (notably the Serine/Threonine phosphatases PP1 and PP2A) failed to show substantial modifications. Both digestive gland and gill tissues responded similarly to OA, although expression modifications were more dramatic in the former, supporting the choice of this tissue for future biomonitoring studies. Discussion. Exposure to OA concentrations within legal limits for safe consumption of shellfish is enough to disrupt important cellular processes in mussels, eliciting sharp transcriptional changes as a result. By combining the study of cDNA libraries and a custom-made OA-specific microarray, our work provides a comprehensive characterization of the OA-specific transcriptome, improving the accuracy of the analysis of expresion profiles compared to single-replicated RNA-seq methods. The combination of our data with related studies helps understanding the molecular mechanisms underlying molecular responses to DSP episodes in marine organisms, providing useful information to develop a new generation of tools for the monitoring of OA pollution.

modifications. Both digestive gland and gill tissues responded similarly to OA, although expression modifications were more dramatic in the former, supporting the choice of this tissue for future biomonitoring studies. Discussion. Exposure to OA concentrations within legal limits for safe consumption of shellfish is enough to disrupt important cellular processes in mussels, eliciting sharp transcriptional changes as a result. By combining the study of cDNA libraries and a custom-made OA-specific microarray, our work provides a comprehensive characterization of the OA-specific transcriptome, improving the accuracy of the analysis of expresion profiles compared to single-replicated RNA-seq methods. The combination of our data with related studies helps understanding the molecular mechanisms underlying molecular responses to DSP episodes in marine organisms, providing useful information to develop a new generation of tools for the monitoring of OA pollution.

INTRODUCTION
Harmful Algal Blooms (HABs) constitute an environmental phenomenon encompassing critical relevance due to their increasing frequency and impact in coastal areas (Anderson, 2009). Diarrhetic Shellfish Poisoning (DSP) blooms represent a major threat in widespread geographic areas comprising the Atlantic coast of Europe, Chile and Japan (Reguera et al., 2014), where natural outbreaks of toxic Dinophysis and Prorocentrum microalgae produce large amounts of DinophysisToXins (DTXs) and Okadaic Acid (OA) biotoxins (Sellner, Doucette & Kirkpatrick, 2003). OA is the primary cause of acute DSP intoxication of human consumers of shellfish, causing strong economic losses for the aquaculture industry. This biotoxin constitutes a well-known phosphatase inhibitor encompassing tumorigenic and apoptotic effects, even at low concentrations (Prego-Faraldo et al., 2015). Indeed, OA is capable of inducing genotoxic and cytotoxic damage, representing a hazard under chronic exposure conditions (Prego-Faraldo et al., 2013;Valdiglesias et al., 2013).
Given the noted risks of OA for human health and marine ecosystems, DSP events represent one of the most important threats for the shellfish aquaculture industry. Consequently, important efforts have been dedicated to develop rapid and sensible DSP biomonitoring methods, most notably using bivalve molluscs (e.g., mussels, oysters, clams, etc.) as sentinel organisms (Manfrin et al., 2010;Fernandez-Tajes et al., 2011;McNabb et al., 2012;Romero-Geraldo, Garcia-Lagunas & Hernandez-Saavedra, 2014;Huang et al., 2015). The choice of these organisms is supported by their wide distribution, sessile and filter-feeding lifestyles as well as their ability to accumulate high amounts of biotoxins, while displaying a particularly strong resilience to their harmful effects (Svensson, Sarngren & Forlin, 2003;Prado-Alvarez et al., 2012;Prado-Alvarez et al., 2013). During the last decade, the increasing availability of genomic resources in bivalves has improved classical biomonitoring approaches (e.g., quantification of biotoxin content in mollusc tissues), notably by developing molecular high-throughput studies evaluating omic (transcriptomic and proteomic) responses to HAB stress and their potential biomarker application (Manfrin et al., 2010;Suarez-Ulloa et al., 2013a;Gerdol et al., 2014;Huang et al., 2015). Nonetheless, while this approach has proven to be a promising venue for pollution biomonitoring (Campos et al., 2012;Suarez-Ulloa et al., 2013b), additional efforts are still required to clarify the cause-effect relationships between environmental stressors and changes in gene expression patterns. In doing so, it will be possible to transform the extraordinary amount of molecular data resulting from omic experiments into a practical tool for marine pollution biomonitoring.
Mussels start accumulating OA in their tissues during early stages of DSP blooms, however, their commercialization is still allowed by the applicable legislation as long as the concentration of this biotoxin does not exceed the legal threshold of 160 µg OA equivalents/kg shellfish meat (European Union legislation). Nonetheless, it has been demonstrated that exposure to low OA concentrations for short periods of time is enough to produce genotoxic and cytotoxic effects in vitro (Prego-Faraldo et al., 2015). The present work aims to provide a better understanding of the molecular mechanisms underlying the environmental responses of bivalve molluscs to sublethal concentrations of OA. For this purpose, Next-Generation sequencing and custom-made microarray technologies were combined to develop an unbiased characterization of the transcriptomic response of bivalve molluscs (mussels) to OA during early stages of a DSP bloom. These analyses build on previous studies (including our own) focused on specific subsets of genes (i.e., chromatin structure/function (Suarez-Ulloa et al., 2013a); oxidative stress, cell cycle regulation and immune response (Romero-Geraldo, Garcia-Lagunas & Hernandez-Saavedra, 2014;Romero-Geraldo & Hernandez-Saavedra, 2014)), as well as on the application of microarray technology to study the OA-specific transcriptome (Manfrin et al., 2010). Our results expand the scope, dimension and methodological approaches of these studies, improving the description of the cellular processes involved in the mussel response to OA toxicity. In doing so, this study generates omic information useful for identifying molecular signatures of marine pollution during DSP blooms. Contrary to quantitative analytical methods (i.e., LC-MS), this approach selectively identifies stressors of very different nature, assessing the magnitude of the toxic effects for organisms and communities. In addition, it provides further insights into the molecular strategies underlying the extraordinary resilience of bivalve molluscs to environmental stress.

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 TM PCR cDNA synthesis kit (Clontech, Mountain View, California, USA) and purified with GeneJET TM 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 TM 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.

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 TM 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 (pvalues) 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  (logFC) with statistical significance (FDR) indicated as a logarithmic scale. Probes highlighted in blue (FDR < 0.05) and purple (FDR < 0.05 and logFC > 2) represent the groups of transcripts displaying largest changes in gene expression between exposed and control treatments. 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 Figure 3 Correlation between paired logFC values calculated for transcripts identified in digestive gland and gill tissues between exposed and control treatments. Overall, a good level of agreement is found for gene expression changes (R 2 ∼ = 0.6). 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
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

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.

Table 3 Enriched GO terms in sets of differentially expressed transcripts in both digestive gland and gill tissues. Data is sorted based on p-value in increasing (p-values are calculated according to the weight algorithm in TopGO
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).

CONCLUSIONS
The present work dissects the gene expression changes in different mussel tissues during early stages of DSP HAB episodes, suggesting that low concentrations of OA (below the legal OA limit established for safe consumption of shellfish) are enough to elicit sharp changes in the expression of genes involved in the response to this biotoxin. Prior to this work, a few studies attempted to investigate the transcriptomic changes in bivalves during HABs using high-throughput methods (Manfrin et al., 2010;Suarez-Ulloa et al., 2013a;Gerdol et al., 2014). However, the combined application of normalized and SSH libraries together with the development of a custom-made OA-specific microarray in the present work, provides a more comprehensive characterization of the OA-specific transcriptome, improving the accuracy of the analysis of expresion profiles compared to single-replicated RNA-seq methods (Suarez-Ulloa et al., 2013a). The custom-made microarray platform generated in this work represents a convenient tool for long-term monitoring projects, offering a good level of standardization with lower requirements in computational resources comparing to the otherwise more informative RNA-seq methodology (Guo et al., 2013). In addition, the transcriptomic coverage of this microarray is comparable to recent estimations for the size of the complete transcriptome in digestive gland of Mytilus galloprovincialis (Gerdol et al., 2014) as well as for the transcriptome of the Pacific oyster Crassostrea gigas (Zhang et al., 2012), thus representing a good approximation to an unbiased tool for expression analysis.
Our results suggest that the response to OA found in mussels is consistent with the model of intracellular response to stress previously reported for bivalve molluscs (Anderson et al., 2015). Accordingly, the activation of energy production mechanisms observed in the present work could be producing potentially harmful Reactive Oxygen Species (ROS), which unless controlled by chaperones or eliminated in the proteasomes, would induce apoptosis. An increase in ROS production has been recently reported for the mussels exposed to saxitoxins (i.e., neurotoxins responsible for the paralytic shellfish poisoning), supporting the applicability of this model to HABs exposure (Astuya et al., 2015). Indeed, our results show an upregulation in important chaperones (Hsp70) and proteases (cathepsins b and d) ( Table 2) consistently with this model. Particularly the strong upregulation of cathepsins, known to be activated in the lysosomes (Kagedal, Johansson & Ollinger, 2001), in conjunction with the activation of transport mechanisms suggested by our results (Table 3), offer support to the lysosomal uptake hypothesis proposed by Svensson, Sarngren & Forlin (2003). In addition, the upregulation of antimicrobial peptides suggests the activation of immunity mechanisms in conjunction with the general environmental stress response. However, it remains unclear whether this immune response is automatically triggered by abiotic factors or whether there is an opportunistic attack of pathogens present in the microbiota of the mussels. Current efforts are directed to clarify this question (De Rijcke et al., 2015).
Further work studying more restricted conditions with shorter periods of exposure and lower concentrations of dinoflagellates would better inform about the sensitivity of the transcriptomic approach for the detection of OA-pollution in the ocean. Complementary, long-term monitoring projects in combination with meta-analysis of publicly available data could provide valuable information on the basal trancriptomic changes constituting a general environmental response as well as on the specific transcriptomic signature of DSP toxicity stress.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
The present work was supported by grants from the Biomolecular Sciences Institute (800005997) and start-up funds from the College of Arts and Sciences at Florida International University (JME-L), and from the Spanish Ministry of Economy and Competitivity (AGL2012-30897, JM). VS-U was supported by a Graduate Assistantship from the Department of Biological Sciences at FIU. VA-P was supported by the College of Engineering and Computing at Florida International University, under the supervision of Dr. Giri Narasimhan. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.