Transcriptional Response in the Digestive Gland of the King Scallop (Pecten maximus) After the Injection of Domoic Acid

Some diatom species of the genus Pseudo-nitzschia produce the toxin domoic acid. The depuration rate of domoic acid in Pecten maximus is very low; for this reason, king scallops generally contain high levels of domoic acid in their tissues. A transcriptomic approach was used to identify the genes differentially expressed in the P. maximus digestive gland after the injection of domoic acid. The differential expression analysis found 535 differentially expressed genes (226 up-regulated and 309 down-regulated). Protein–protein interaction networks obtained with the up-regulated genes were enriched in gene ontology terms, such as vesicle-mediated transport, response to stress, signal transduction, immune system process, RNA metabolic process, and autophagy, while networks obtained with the down-regulated genes were enriched in gene ontology terms, such as response to stress, immune system process, ribosome biogenesis, signal transduction, and mRNA processing. Genes that code for cytochrome P450 enzymes, glutathione S-transferase theta-1, glutamine synthase, pyrroline-5-carboxylate reductase 2, and sodium- and chloride-dependent glycine transporter 1 were among the up-regulated genes. Therefore, a stress response at the level of gene expression, that could be caused by the domoic acid injection, was evidenced by the alteration of several biological, cellular, and molecular processes.


Introduction
Some diatom species of the genus Pseudo-nitzschia produce domoic acid, a toxin that can cause amnesic shellfish poisoning (ASP) in humans [1][2][3][4][5]. During harmful algae blooms the bivalves accumulate the toxins in their tissues and therefore they can act as vectors of ASP [1][2][3]. The accumulation of biotoxins in shellfish can cause harvesting closures and thus it has adverse economic impacts. In recent years there has been an increase in the number of toxigenic Pseudo-nitzschia blooms worldwide [2,6]. Domoic acid, an amino acid structurally similar to glutamate and kainic acid, is a glutamate receptor agonist that binds mainly to two subtypes of ionotropic receptors (α-amino-3-hydroxy-5-methylisoxazole-4-propionate, AMPA, and kainate receptors [4,7,8]) and exerts excitotoxic effects in the central nervous system of vertebrates [4,7,8].
The king scallop Pecten maximus is a valuable fisheries resource in Europe [9]. Unlike mussels and oysters [10][11][12][13], with high domoic acid depuration rates, in the king scallop (Pecten maximus) the depuration rate of domoic acid is very low [14,15]. Due to the blooms of domoic acid-producing species of the genus Pseudo-nitzschia and the low depuration rate of P. maximus [14,15], the domoic acid concentration in this scallop is usually above the regulatory limits (20 µg of domoic acid g −1 ) in many areas [16]. There are other scallop species (Placopecten magellanicus and Argopecten purpuratus) that do not show this slow domoic acid depuration [17,18].
Domoic acid in bivalves is mostly unmetabolized and is excreted unchanged [10]. In many bivalves, including the king scallop, the digestive gland is the main organ of accumulation of domoic acid [10,14,15,17,19,20]. In the digestive gland of P. maximus, domoic acid was found to be present in the cytosol in soluble form [21], therefore Mauriz and Blanco [21] suggested that the lack of an efficient membrane transporter could be the cause of the low depuration rate. In a recent work, Blanco et al. [16] showed that this toxin, in the digestive gland, is mainly accumulated in large cells (digestive cells) and the concentration was lower in small cells (secretory cells).
Several publications have described the physiological effects of domoic acid and domoic acid-producing organisms on marine bivalves: transient DNA damage in Mytilus edulis after de injection of domoic acid [22], mild respiratory alkalosis in Mytilus californianus after the exposure to domoic acid-producing Pseudo-nitzschia [23], transient respiratory acidosis, hypoxia, and increase in the activity and number of haemocytes in Crassostrea gigas following the exposure to domoic acid-producing Peseudo-nitzschia [23,24], negative impacts on survival and growth rate in juvenile Pecten maximus exposed to domoic acidspiked feed [25], negative effects on growth and survival in Pecten maximus larvae exposed to dissolved domoic acid [26], impairment of immune functions and oxidative stress in Argopecten irradians exposed to dissolved domoic acid [27,28]. In two previous works, we studied the transcriptional effects of domoic acid-containing Pseudo-nitzschia in the digestive gland of the mussel Mytilus galloprovincialis [29] and the queen scallop Aequipecten opercularis [30] by means of RNA-seq technology. The results obtained showed the differential expression of genes involved in protection against oxidative stress, the metabolism of xenobiotics (detoxification), transmembrane transport, and immunological processes [29,30]. Oxidative stress is one of the principal effects caused by exposure to domoic acid in vertebrates [31][32][33][34], and also in invertebrates [35,36], including marine bivalves (Argopecten irradians [28] and Aequipecten opecularis [30]), and can play an important role in domoic acid induced toxicity. In the present work, a transcriptomic approach (RNA-seq) has been used to identify the genes differentially expressed in the digestive gland of the king scallop P. maximus after the injection of domoic acid, in order to uncover dysregulated biological and molecular processes and contribute to the knowledge of detoxification mechanisms in bivalves. Table 1 and Figure 1 display the domoic acid concentration (µg g −1 wet weight) and burden (µg) in different P. maximus tissues. Due to the low depuration rate of P. maximus [14,15], and the blooms of toxigenic Pseudo-nitzschia, the king scallops collected in Galicia (NW Spain) have high levels of domoic acid in their body tissues. This explains why domoic acid concentration in the digestive gland (DG) is high in control and treated scallops (2222.9 and 1987.0 µg g −1 wet weight, respectively). The domoic acid burden in the digestive gland of control and treated animals was 3690.2 µg and 3884.8 µg, respectively.

Domoic Acid Content in the Tissues of P. maximus
The domoic acid content in the digestive gland accounted for 84% of the total domoic acid burden in the body tissues of P. maximus, despite the fact that the digestive gland only represented about 7% of the total weight of soft tissues. The kidney (Table 1 and Figure 1) is the only organ that shows a clear increase in the concentration of domoic acid (and in the domoic acid burden) in treated animals (1226.34 µg g −1 wet weight) in relation to the control (370.70 µg g −1 wet weight). In absolute terms, the greatest increases in the domoic acid burden were found in the kidney followed by the digestive gland (Table 1). Significant differences between control and treated animals were found only in the concentration of domoic acid in the kidney (log10 transformed to meet the assumptions of normality and homoscedasticity; t-test, p = 0.041).
Although the experiment was not designed to see the effect of domoic acid on the behavior of scallops, we have observed that the scallops treated with domoic acid maintained a larger valve opening than those injected with seawater, but no mortality was recorded.

Sequencing and De Novo Assembly
The reads obtained by means of paired-end sequencing in the Illumina Hiseq 2000 platform were de novo assembled with Trinity and Oasis. Then, to reduce redundancy, the assembled transcripts were clustered (homology > 90%), thus 72,673 unigenes were obtained (Table 2). Mean contig length and N50 contig length were 1481 bp and 1883 bp respectively ( Table 2). The raw data are accessible from the NCBI Sequence Read Archive (BioProject ID PRJNA704533, BioSample accessions: SAMN18043529 to SAMN18043540). A BUSCO analysis of completeness of the de novo assembly identified 82.29% complete orthologs (45.6% single-copy and 36.69% duplicated), 3.88% fragmented orthologs and 13.83% missing orthologs (Table S1). The BUSCO results indicate that our de novo transcriptome assembly is of high quality, and although we only sequenced one tissue (digestive gland) the assembly has a completeness of 86.2%.

Differential Expression and Functional Annotation
A representation of the differentially expressed genes (DEGs) is displayed in a MA plot (Supplementary Figure S1) and a volcano plot (Supplementary Figure S2). The number of genes with absolute fold change > 1.5 and adjusted p-value < 0.05 (DEGs) was 535, of which 226 were up-regulated and 309 down-regulated (Supplementary Files S1 and S2). Tables 3 and 4 display the top 20 significantly up-regulated genes and the top 20 significantly down-regulated genes, respectively.  Among the top up-regulated genes ( Table 3, Supplementary File S1) were genes involved in the metabolism of xenobiotics (cytochromes P450, glutathione S-transferase) a gene (SLC6A9) coding for a sodium-and chloride-dependent glycine transporter (a solute carrier of the SLC6 family [37]) and two genes coding for enzymes involved in the metabolism of glutamate and proline (glutamine synthetase and pyrroline-5-carboxylate reductase 2). Pinin, serine protease inhibitor Cvsi-1, innexin UNC-7, serine/threonineprotein kinase PINK1 and apoptosis-stimulating of p53 protein 1 were among the top down-regulated genes ( Table 4). Many of the top down-regulated genes code for uncharacterized proteins (Table 4).
It is worth pointing out that we have found 19 transcripts that code for putative glutamate receptors in the digestive gland of P. maximus, 10 of them for ionotropic receptors and 9 for metabotropic receptors (Table 5 and Supplementary File S3), but none of these genes were differentially expressed.  The functional annotation results are shown in Table 6 and in Supplementary Files S4 and S5.

Protein Network Analysis
Differentially expressed protein-coding genes can be grouped by means of the proteinprotein interactions [38]. A Blastx search of the sequences of the 226 up-regulated and the 309 down-regulated genes found 88 and 93 human homologs in the STRING database, respectively. The networks obtained with up-regulated ( Figure 2 and Supplementary File S6) and down-regulated ( Figure 3 and Supplementary File S7) genes were enriched in interac-tions (p-values 5.31 × 10 −6 and 0.0152, respectively). Up-regulated genes ( Figure 2 and Supplementary File S6) were involved in biological processes, such as the metabolism of xenobiotics, stress response, immune response, lysosomal degradation, autophagy, apoptosis, RNA processing, and exocytosis. On the other hand, down-regulated genes ( Figure 3 and Supplementary File S7) were involved in biological processes such as RNA processing, ribosome formation, apoptosis, immune and inflammatory responses.
The up-regulated genes that code for proteins that showed interactions in the protein network analysis were enriched in gene ontology (GO) terms (Supplementary File S8) such as vesicle-mediated transport, response to stress, signal transduction, immune system process, RNA metabolic process, autophagy, lysosome and oxidoreductase activity. The down-regulated genes that code for proteins that showed interactions in STRING analysis were enriched in GO terms (Supplementary File S8) such as response to stress, immune system process, ribosome biogenesis, signal transduction, mRNA processing, and oxidoreductase activity.

Real Time RT-qPCR
Four candidate reference genes (GAPDH, EF1A, COX1 and NDUFA7, Table 7) were selected based on a previous work on P. maximus [39]. Table 8 shows the results obtained with the three algorithms (geNorm, NormFinder and BestKeeper) employed to test the stable expression of these genes. GAPDH, and COX1 were the best ranked genes and therefore selected as reference genes.
The comparison of gene expression results of the six target genes (MRP7, CYP2B4, P5CR, SLC6A9, FERRITIN, CYP2U1, Table 7) obtained with RT-qPCR and RNA-seq is displayed in Figure 4. The log2 FC values achieved by the two methods were very similar and showed good correlation (r = 0.969; r 2 = 0.939).

Discussion
Several authors have shown [14][15][16] that the digestive gland in the king scallop P. maximus accumulates most of the domoic acid. This agrees with the results obtained in the present work that found that domoic acid content in the digestive gland accounts for 84% of the total domoic acid burden in the body tissues of P. maximus ( Table 1). The total amount of domoic acid injected was 3000 µg per animal, but only a mean difference of 339 µg between treated and control animals was found. These results suggest that part of the administered domoic acid was depurated. After the injection into the adductor muscle, the domoic acid enters into the circulatory system, through the hemolymph sinuses [40]. In scallops, urine formation takes place through hemolymph filtration into the pericardial cavity, which in turn is connected to the kidneys through the reno-pericardial ducts [40]. Therefore, part of the domoic acid is probably filtered from the hemolymph into the kidneys and only a part of the injected domoic acid was distributed to the digestive gland. In kidneys, we have found a significant difference in domoic acid concentration between control and treated animals (Figure 1). It is important to emphasize that in the natural environment the scallops obtain the toxin from the food-chain through the digestive system.
We found 535 DEGs (226 up-regulated and 309 down-regulated, Supplementary File S1) in the digestive gland of P. maximus after the injection of domoic acid, therefore the toxin might have an effect on the gene expression in the digestive gland. The analysis of gene expression found the alteration of some processes at the biological, cellular, and molecular level (Figures 2 and 3, Supplementary File S8), that could be due to the effects of domoic acid. Thus, genes involved in xenobiotic metabolism, immune response, response to stress, signal transduction, apoptosis, RNA processing, ribosome biogenesis, lysosomal degradation, autophagy, and regulated exocytosis were differentially expressed. Scallops injected with domoic acid showed a behavior different (maintained a clearly larger valve opening) to the ones to which only seawater was injected. This reveals that domoic acid probably had an effect at the level of the central nervous system. There are three ways, not mutually exclusive, that can explain how domoic acid exerts its effects on the digestive gland of P. maximus: (a) the interaction of domoic acid with different biomolecules (mainly proteins) after the entry of domoic acid into the cells of the digestive gland; (b) the binding of domoic acid to glutamate receptors present on the plasma membrane of the cells of the digestive gland (we found that there is expression of mRNA that codes for glutamate receptors, both ionotropic and metabotropic, in the digestive gland of P. maximus, Table 5 and Supplementary File S3); and (c) the binding of domoic acid to glutamate receptors present in the central nervous system (cerebral, pedal, and parietovisceral ganglia) of the scallop (the actions triggered in these ganglia could be transmitted to the digestive gland through the nervous and neuroendocrine systems). Dizer et al. [22] also found an effect (increased DNA damage) of the intramuscular injection of domoic acid on digestive gland cells of a bivalve, Mytilus edulis.
In vertebrates, domoic acid is a potent neurotoxin [4,8,[41][42][43], and the response to domoic acid includes genes involved in transcription (transcription factors), signal transduction, ion transport, generalized response to stress, mitochondrial function, inflammatory response, response to DNA damage, apoptosis, neurological function and neuroprotection [31,41,44,45]. Although there are fewer studies on the effects of domoic acid on invertebrates than vertebrates, this toxin also exerts effects on marine bivalves at the physiological and gene expression levels [22][23][24][25][26][28][29][30]. In two previous studies [29,30] we showed that exposure to domoic acid containing Pseudo-nitzschia alters the transcriptomic profile of the digestive gland of the mussel Mytilus galloprovincialis and of the queen scallop Aequipecten opercularis. The results obtained by Ventoso et al. [30] suggest that exposure to domoic acid-producing organisms causes oxidative stress and mitochondrial dysfunction in A. opercularis, thus the transcriptional response of the queen scallop is involved in the protection against oxidative stress. This agrees with the results obtained by Song et al. [28] that showed that domoic acid induces oxidative stress and impairs defence mechanisms in bay scallops (Argopecten irradians). The contribution of oxidative stress to the effects and toxicity of domoic acid has been highlighted by several authors [6,28,[31][32][33]35,36].
One of the effects of domoic acid in P. maximus was the down-regulation of genes involved in RNA processing, ion transport, immune response, metabolic process and signal transduction (Supplementary File S8); this agrees with the results of Lefebvre et al. [41] with zebrafish, after low-level domoic acid exposures, that found the down-regulation of genes involved in those same biological processes.
Genes coding for several phase I (cytochromes P450) and phase II (glutathione Stransferases and sulfotransferases) drug metabolizing enzymes were up-regulated in P. maximus (Table 3, Figure 2 and Supplementary File S1), these types of genes were also up-regulated in A. opercularis [30] and M. galloprovincialis [29] following exposure to domoic acid-producing Pseudo-nitzschia.
Several authors have shown that glutamate receptors are expressed not only in the central nervous system but also in other types of tissues or organs (intestine, liver, kidney, stomach, etc.) [51][52][53]. Therefore, glutamate and glutamate receptor agonists could participate in the regulation of several physiological processes in peripheral organs [51][52][53]. We have found 19 genes that code for possible glutamate receptors in the digestive gland of P. maximus, 10 of them for ionotropic receptors and 9 for metabotropic receptors (Table 5 and Supplementary File S3). Some of the effects of domoic acid on the cells of the digestive gland may be mediated by these receptors. None of the genes coding for these receptors were differentially expressed in P. maximus (Table 5 and Supplementary File S3). In A. opercularis, some genes coding for glutamate ionotropic receptors were down-regulated [30] in the digestive gland of animals exposed to domoic acid-containing Pseudo-nitzschia. This may be due to a compensatory response to elevated glutamatergic activity, thus Hiolski et al. [31] found this type of compensatory response in zebrafish after domoic acid exposure.
Glycine, in addition to acting as an inhibitory neurotransmitter, is also a co-agonist at N-methyl-D-aspartate (NMDA) glutamate receptors [54]. In the central nervous system of vertebrates, the glycine transporter 1 (sodium-and chloride-dependent glycine transporter 1) regulates the binding of glycine to NMDA receptors [54], because the action of glycine is terminated through the reuptake mediated by sodium-and chloride-dependent glycine transporters [55]. The up-regulation of the SLC6A9 gene (coding for sodium and chloride-dependent glycine transporter 1) could prevent or reduce NMDA receptor activation. The SLC6A9 gene was among the top up-regulated genes in P. maximus (Table 3). There was another gene of this family (SLC6) that was downregulated in P.maximus (Supplementary File S1). Although both genes code for putative sodium-and chloridedependent glycine transporters, they share only 52% sequence identity at the amino acid level. Genes of this family (SLC6) were up-regulated in M. galloprovincialis [29] and downregulated in A. opercularis [30] after exposure to domoic acid-producing Pseudo-nitzschia. A gene of the SLC6 family was up-regulated in Pseudo-nitzschia multiseries under toxinproducing conditions [56], and this gene was also up-regulated in a domoic acid-producing Pseudo-nitzschia species in relation to two Pseudo-nitzschia species that do not produce domoic acid [57].
The SLC6 family is expanded in the genome of the scallops Chlamys farreri and Patinopecten yessoensis [58,59], in relation to other bivalves. In the A. opercularis [30] and in P. maximus digestive gland transcriptome, the number of transcripts belonging to this family is also very high (we found 58 in P. maximus).
One of the up-regulated genes in P. maximus, glutamine synthetase (Supplementary File S1), might play a neuroprotective role against glutamate neurotoxicity in neural tissues [60,61], because it catalyzes the transformation of glutamate to glutamine. Glutamate and glutamate receptor agonists increased glutamine synthetase expression and glutamine synthetase activity in cultured astrocytes [62,63]. Glutamine synthetase also participates in the production of GABA (gamma-aminobutyric acid), an inhibitory neurotransmitter. GABA has been shown to be able to prevent, at least partially, the effects of domoic acid in rat glial cells [64]. Therefore, the overexpression of this gene could have a protective effect against domoic acid. Another gene involved in the metabolism of amino acids (glutamate and proline) is up-regulated in P. maximus. This gene codes for the enzyme pyrroline-5-carboxylate reductase 2 that catalyzes the conversion of pyrroline-5-carboxylate to proline, and proline has a protective effect against oxidative stress [65]. Kenny et al. [66] suggested that mutations in the sodium channel gene, Neuron Navigator 1 (Nav1), could protect against the effects of domoic acid in P. maximus. We have not found transcripts of the Nav1 gene in the digestive gland of P. maximus, therefore the sodium channel is likely expressed in nervous tissue but not in the digestive gland.
The immune system of marine bivalves is sensitive to toxins and harmful algae [67,68]. Several harmful algae can provoke a stimulation of immune functions, while others cause inhibition [68]. Chi et al. [27] found that exposure to domoic acid impaired immune functions in the bay scallop A. irradians. In P. maximus, the immune system process was one of the enriched terms in the proteins coded by DEGs that showed interaction in STRING (Figures 2 and 3, Supplementary Files S6-S8). C-type lectins are calcium-dependent carbohydrate-binding proteins and participate in innate immunity in bivalves [67,69]. In bivalve mollusks there is a high number of genes coding for C-type lectins [69,70]. Three transcripts coding for putative C-type lectins were up-regulated in the P. maximus digestive gland (Supplementary File S1). This agrees with the results obtained with A. opercularis [30] and M. galloprovincialis [29] after exposure to domoic acid-producing Pseudo-nitzschia, where most of the genes coding for C-type lectins were up-regulated.
Heat shock proteins (HSP) can be induced by several types of stress (high temperature, hypoxia, toxins, or pathogens) and they are involved in protein folding [71]. There is an expansion of Hsp70 (heat shock protein 70 kDa) genes from the Hspa12 subfamily in Mizuhopecten yessoensis [71], and a gene of this subfamily (heat shock 70 kDa protein 12A-like) was up-regulated in P. maximus (Supplementary File S1). A up-regulation of HSP genes after exposure to harmful algae toxins (including domoic acid) has been found in several bivalves [27,28,[71][72][73][74][75], although in A. opercularis, after exposure to domoic acid-producing Pseudo-nitzschia, HSP genes were both up-and down-regulated [30].
The digestive gland of Pecten maximus contains principally two types of cells, secretory cells and digestive cells [16,40,76,77]. Vesicle-mediated transport and regulated exocytosis were enriched biological processes identified in protein networks obtained with the upregulated genes (Figure 2 and Supplementary File S8), therefore the activity of secretory cells in the digestive gland may be stimulated by domoic acid. These cells may be involved in the secretion of digestive enzymes [40,76,77].
Collagens are structural components of the extracellular matrix [78]. Components of collagen and proteins involved in cytoskeleton dynamics were among the proteins that appeared in the network obtained with the down-regulated genes in A. opercularis after exposure to domoic acid-producing Pseudo-nitzschia [30]. In the present work, several collagen genes were down-regulated (File S1) and cytoskeleton (Supplementary File S8) is one of the enriched cellular components in the down-regulated genes coding for proteins that showed interactions in the protein network analysis.

Conclusions
The domoic acid injected might have an effect on the gene expression in the digestive gland as reflected in the 535 DEGs found (226 up-regulated and 309 down-regulated). Some genes that code for putative glutamate receptors were expressed in the digestive gland of P. maximus, therefore part of the effects of domoic acid may be mediated by these receptors. A stress response at the level of gene expression, that could be caused by the domoic acid injection, was evidenced by the alteration of several biological, cellular, and molecular processes. Thus, protein networks obtained with the up-regulated genes were enriched in gene ontology (GO) terms, such as vesicle-mediated transport, response to stress, signal transduction, immune system process, RNA metabolic process, autophagy, lysosome, and oxidoreductase activity. On the other hand, networks obtained with the down-regulated genes were enriched in terms, such as response to stress, immune system process, ribosome biogenesis, signal transduction, mRNA processing, and oxidoreductase activity. In future research, it would be interesting investigate the domoic acid effects on gene expression in the kidneys and in the central nervous system (cerebral, pedal, and parietovisceral ganglia) of P. maximus.

Materials and Methods
The methods employed were similar as those previously described [29,30] except for minor modifications.

Animals
King scallops were obtained from the Ría de Arousa (Galicia, NW Spain) and maintained for a week in a 500 L tank, with a continuous unfiltered seawater flow (approximate) of 1200 L·h −1 . No domoic acid was detected in the routine monitoring of mollusks from the area and no toxic Pseudo-nitzschia cells were present in the area neither during the experiment nor during the previous month. All scallops, notwithstanding, contained domoic acid. It was impossible to obtain individuals free of toxin from the study area (and even from other European Countries) because the very low depuration rate of this species [14] makes easy their re-contamination. In the data obtained by the monitoring system run by Intecmar during the last 25 years in Galicia, it has been observed that the prevalence of domoic acid in the king scallop was 100%, and the same happened in other areas, such as Scotland [79]. The experimental approach was conditioned by this limitation, so we try to induce a response in the scallops by increasing the domoic acid levels that, in natural conditions, only undergo the progressive and slow changes derived from the depuration process. Fourteen scallops, with average height 10.98 ± 0.16 cm and 24.9 ± 1.3 g of soft tissues weight, were randomized into two groups: seven scallops as a control group and the other seven as a treated (with toxin injections) group. The scallops in the treated group were subjected to repeated injections of domoic acid (ABCAM) dissolved in filtered seawater and those in the control group, to equivalent injections of filtered seawater. One injection was made into the adductor muscle of each scallop every other day, for a period of 12 days. A volume of 62.5 µL of filtered seawater with a concentration of 8 µg domoic acid µL −1 was used for each injection (500 µg domoic acid in each injection). The total amount of domoic acid injected was 3000 µg per animal. The same volume of filtered seawater was injected into the scallops in the control group. When the scallops spontaneously opened the valves, a silicone stopper was placed between them to maintain them opened, and the injection was carried out.
After the injection, each scallop (treatment and control) was placed into a 5 L beaker, with 4 L of seawater water, with aeration, and maintained there for 24 h. After that period, the scallops were transferred to a 500 L tank with running seawater for one day. Following that, the injection process was repeated. Twenty-four hours after the last injection the scallops were dissected into digestive gland, gill, kidney, gonad, adductor muscle and remaining tissues. The dissected tissues were used for the determination of domoic acid content using liquid chromatography-tandem mass spectrometry (LC-MS/MS). A portion of the digestive gland was treated with RNAlater (ref. AM7021, Ambion, Life Technologies) and stored at −80 • C until the RNA extraction.

Determination of the Domoic Acid Content
Methanol for HPLC and formic acid were purchased from Labscan (Bangkok, Thailand) and Sigma Aldrich (St. Louis, MO, USA), respectively. Ultrapure water was obtained using a Milli-Q Gradient system, coupled to an Elix Advantage 10, both from Millipore (Merck Millipore, Darmstadt, Germany).
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 Ultraturax T25 (IKA, Staufen, Germany). The extract was clarified using centrifugation at 18,000× g at 4 • C for 10 min, retaining the 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 column (Phenomenex, Torrance, CA, USA). An elution gradient, with a flow of 280 µL min −1 , 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 one minute, linearly changed until reaching 55% B in minute 5, held for 2 min, and reverted to the initial conditions to equilibrate before the next injection. Five microliters 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 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
Twelve scallops (six obtained from the control group and six from the treated group) were subjected to RNA-seq analysis.

Library Preparation and Sequencing
Twelve cDNA libraries were generated. 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 oligo-attached 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 (Agilent Technologies, Santa Clara, CA, USA) high sensitivity assay; the quantity of the libraries was determined by real-time PCR in LightCycler 480 (Roche Diagnostics, Mannheim, Germany). Prior to cluster generation in cbot (Illumina), an equimolar pooling of the libraries was performed. The pool of the cDNA libraries was sequenced by paired-end sequencing (100 × 2 bp) on an Illumina HiSeq 2000 sequencer (Illumina, San Diego, CA, USA).

De Novo Assembly
Quality control checks of raw sequencing data were performed with FastQC. The technical adapters were eliminated using Trimgalore software version 0.3.3 (Trim Galore. Available Online: http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ (accessed on 10 March 2021)). Additionally, the reads with a mean Phred score > 30 were selected. Subsequently, the twelve samples were combined, and the complexity of the reads was reduced by removing duplicates. Then, a de novo assembly was performed using the programs Oasis, version 0.2.09 [80] and Trinity, version 2.1.1 [81]. 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. The completeness of the de novo assembly was evaluated with BUSCO [82] in OmicsBox software (BioBam Bioinformatics-2019, Valencia, Spain. OmicsBox-Version 1.4.11.) [83,84], using the Metazoa database (metazoa_odb10) [85], with a total of 954 orthologs (Blast E-value < 10 −3 ).
Each sample was then mapped with Bowtie2, version 2.2.6 [86] 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. Differential gene expression analysis of treated versus control samples was performed with DESeq2 algorithm version 1.8.2 (DESeq2. Available online: http://www.bioconductor.org/packages/devel/bioc/html/DESeq2.html (accessed on 10 March 2021)). The genes with a fold change of less than −1.5 or greater than 1.5, and a p-value adjusted using the Benjamini and Hochberg [87] method for controlling false discovery rate (FDR) < 0.05 were considered differentially expressed.

Protein Network Analysis
To search for the protein-protein interactions, network analyses using the String 10.5 algorithm [91] were performed. The putative human homologs of proteins coded by the up-and down-regulated genes in the P. maximus digestive gland were identified by means of a Blastx search [92] 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.
The genes that code for proteins that showed protein-protein interactions were subjected to GO enrichment analysis with OmicsBox using Fisher's exact test [93] (up-and down-regulated genes were analyzed separately). The false discovery rate (FDR) adjusted p-value [87] was set at a cutoff of 0.05. A normalization step using reference genes was performed for the relative expression of gene expression by means of RT-qPCR [94][95][96][97]. Only genes which show stable expression must be employed [39,98].
Four reference gene candidates (Table 7), GAPDH, EF1A, COX1, NDUFA7, and six target genes randomly selected (Table 7), MRP7, CYP2B4, P5CR, SLC6A9, FERRITIN, CYP2U1, were used in the gene expression study. The candidate reference genes have been successfully employed previously in P. maximus [39]. Oligonucleotide primers (Table 7) were synthesized by Integrated DNA Technologies. 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 (Real-time PCR Miner. Available online: http://www.miner.ewindup.info/ (accessed on 15 March 2021) [99]). The mean amplification efficiency (E) of each amplicon (Table 7)  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 plates using Bio-Rad iCycler iQ software V3.1 (IQ™ Real-Time PCR Detection System). The threshold value was set manually at 100 RFU to calculate the Cq values. Non-reverse transcriptase controls and non-template controls (NTC) were also included in each run.
Gene expression was normalized to reference genes that had stable expression levels [94][95][96][97]. The gene expression stability of candidate reference genes was analyzed using three Microsoft Excel based software applications, geNorm V3.5 [97], NormFinder V0.953 [94] and BestKeeper V1 [96]. 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) [39].

Statistical Analyses
The data were log-transformed to meet the requirements of normality and homogeneity of variances. The domoic acid concentration and domoic acid burden in control and treated scallops was compared using Student's t-test. The normalized expression of target genes (log2-transformed) in treated scallops, in relation to the control group, was also compared using Student's t-test. p < 0.05 was considered statistically significant. Statistical analyses were carried out with the IBM SPSS Statistics 24.0 package.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/toxins13050339/s1, Table S1: Summary of BUSCO analysis results obtained in the transcriptome of Pecten maximus digestive gland using the Metazoa database (metazoa_odb10), Figure  S1: MA plot showing log2 fold-change as a function of mean log expression level. The red dots represent genes with adjusted p-value < 0.05 and FC > 1.5 or <−1.5 (DEGs); the grey dots represent non-DEGs, Figure S2. Volcano plot. The red dots represent genes with adjusted p-value < 0.05 and FC > 1.5 or <−1.5 (DEGs); the grey dots represent non-DEGs, Figure S3: Species distribution of the top Blastx hits, File S1: List of differentially expressed genes. Sequence name, description, fold change (FC), FDR adjusted p-value (padj) and annotation results are shown, File S2: Nucleotide sequences of differentially expressed genes (in fasta format), File S3: List of transcripts coding for putative glutamate receptors in the digestive gland of Pecten maximus, File S4: Annotation of the transcripts expressed in the digestive gland of P. maximus. File S5: List of KO (KEGG Orthologs) for the differentially expressed genes and for all 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, File S8: Significantly enriched GO terms in the genes that code for proteins that showed interactions in the protein network analysis. The first two spreadsheets list the enriched GO terms for the up-regulated and the down-regulated genes respectively, the next two spreadsheets show the enriched GO terms after GO slim.