RNA-Seq Transcriptome Profiling of the Queen Scallop (Aequipecten opercularis) Digestive Gland after Exposure to Domoic Acid-Producing Pseudo-nitzschia.

Some species of the genus Pseudo-nitzschia produce the toxin domoic acid, which causes amnesic shellfish poisoning (ASP). Given that bivalve mollusks are filter feeders, they can accumulate these toxins in their tissues. To elucidate the transcriptional response of the queen scallop Aequipecten opercularis after exposure to domoic acid-producing Pseudo-nitzschia, the digestive gland transcriptome was de novo assembled using an Illumina HiSeq 2000 platform. Then, a differential gene expression analysis was performed. After the assembly, 142,137 unigenes were obtained, and a total of 10,144 genes were differentially expressed in the groups exposed to the toxin. Functional enrichment analysis found that 374 Pfam (protein families database) domains were significantly enriched. The C1q domain, the C-type lectin, the major facilitator superfamily, the immunoglobulin domain, and the cytochrome P450 were among the most enriched Pfam domains. Protein network analysis showed a small number of highly connected nodes involved in specific functions: proteasome components, mitochondrial ribosomal proteins, protein translocases of mitochondrial membranes, cytochromes P450, and glutathione S-transferases. The results suggest that exposure to domoic acid-producing organisms causes oxidative stress and mitochondrial dysfunction. The transcriptional response counteracts these effects with the up-regulation of genes coding for some mitochondrial proteins, proteasome components, and antioxidant enzymes (glutathione S-transferases, thioredoxins, glutaredoxins, and copper/zinc superoxide dismutases).


Introduction
The amnesic shellfish poisoning (ASP) toxin, domoic acid, is produced by some species of the genera Pseudo-nitzschia and Nitzschia [1,2]. The prevalence of domoic acid and toxic diatoms seems to have increased worldwide [1]. Domoic acid is a tricarboxylic amino acid that resembles glutamic acid and is a potent glutamate receptor agonist [3,4]. Bivalve mollusks are filter feeders, and during harmful algal blooms, they can accumulate toxins in their tissues. This is why they are the primary vectors for causing ASP in humans [1,5,6].
Domoic acid depuration time in bivalves is species-specific and can differ largely from one species to another [7]. Mussels of the genus Mytilus [8][9][10][11] and the oyster Crassostrea gigas [9] rapidly eliminate domoic acid, while the king scallop Pecten maximus [7] and the razor clam Siliqua patula [12] are very slow domoic acid depurators. In mussels and scallops, the digestive gland is the tissue with the highest domoic acid concentration [7,10,13,14]. Mauriz and Blanco [15] suggested that the very low depuration rate of P. maximus could be due to the lack of an efficient transmembrane transporter. Unlike the king scallop, in the queen scallop (Aequipecten opercularis) the depuration rate is fast [15].
Domoic acid is excitotoxic in the central nervous system of mammals and other vertebrates [3,4], but its putative effects on invertebrates have been less studied. Although it seems that domoic acid-producing organisms are not toxic to shellfish (or at least not highly toxic), they can exert several physiological and sublethal effects on marine bivalves [16][17][18][19][20]. Some of these effects include DNA damage in mussels [16], stress response characterized by shell closure, hemolymph acidosis, hypoxia, an increase in the number and activity of hemocytes in the oyster C. gigas [17,18], reduced larval growth in P. maximus [19], and negative impacts on growth rate and survival in juvenile king scallops [20]. Some authors have suggested that although several harmful algae toxins do not affect the survival of bivalve mollusks, they provoke oxidative stress [21][22][23]. However, the molecular mechanisms that cause oxidative stress are poorly understood, and furthermore, domoic acid causes oxidative stress in the vertebrate central nervous system [24][25][26][27]. The works cited above [16][17][18][19][20][21] analyzed several physiological and biochemical parameters after exposure to domoic acid but did not study the gene expression patterns in both exposed and non-exposed bivalves.
There are some publications regarding the gene expression changes associated with exposure to domoic acid in vertebrates [24,[28][29][30]. The transcriptome response was dependent on the dose and the exposure duration; among the differentially expressed genes were those involved in transcription (transcription factors), signal transduction, ion transport, generalized stress response, mitochondrial function, inflammatory response, DNA damage, apoptosis, neurological function, and neuroprotection [24,[28][29][30]. In a previous work, we studied (by RNA-seq) the effects of domoic acid-containing Pseudo-nitzschia on gene expression in the mussel Mytilus galloprovincialis [31], and to our knowledge, this is the only published work about the transcriptional effects of domoic acid exposure on mollusks. As stated before, among the bivalves there are large interspecific differences in the domoic acid depuration rate [7][8][9][10][11][12]. It is therefore necessary to carry out gene expression studies on different species.
Understanding the molecular mechanisms of domoic acid uptake and elimination in bivalves and how the toxin (and the toxin-producing species) affects gene expression are two knowledge gaps in this field. The aim of the present work is to contribute to filling these gaps by means of a transcriptomic approach. First, the whole transcriptome of the A. opercularis digestive gland was de novo assembled, and then, we analyzed by RNA-seq the transcriptional changes after exposure to domoic acid-producing Pseudo-nitzschia. This approach can provide some clues regarding the biological and molecular processes altered by domoic acid. The transcriptomic approach has been successfully employed to uncover the genetic response of bivalves to diarrhetic shellfish poisoning (DSP) and paralytic shellfish poisoning (PSP) toxins and also to identify the genes putatively involved in detoxification processes [32][33][34][35][36].

Domoic Acid Content in the Digestive Gland of A. opercularis
A group of six scallops sampled on April 9 from the tank (group DB) had an average domoic acid content of 1361 ± 804 ng/g digestive gland wet weight, while a group of six scallops sampled on April 17 from the raft (group DA) had an average domoic acid content of 6680 ± 1661 ng/g digestive gland wet weight (Table 1). In the six control scallops sampled on May 12 from the tank (group C), the domoic acid levels were below the limit of quantification (BLOQ, Table 1). The total scallop wet weights and digestive gland (DG) wet weights are shown in Table 1 and Table S1. Table 1. Domoic acid concentration (ng/g digestive gland wet weight), wet weight (g) of the soft tissues (Total weight), and wet weight (g) of the digestive gland (DG weight) in sampled scallops (Aequipecten opercularis).

Group
Sampling

Sequencing and de novo Assembly
After the de novo assembly, the transcripts were clustered (homology >90%) to reduce redundancy. Thus, 142,137 unigenes were obtained ( Table 2). The minimum, maximum, and mean contig lengths were 200, 17,867, and 1343.9 bp, while the N50 contig length was 1845 bp ( Table 2). The raw data are accessible from the NCBI Sequence Read Archive (Project PRJNA508885, sample accession numbers from SAMN10537388 to SAMN10537405). Table 2.
Summary of Illumina transcriptome sequencing and assembly for A. opercularis digestive glands.

Differential Expression, Functional Annotation, and Functional Enrichment Analysis
A total of 26,932 and 20,608 differentially expressed genes (DEGs) were detected in group DA and in group DB, respectively, when compared to the control (C) group ( Figure 1). Genes that were differentially expressed in both groups (the groups that had accumulated domoic acid) and in the same direction (either down-or up-regulation) were selected for further study: 10,144 genes, including 4913 up-regulated and 5231 down-regulated ( Figure 1; Tables 3 and 4; Files S1 and S2). The top 25 significantly up-regulated genes are listed in Table 3. Genes coding for fatty acid-binding proteins and cytosolic sulfotransferases were among the top up-regulated genes ( proteins and cytosolic sulfotransferases were among the top up-regulated genes ( Table 3). Most of  the top 25 down-regulated genes do not have a Blastx hit (Table 4; File S2). Figure 1. Scheme of the differential expression results obtained in A. opercularis digestive glands after exposure to domoic acid-producing Pseudo-nitzschia.
A summary of the functional annotation results is shown in Table 5. After functional enrichment performed using the Pfam annotations, 374 domains were found to be significantly (false discovery rate (FDR)-adjusted p-value <0.05) enriched in the DEGs (Table 6; File S3). The C-type lectin, the RNA recognition motifs, the major facilitator superfamily, and the cytochrome P450 were among the most enriched Pfam domains whose genes were mostly up-regulated (Table 6). In addition to the cytochromes P-450, several Pfam domains involved in biotransformation (phase I and phase II metabolism of xenobiotics) were functionally enriched and up-regulated: glutathione S-transferases, sulfotransferases, methyltransferases, and aldehyde dehydrogenases (Table 6; File S3). By contrast, most of the genes coding for proteins with C1q domains, immunoglobulin domains, tetratricopeptide repeats, and collagen triple helix repeats were down-regulated (Table 6; File S3).
Significantly enriched gene ontology (GO) terms (Fisher's exact test, FDR-adjusted p-value <0.05) in the biological process (BP), molecular function (MF), and cellular component (CC) categories are displayed in Table 7 (up-regulated DEGs),  (Table 7) metabolic process, oxidation-reduction process, and organic substance catabolic process (in the BP category); catalytic activity, oxidoreductase activity and threonine-type peptidase activity (in the MF category); and cytoplasm, proteasome complex, and endopeptidase complex (in the CC category).
On the other hand, the number of enriched GO terms for the down-regulated DEGs (File S4) were 86 (BP), 111 (MF), and 32 (CC). Table 8 shows that the top enriched GO terms were as follows: neurotransmitter transport, regulation of cellular process, and cell communication (in the BP category); neurotransmitter transporter activity, neurotransmitter/sodium symporter activity, and solute/sodium symporter activity (in the MF category); and transcription factor complex, collagen trimer, and cytoskeleton (in the MF category). A summary of the functional annotation results is shown in Table 5. After functional enrichment performed using the Pfam annotations, 374 domains were found to be significantly (false discovery rate (FDR)-adjusted p-value <0.05) enriched in the DEGs (Table 6; File S3). The C-type lectin, the RNA recognition motifs, the major facilitator superfamily, and the cytochrome P450 were among the most enriched Pfam domains whose genes were mostly up-regulated (Table 6). In addition to the cytochromes P-450, several Pfam domains involved in biotransformation (phase I and phase II metabolism of xenobiotics) were functionally enriched and up-regulated: glutathione S-transferases, sulfotransferases, methyltransferases, and aldehyde dehydrogenases (Table 6; File S3). By contrast, most of the genes coding for proteins with C1q domains, immunoglobulin domains, tetratricopeptide repeats, and collagen triple helix repeats were down-regulated (Table 6; File S3).
Significantly enriched gene ontology (GO) terms (Fisher's exact test, FDR-adjusted p-value <0.05) in the biological process (BP), molecular function (MF), and cellular component (CC) categories are displayed in Table 7 (up-regulated DEGs),  (Table 7) metabolic process, oxidation-reduction process, and organic substance catabolic process (in the BP category); catalytic activity, oxidoreductase activity and threonine-type peptidase activity (in the MF category); and cytoplasm, proteasome complex, and endopeptidase complex (in the CC category).
On the other hand, the number of enriched GO terms for the down-regulated DEGs (File S4) were 86 (BP), 111 (MF), and 32 (CC). Table 8 shows that the top enriched GO terms were as follows: neurotransmitter transport, regulation of cellular process, and cell communication (in the BP category); neurotransmitter transporter activity, neurotransmitter/sodium symporter activity, and solute/sodium symporter activity (in the MF category); and transcription factor complex, collagen trimer, and cytoskeleton (in the MF category).     Among the level-2 enriched GO terms ( Figure 2; Table S2), the genes in the categories of metabolic process, cellular process, catalytic activity, structural molecule activity, and transporter activity were mainly up-regulated, while most of the genes in the categories of biological regulation, signaling, immune system process, response to stimulus, and transcription regulator activity were down-regulated. File S5 shows the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologues (KO) of DEGs and of all unigenes. Among the level-2 enriched GO terms ( Figure 2; Table S2), the genes in the categories of metabolic process, cellular process, catalytic activity, structural molecule activity, and transporter activity were mainly up-regulated, while most of the genes in the categories of biological regulation, signaling, immune system process, response to stimulus, and transcription regulator activity were down-regulated. File S5 shows the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologues (KO) of DEGs and of all unigenes.

Protein Network Analysis
Protein-protein interactions can be employed to group and organize all the protein-coding genes in a genome [37]. From the 4913 up-regulated genes, a Blastx search found 931 human homologs in the STRING database. The network obtained in the highest confidence (0.9) mode is enriched in interactions (p-value <1 × 10 −16 ). The results obtained with the up-regulated DEGs showed a small number of highly connected protein nodes. Each group of proteins is involved in specific biological processes (Figure 3; Figure S1; File S6): degradation of proteins (proteasome components), synthesis of mitochondrial proteins (mitochondrial ribosomal proteins), translocation of cytosolically synthesized mitochondrial preproteins (translocases of outer and inner mitochondrial membrane, TOMMs, and TIMMs), splicing of mRNA (spliceosome components), and phase I and phase II metabolism of xenobiotics (cytochromes P450 and glutathione S-transferases).

Protein Network Analysis
Protein-protein interactions can be employed to group and organize all the protein-coding genes in a genome [37]. From the 4913 up-regulated genes, a Blastx search found 931 human homologs in the STRING database. The network obtained in the highest confidence (0.9) mode is enriched in interactions (p-value <1 × 10 −16 ). The results obtained with the up-regulated DEGs showed a small number of highly connected protein nodes. Each group of proteins is involved in specific biological processes ( Figure 3; Figure S1; File S6): degradation of proteins (proteasome components), synthesis of mitochondrial proteins (mitochondrial ribosomal proteins), translocation of cytosolically synthesized mitochondrial preproteins (translocases of outer and inner mitochondrial membrane, TOMMs, and TIMMs), splicing of mRNA (spliceosome components), and phase I and phase II metabolism of xenobiotics (cytochromes P450 and glutathione S-transferases).
From the 5231 down-regulated DEGs, 855 human homologs were found in the STRING database. The network is enriched in interactions (p-value <1 × 10 −16 ). Components of different types of collagen, heat shock proteins, and proteins involved in cytoskeleton dynamics ( Figure S2; File S7) were among the proteins that appeared in the network obtained with the down-regulated DEGs. Network was constructed using the String 10.5 algorithm and obtained in the highest confidence (0.9) mode. Some highly connected protein nodes are highlighted. Proteins were named according to the human protein name. A full list of protein names is available in Figure S1 and File S6.
The normalized gene expression of the five target genes, and the non-selected reference gene (MYH9) is displayed in Figure 5. There was good agreement between RT-qPCR ( Figure 5, upper panel) and RNA-seq ( Figure 5, lower panel). The RNA-seq results showed that CYP2C14, SLC16A12, ANT1, and SLC16A13 were up-regulated in groups DB and DA in relation to the control; these genes were also up-regulated when the RT-qPCR data were analyzed. The SLC6A9 gene was down-regulated in groups DB and DA in relation to the control group ( Figure 5, lower panel), but the RT-qPCR data showed significant differences only between group DA and the control. The candidate reference gene MYH9 is not differentially expressed. Network was constructed using the String 10.5 algorithm and obtained in the highest confidence (0.9) mode. Some highly connected protein nodes are highlighted. Proteins were named according to the human protein name. A full list of protein names is available in Figure S1 and File S6.
The normalized gene expression of the five target genes, and the non-selected reference gene (MYH9) is displayed in Figure 5. There was good agreement between RT-qPCR ( Figure 5, upper panel) and RNA-seq ( Figure 5, lower panel). The RNA-seq results showed that CYP2C14, SLC16A12, ANT1, and SLC16A13 were up-regulated in groups DB and DA in relation to the control; these genes were also up-regulated when the RT-qPCR data were analyzed. The SLC6A9 gene was down-regulated in groups DB and DA in relation to the control group ( Figure 5, lower panel), but the RT-qPCR data showed significant differences only between group DA and the control. The candidate reference gene MYH9 is not differentially expressed.   Figure 4. Determination of the optimal number of reference genes for normalization. The pairwise variation (Vn/n+1) was calculated between the normalization factors NFn and NFn+1 (using n or n+1 reference genes respectively) by geNorm software.

Discussion
There are many studies about the mechanisms of the neurotoxicity of domoic acid in vertebrates (especially mammals), but there is very little knowledge about the putative effects of domoic acid on bivalve mollusks. Several publications showed that domoic acid can exert physiological and sublethal effects on marine bivalves [16][17][18][19][20]. Dizer et al. [16] found that in Mytilus edulis, DNA damage was significantly increased after the injection of domoic acid and suggested the existence of genotoxic responses in the cells of digestive glands. It is interesting to point out that DNA repair, cellular response to DNA damage stimulus, and cellular response to stress were among the enriched GO terms in the present work for the up-regulated DEGs (File S4). This could be the transcriptomic response of A. opercularis to the putative DNA damage provoked by domoic acid. In C. gigas, domoic acid provoked a generalized stress response [17] and an increase in the number and activity of hemocytes [18]. Domoic acid induces oxidative stress in the central nervous system and spinal cord in vertebrates [24][25][26][27]. Furthermore, harmful algae toxins sometimes provoke oxidative stress in bivalves [21][22][23], and Prego-Faraldo et al. [40] found that exposure to the toxic dinoflagellate Prorocentrum lima induces the differential expression of genes coding for antioxidant enzymes. The glutathione S-transferase, thioredoxin, glutaredoxin, and copper/zinc superoxide dismutase Pfam domains were functionally enriched in queen scallops (File S3), and these genes were, for the most part, up-regulated; these domains are found in proteins involved in protection against reactive oxygen species (ROS).
In a recent work about the effects of environmental stress on gene transcription in oysters, Anderson et al. [41] proposed a consensus model of sub-cellular stress responses in oysters with the involvement of mitochondria and reactive oxygen species (ROS) production. If the anti-oxidant enzymes and molecular chaperones cannot limit the damage caused by ROS, then the consequences are probably cellular dysfunction and apoptosis [41]. In vertebrates, domoic acid causes mitochondrial dysfunction as a consequence of oxidative stress [4,24,26]. Hiolski et al. [24] suggested the existence of compensatory mitochondrial biogenesis in response to mitochondrial dysfunction. Our results (Figure 3; Figure S1; File S2,S6) showed an up-regulation of genes coding for mitochondrial ribosomal proteins and translocases of the outer and inner mitochondrial membrane (proteins involved in mitochondrial biogenesis); in the protein interaction network, these proteins form highly connected protein nodes (Figure 3; Figure S1) Up-regulation of proteasome subunits (Figure 3; Figure S1; File S2) is also a possible consequence of oxidative stress [42], because the proteasome is responsible for the selective degradation of oxidized proteins [43]. The 26S proteasome is a protease complex, which is responsible for the regulated degradation of proteins in eukaryotic organisms [42]. The proteasome system can be activated to accomplish the destruction of proteins altered by stress conditions [44]. The proteasome complex and proteasome core complex were two of the most enriched GO terms in the cellular component category for the up-regulated DEGs (Table 7), and the proteasome Pfam domain (PF00227) was also enriched (File S3); the genes coding for proteins with this domain were up-regulated in Pseudo-nitzschia-exposed scallops (Files S2 and S3). Proteasome proteins form a group of highly connected nodes in the protein-protein interaction network (Figure 3, Figure S1). In M. galloprovincialis exposed to the toxin okadaic acid, there is an up-regulation of several mRNAs involved in proteasome activity [45]. Therefore, the results support the hypothesis that exposure to domoic acid-producing Pseudo-nitzschia causes oxidative stress and the impairment of the mitochondrial function in A. opercularis and that the transcriptional changes are directed, at least in part, to counteract the stress effects.
The metabolism of xenobiotics (such as toxins) has three phases; phase I (functionalization) and phase II (conjugation) are catalyzed by metabolizing enzymes, while phase III consists of the export from the cell by transmembrane transporter proteins. We found that the Pfam domains of some phase I (cytochromes P450 and aldo-keto reductases) and phase II (glutathione S-transferases and sulfotransferases) drug metabolizing enzymes were functionally enriched, and the genes coding for these enzymes were mostly up-regulated. Cytochromes P450 and glutathione S-transferases constituted a group of highly connected nodes in the protein interaction network (Figure 3, Figure S1). Genes of these families were also up-regulated in mussels (M. galloprovincialis) exposed to domoic acid-containing Pseudo-nitzschia [31]. Peña-Llopis et al. [46] showed that the treatment of scallops (P. maximus) with N-acetylcysteine increased glutathione S-transferase (GST) activity and enabled the scallops to eliminate domoic acid more efficiently. Therefore, it is possible that glutathione S-transferases play a role in domoic acid detoxification. Li et al. [35] found an up-regulation of sulfotransferase genes in the kidney of the Zhikong scallop Chlamys farreri after exposure to paralytic shellfish toxin-producing Alexandrium minutum. Furthermore, the family of sulfotransferases is significantly expanded in the C. farreri genome [35]; the high number of transcripts coding for sulfotransferases in the A. opercularis transcriptome is indicative of an expansion of this family in the queen scallop.
The molecular mechanisms of domoic acid absorption and excretion in bivalve mollusks are poorly understood [31]. Mauriz and Blanco [15] found that in the king scallop P. maximus, domoic acid is free in the cytosol of the digestive gland and suggested that the low depuration rate of domoic acid in this species could be due to the lack of membrane transporters. Domoic acid is a charged compound and probably needs a transport protein to pass through the plasma membrane, as is the case with glutamic acid [47,48]. This putative transmembrane transporter(s) could therefore play an important role in the absorption and/or the excretion of domoic acid. The results of Kimura et al. [47] suggest that anion exchange transporters are responsible for the transmembrane transport of domoic acid in Caco-2 cell monolayers (which represent the intestinal barrier of mammals); these transporters belong to the solute carrier (SLC) superfamily. In A opercularis, we found that transmembrane transport and transmembrane transporter activity were two of the enriched GO terms (File S4); furthermore, the major facilitator superfamily (MFS) was one of the most significantly enriched Pfam families ( Table 6). Most of the genes belonging to these categories were up-regulated (Table 6; Files S2-S4). MFS is a clan of the SLC superfamily [49], and Hediger et al. [50] reported that the SLC gene series included 52 families in the human genome, although it has recently been updated to 65 families [51]. A total of eight SLC families (SLC5, SLC16, SLC17, SLC21, SLC22, SLC26, SLC39, and SLC49) contain up-regulated genes in A. opercularis (File S8) and in M. galloprovincialis [31] exposed to domoic acid-containing Pseudo-nitzschia. The transporter protein(s) putatively involved in the uptake and/or elimination of domoic acid in the digestive gland of bivalve mollusks could be encoded by a gene from one of those families. The families with a higher number of up-regulated genes in both A. opercularis and M. galloprovincialis were SLC16 (the monocarboxylate transporters family) and SLC22 (organic cation/anion/zwitterion transporters). A total of four members of the human SLC16 gene family encode monocarboxylate transporters, but the substrates of several members are unknown [52]. The SLC22 family [53] comprises organic cation, zwitterion, and anion transporters (OCTs, OCTNs, and OATs), which participate in the absorption (in the small intestine) and excretion (in the liver and kidney) of xenobiotics and endogenous substances [53]. Unfortunately, the lack of knowledge about the identity of the domoic acid transmembrane transporter(s) in mammals makes it difficult to identify them in bivalve mollusks. Schultz et al. [54] suggested that ATP-binding cassette (ABC) transporters are responsible for the absorption of domoic acid in Dungeness crabs, but in A. opercularis, we found only five up-regulated ABC transporters. A similar result was reported by Pazos et al. [31] in M. galloprovincialis (with two up-regulated ABC transporters). This contrasts with the high number of up-regulated SLC genes found in both bivalves.
Although most of the SLC genes differentially expressed in A. opercularis were up-regulated, the SLC6 family (the sodium-and chloride-dependent neurotransmitter transporter family) is an exception, with 48 down-regulated unigenes (File S8). Furthermore, neurotransmitter/sodium symporter activity is one of the most enriched GO terms for the down-regulated genes (Table 8). On the contrary, in M. galloprovincialis, two genes from this family were up-regulated and none down-regulated [31]. In A. opercularis, the number of genes in this family is very high, and this agrees with results from Li et al. [35], who found that the SLC6 family is expanded in the C. farreri (a scallop) genome in relation to other bivalves.
Harmful algae and biotoxins exert different effects on the immune systems of bivalve mollusks [23,32,55]. Immune response and immune system process were two of the most enriched GO terms for the down-regulated DEGs (File S4), and Pfam families involved in immunological processes were significantly enriched in A. opercularis (Table 6; File S3): the C1q domain-containing proteins, the C-type lectin, the fibrinogen beta and gamma chains C-terminal globular domain, the immunoglobulin domain, and tumor necrosis factors (TNF). Except for the C-type lectins, the genes in these categories were mainly down-regulated (Table 6; Files S2 and S3). Differentially expressed genes from these families have been found in several bivalve mollusks after exposure to different biotoxins [31][32][33][34]36,56,57]. Hégaret et al. [23] found that some harmful algae provoked a stimulation of immune function of bivalve hemocytes, while others were immunosuppressive. The C1q domain containing proteins and the C-type lectins are particularly abundant in the digestive glands of bivalves [58,59]. The C1q domain-containing proteins are indispensable in the innate immune systems of invertebrates [60] and could be involved in several functions, such as activation of the complement pathway, cell adhesion, pathogen recognition, response to pollutants, and apoptosis [58,60,61]. An expansion of the genes coding for proteins containing the C1q domain was found in several bivalves [58,[61][62][63]. For example, 321 C1q domain-containing proteins are encoded by the C gigas genome [62]; this represents approximately 10-fold more than the Ciq proteins encoded by the Homo sapiens genome [62]. Some genes coding for C1q domain-containing proteins were down-regulated in M. galloprovincialis fed with toxigenic strains of Alexandrium minutum [34]. The C-type lectins are characterized by a calcium-dependent carbohydrate recognition domain and participate in pathogen recognition and in innate immunity in bivalves [59], but they can also perform non-immune functions; for example, a role in efficient food particle sorting (food recognition) was found in the oyster Crassostrea virginica [64]. There is a high number of genes coding for C-type lectins in bivalve mollusks [59,65]. Most of the genes coding for C-type lectins were up-regulated in A opercularis (Table 6; Files S2 and S3) and M. galloprovincialis [31] after exposure to domoic acid-producing Pseudo-nitzschia; this agrees with the up-regulation of C-type lectins in M. chilensis after exposure to saxitoxin [33,56,57]. On the contrary, these genes were down-regulated in Argopecten irradians in response to okadaic acid [32].
One of the most enriched GO terms for the down-regulated DEGs in A. opercularis was collagen trimer (Table 8), and collagen triple helix repeat was among the enriched Pfam domains (Table 6). Furthermore, several collagen components form a group of highly connected protein nodes in the network obtained with the down-regulated genes ( Figure S2). In M. galloprovincialis, after exposure to domoic acid-containing Pseudo-nitzschia [31], some collagen genes (7) were down-regulated, although the number of induced genes was greater (13). Collagens are components of the extracellular matrix characterized by the presence of at least one triple-helical domain [66]. They are among the most abundant proteins and have mainly a structural function [66].
Another group of predominantly down-regulated genes (14 up-regulated and 35 down-regulated, File S2) were those coding for heat shock proteins (HSPs). Half of the induced HSP genes were mitochondrial forms, and among the repressed ones, the HSP70 genes predominated. Heat shock proteins are involved in protein folding and can be induced by several types of stress, including high temperature, toxins, pathogens, and hypoxia [67]. Several publications have reported the increased expression of heat shock protein genes in bivalves after exposure to harmful algae toxins [32,57,[67][68][69]. Cheng et al. [67] found an expansion of Hsp70 (heat shock protein 70 kDa) genes from the Hspa12 subfamily in Mizuhopecten yessoensis. Several of these genes were differentially expressed in response to Alexandrium catenella exposure (most of them were induced, but there were also some Hsp70 genes down-regulated [67]). However, Ryan et al. [29] reported the down-regulation of Hsp68 (a member of the Hsp70 family) after domoic acid exposure in mouse brain. Furthermore, the exposure to domoic acid-producing Pseudo-nitzschia provoked a down-regulation of HSPs in M. galloprovincialis [31]; these results were coincident with those obtained in the present work.
Another result worth highlighting is that some genes coding for glutamate ionotropic receptors, two genes coding for NMDA receptors, and five coding for kainate (KA) receptors, were all down-regulated in the present study (File S2). The zebrafish gria2 gene (glutamate ionotropic receptor AMPA 2) was down-regulated after two weeks of low-level domoic acid exposure [24], and the authors suggested that this down-regulation is a compensatory response to elevated glutamatergic activity [24]. It is possible that a similar compensatory mechanism takes place in queen scallops after exposure to domoic acid.
The RT-qPCR results confirm the differential gene expression obtained by RNA-seq. The gene expression changes and expression levels (fold change in relation to control) assessed by the two methods ( Figure 5) were very similar. In the determination of gene expression by means of RT-qPCR, the validation of the reference genes for each experimental situation is very important [70,71]. The utilization of RNA-seq expression data allowed us to find more suitable candidate reference genes. Thanks to this, their stability values (Table 10), calculated by geNorm and NormFinder, were low (which means that the expression was stable). The selected reference genes performed better than some traditional reference genes, such as 18S rRNA, ACTB, and EF1A [70,71].

Conclusions
RNA-seq technology was employed to elucidate the transcriptional response triggered by exposure to domoic acid-producing Pseudo-nitzschia in the queen scallop A. opercularis. A total of 10,144 genes were differentially expressed in the two toxin-exposed groups of scallops in relation to the control group (4913 up-regulated and 5231 down-regulated).
The results obtained are compatible with the hypothesis that exposure to domoic acid-producing Pseudo-nitzschia causes oxidative stress in A. opercularis. Some consequences of oxidative stress are the impairment of mitochondrial function and oxidation of proteins; therefore, the transcriptional response of the queen scallop tries to counteract these effects with the up-regulation of genes coding for proteins involved in the following: degradation of oxidized proteins (proteasome components), mitochondrial biogenesis (mitochondrial ribosomal proteins, TOMMs, and TIMMs) and antioxidant enzymatic activity (glutathione S-transferases, thioredoxins, glutaredoxins, and copper/zinc superoxide dismutases). The results of the present work and those cited in the literature show that oxidative stress is one of the most common effects of the exposure to toxins and toxin-producing algae, and a part of the harmful effects of the toxins are due to oxidative stress.
A great number of up-regulated genes code for proteins involved in the metabolism of xenobiotics (cytochromes P450, aldo-keto reductases, glutathione S-transferases, and sulfotransferases) and transmembrane transport (solute carriers), while the genes coding for proteins with domains involved in immunological processes (C1q domain, C-type lectins, immunoglobulin domain, fibrinogen beta and gamma chains C-terminal globular domain, and tumor necrosis factors) were mainly down-regulated, with the exception of the C-type lectins.

Materials and Methods
The methods employed were the same as those previously described [31] except for minor modifications.

Animals
Queen scallops (A. opercularis) were obtained from a natural bed in the Ría de Arousa in December 2014 and maintained in a 500-L tank, in the Centro de Investigacións Mariñas, (CIMA, Pedras de Corón, Vilanova de Arousa, Spain), with a continuous unfiltered seawater flow (approximate) of 1200 L/h. On April 9, 2015, 2 random samples of the scallops were obtained, and the remaining scallops (control, group C) were maintained in the tank. The scallops in 1 of the samples (group DB) were analyzed to determine their individual content and concentration of domoic acid. The scallops in the other sample (group DA) were placed in culture baskets and transferred to a raft in the culture area Grove C2 in the Ría de Arousa, where a bloom of Pseudo-nitzschia was taking place. The recorded levels of domoic acid in the mussels from that raft showed a maximum of 22 mg DA/kg on April 9 (data obtained from Intecmar, www.intecmar.gal [72]). The scallops were maintained on the raft until April 17, 2015 in order to be exposed to domoic acid-containing Pseudo-nitzschia. On that date, the scallops (group DA) were brought back to the laboratory to determine their domoic acid content. The scallops from the control group were sampled on May 12, 2015, after the end of the toxic episode caused by Pseudo-nitzschia. From April 9 to May 12, the main characteristics of the seawater, temperature, salinity, light transmission (index of suspended solids), O 2 , and fluorescence (index of phytoplankton abundance) in GROVE and in the area of CIMA were very similar ( Figure S3). As previously explained by Pazos et al. [31], the experimental approach (animals naturally exposed to domoic acid-producing Pseudo-nitzschia) was chosen because of the difficulty of supplying toxic Pseudo-nitzschia under controlled conditions in the laboratory due to the relatively low absorption efficiency of the scallops and to the loss of toxicity of the Pseudo-nitzschia cultures.
Digestive glands, gills, and the remaining tissues were obtained by dissecting the scallops. Then, 1 part of each digestive gland was used in the determination of the domoic acid content. The second part was stored in RNAlater (ref. AM7021, Ambion, Life Technologies, Carlsbad, CA, USA) at −80 • C until RNA extraction.

Chemicals and Reagents for Toxin Extraction and LC-MS/MS
Methanol for HPLC and formic acid were purchased from RCI Labscan Limited (Bangkok, Thailand) and Sigma-Aldrich (St. Louis, MO, USA), respectively. Ultrapure water was obtained using a Milli-Q Gradient system, coupled with an Elix Advantage 10, both from Millipore (Merck Millipore, Darmstadt, Germany).

Determination of the Domoic Acid Content
To extract the toxin, each digestive gland was placed in aqueous methanol (50%) in a proportion of 1:2 (w/v) and homogenized with an Ultra-Turrax T25 (IKA, Staufen, Germany). The extract was clarified using centrifugation at 18,000 g at 4 • C for 10 min, retaining a supernatant that was immediately analyzed.
Domoic acid in the obtained extracts was analyzed using LC-MS/MS. The chromatographic separation was carried out using a Thermo Accela chromatographic system (Thermo Fisher Scientific, Waltham, MA, USA), with a high-pressure pump and autosampler. The stationary phase was a solid core Kinetex C18, 50 × 2.1 mm, 2.6 µm particle size, column (Phenomenex, Torrance, CA, USA). An elution gradient, with a flow of 280 µL/min, was used with mobile phase A (formic acid 0.2%) and B (50% MeOH with formic acid 0.2%). The gradient started at 100% A, maintained this condition for 1 min, linearly changed until reaching 55% B in minute 5, was held for 2 min, and then reverted to the initial conditions in order to equilibrate before the next injection. Next, 5 µL of extract, previously filtered through a PES 0.2-µm syringe filter (MFS), were injected.
After the chromatographic separation, domoic acid was detected and quantified by means of a Thermo TSQ Quantum Access MAX triple quadrupole mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA), equipped with a HESI-II electrospray interface, using positive polarization and SRM mode. The transition 312.18 > 266.18 m/z was used to quantify the response and 312.18 > 248.18 was used for confirmation. The spectrometer was operated under the following conditions: spray voltage 3400 V, capillary temperature 270 • C, HESI-II temperature 110 • C, sheet gas (nitrogen) 20 (nominal pressure), auxiliary gas (nitrogen) 10 (nominal pressure), collision energy of 15 V, and collision gas (argon) pressure of 1.5 mTorr.
Concentrations of domoic acid were obtained by comparing the response of the quantification transition in the sample extracts with that of a reference solution obtained from NRC Canada. The quantification limit of the method for tissue analysis is less than 20 ng/mL of extract.

RNA Extraction
For digestive gland total RNA isolation, the NucleoSpin RNA kit (ref. 740955, Macherey-Nagel, Düren, Germany) was used following the manufacturer's protocol. Then, RNA was precipitated with 0.5 volumes of Li CL 7.5 M, and the RNA pellet was dissolved in 50 µL of RNA storage solution (ref. AM7000, Ambion, Life Technologies, Carlsbad, CA, USA). To remove DNA contamination, total RNA was treated with DNA-free (ref. AM1907M, Ambion, Life Technologies, Carlsbad, CA, USA). The integrity and quality of the RNA samples were measured using agarose gel electrophoresis, an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and a Nanodrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The quantity of the total RNA was determined using Qubit 2.0 (Invitrogen, Carlsbad, CA, USA).

Library Preparation and Sequencing
A total of 18 cDNA libraries were generated ( Figure 6) from the digestive gland of the scallops (6 obtained from each group: DB, DA, and control). The poly(A) + mRNA fraction was isolated from total RNA, and cDNA libraries were obtained following Illumina's recommendations. Briefly, poly(A) + RNA was isolated on poly-T oligoattached magnetic beads and chemically fragmented prior to reverse transcription and cDNA generation. The cDNA fragments then went through an end repair process, the addition of a single 'A' base to the 3' end, and afterwards, the ligation of the adapters.

Library Preparation and Sequencing
A total of 18 cDNA libraries were generated ( Figure 6) from the digestive gland of the scallops (6 obtained from each group: DB, DA, and control). The poly(A) + mRNA fraction was isolated from total RNA, and cDNA libraries were obtained following Illumina's recommendations. Briefly, poly(A) + RNA was isolated on poly-T oligoattached magnetic beads and chemically fragmented prior to reverse transcription and cDNA generation. The cDNA fragments then went through an end repair process, the addition of a single 'A' base to the 3' end, and afterwards, the ligation of the adapters. Finally, the products were purified and enriched with PCR to create the indexed final double-stranded cDNA library. The quality of the libraries was analyzed using a Bioanalyzer 2100 high sensitivity assay; the quantity of the libraries was determined by real-time PCR in a LightCycler

de novo Assembly
Quality control checks of the raw sequencing data were performed with FastQC. The technical adapters were eliminated using Trimgalore software version 0.3.3 (Babraham Bioinformatics, Cambridge, UK) (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Additionally, Figure 6. Scheme of the methods employed for sequencing and de novo transcriptome assembly.

de novo Assembly
Quality control checks of the raw sequencing data were performed with FastQC. The technical adapters were eliminated using Trimgalore software version 0.3.3 (Babraham Bioinformatics, Cambridge, UK) (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Additionally, the reads with a mean Phred score >30 were selected. Subsequently, all the samples were combined, and the complexity of the reads was reduced by removing duplicates. Then, a de novo assembly was performed using the programs Oases, version 0.2.09 [73] and Trinity, version 2.1.1 [74]. The assembled transcripts were clustered (>90% homology) to reduce redundancy using cd-hit software version 4.6. For each sequence, the potential ORFs were detected using Transdecoder software, version 2.0, with standard parameters.
Each sample was then mapped with Bowtie2, version 2.2.6 [75] against the reference transcriptome obtained in the previous step. The good quality reads (Mapping Quality ≥20) were selected to increase the resolution of the count expression. Finally, the expression inference was evaluated by means of the counts of properly paired reads in each transcript.

Differential Expression
The transcriptome expression for each sample was normalized by library size, following the DESeq2 protocols. Considering the whole normalized transcriptome, a study of correlation and Euclidean distance between samples was performed using the statistical software R, version 3.2.3 (www.r-project.org), for identifying possible samples outliers.
Differential gene expression analysis was performed with DESeq2 algorithm, version 1.8.2 (http://www.bioconductor.org/packages/devel/bioc/html/DESeq2.html). The genes with a fold change of less than −2 or greater than 2 and a p-value adjusted using the Benjamini and Hochberg [76] method for controlling false discovery rate (FDR) <0.05 were considered differentially expressed (Figure 7).  [74]. The assembled transcripts were clustered (>90% homology) to reduce redundancy using cd-hit software version 4.6. For each sequence, the potential ORFs were detected using Transdecoder software, version 2.0, with standard parameters. Each sample was then mapped with Bowtie2, version 2.2.6 [75] against the reference transcriptome obtained in the previous step. The good quality reads (Mapping Quality ≥20) were selected to increase the resolution of the count expression. Finally, the expression inference was evaluated by means of the counts of properly paired reads in each transcript.

Differential Expression
The transcriptome expression for each sample was normalized by library size, following the DESeq2 protocols. Considering the whole normalized transcriptome, a study of correlation and Euclidean distance between samples was performed using the statistical software R, version 3.2.3 (www.r-project.org), for identifying possible samples outliers.
Differential gene expression analysis was performed with DESeq2 algorithm, version 1.8.2 (http://www.bioconductor.org/packages/devel/bioc/html/DESeq2.html). The genes with a fold change of less than −2 or greater than 2 and a p-value adjusted using the Benjamini and Hochberg [76] method for controlling false discovery rate (FDR) <0.05 were considered differentially expressed (Figure 7).    The contigs that had a lower E-value versus Pseudo-nitzschia compared with M. yessoensis were discarded (approximately 0.7% of the sequences).

Functional Annotation
The genes were annotated using Blastx [77] against the Uniprot database and Blastn [77] against the NCBI nucleotide database (E-value threshold of 10 −2 ). Then, the annotation was expanded by incorporating information from the species, gene name, and functions using gene ontology and protein structure domains associated with the transcript using InterPro (https://www.ebi.ac.uk/interpro/). The genes were also annotated with Blast2GO software version 4.1.9 (BioBam Bioinformatics S.L, Valencia, Spain) [78,79] Orthologue assignment and pathway mapping were performed on the KEGG Automatic Annotation Server (KAAS, [80]) using Blast and the bi-directional best hit (BBH) method (http://www.genome.jp/tools/kaas/).

Functional Enrichment
A functional enrichment study was performed using the Pfam [81] functional information. This study is based on hypergeometric distribution [82] using the statistical software R version 3.2.3 (www.r-project.org). The differentially expressed genes were also subjected to GO enrichment analysis with Blast2GO version 4.1.9. (BioBam Bioinformatics S.L, Valencia, Spain) using Fisher's exact test [83] (up-and down-regulated genes were analyzed separately). The false discovery rate (FDR) adjusted p-value [76] was set at a cutoff of 0.05.

Protein Network Analysis
To search for the protein-protein interactions, network analyses using the String 10.5 algorithm [84] were performed. The putative human homologues of proteins coded by the up-regulated and the down-regulated genes in the A. opercularis digestive gland were identified by means of a Blastx search [85] against the STRING human protein database (9606.protein.sequences.v10.fa), with an E-value threshold of 10 −5 . The top Blastx search results were used as input in the String program. The up-regulated and the down-regulated genes were analyzed separately. For the relative quantification of gene expression by means of RT-qPCR, a normalization step must be performed using internal reference genes, whose expression levels are stable [38,[86][87][88]. Suitable reference genes should be selected for each experimental condition to ensure their stable expression [70,89].
A total of 6 reference gene candidates (Table 9), VAMP7, RPS4, MYH9, EIF4EBP2, DNAJ, and RAP1B and 5 target genes (Table 9), CYP2C14, SLC16A12, ANT1, SLC16A13, and SLC6A9, were used in the gene expression study. The candidate reference genes were selected for their stable expression based on the RNA-seq data. Oligonucleotide primers were designed with OligoAnalyzer 3.1 (http://eu.idtdna.com/analyzer/Applications/OligoAnalyzer/; Integrated DNA Technologies, Leuven, Belgium) from the sequences in Table 9 and were synthesized by Integrated DNA Technologies (Leuven, Belgium). The primer sequences and amplicon lengths are listed in Table 9. The specificity of the primers was confirmed by the presence of a single peak in the melting curve and by the presence of a single band of the expected size when PCR products were run in a 2% agarose gel. The PCR amplification efficiency (E) of each transcript was determined by means of Real-Time PCR Miner software (Version 4.0; http://www.miner.ewindup.info/ [90]). The mean amplification efficiency (E) of each amplicon (Table 9)  The cycling conditions were as follows: 30 s at 95 • C (initial template denaturation) and 40 cycles of 5 s at 95 • C (denaturation) followed by 10 s at 60 • C (annealing and elongation) and 10 s at 75 • C for fluorescence measurement. At the end of each run, a melting curve was carried out: 95 • C for 20 s and 60 • C for 20 s followed by an increase in temperature from 60 to 100 • C (with temperature increases in steps of 0.5 • C every 10 s). Baseline values were automatically determined for all the plates using Bio-Rad iCycler iQ software V3.1 (IQ™ Real-Time PCR Detection System). The threshold value was set manually at 100 RFU (relative fluorescence units) to calculate the Cq values. Non-reverse transcriptase controls and non-template controls (NTC) were also included in each run.
The gene expression was normalized to reference genes that had stable expression levels [38,[86][87][88]. The gene expression stability of candidate reference genes was analyzed using 3 Microsoft Excel-based software applications, geNorm V3.5 [38], NormFinder V0.953 [86], and BestKeeper V1 [88]. The non-normalized expression (Q) was calculated using the equation Q = (1 + E) -Cq . Then, the expression was normalized by dividing it by the normalization factor (the geometric mean of the non-normalized expression of the selected reference genes) [89].
The statistical analyses were performed with the IBM SPSS Statistics 24.0 package (IBM SPSS, Chicago, IL, USA). The data were tested for normality (Shapiro-Wilk test) and for homogeneity of variance (Levene's test). The gene expression was log-transformed (base 2) to meet the requirements of normality and homogeneity of variances. The expression of target genes in domoic acid-exposed scallops (groups DB and DA) in relation to the control group was compared using ANOVA and post hoc Dunnett's t test. p <0.05 was considered statistically significant.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6651/11/2/97/s1, Table S1: Wet weight and domoic acid content of the queen scallops (A. opercularis) in the three groups of the study; Table S2: List of level-2 enriched gene ontology (GO) terms for differentially expressed genes in biological process (BP), molecular function (MF), and cellular component (CC) categories; Figure S1: Network showing interactions (confidence view) of proteins coded by genes up-regulated in the present study; Figure S2: Network showing interactions (confidence view) of proteins coded by genes down-regulated in the present study; Figure S3: Fluorescence (relative units), dissolved oxygen (mL/L), salinity, temperature ( • C), and light transmission (%) of the seawater between 1-and 5-m depth between April 9 and May 12, in the area of CIMA (laboratory), and GROVE (transplanted scallops); File S1: Nucleotide sequences of differentially expressed genes (in fasta format); File S2: List of differentially expressed genes in groups DA and DB (in relation to the control group); File S3: Significantly enriched Pfam families among the differentially expressed genes; File S4: Significantly enriched GO terms; File S5: List of KO (KEGG Orthologues) for the differentially expressed genes and for all the genes; File S6: Results of a Blastx search of up-regulated genes against the STRING human protein database (9606.protein.sequences.v10.fa), and list of input proteins in STRING network analysis; File S7: Results of a Blastx search of down-regulated genes against the STRING human protein database (9606.protein.sequences.v10.fa), and list of input proteins in STRING network analysis; and File S8: List of the differentially expressed genes belonging to the solute carriers (SLC) superfamily.