Next Article in Journal
Spatial and Temporal Variability in the Development and Potential Toxicity of Phormidium Biofilms in the Tarn River, France
Next Article in Special Issue
Impact of Dinophysis acuminata Feeding Mesodinium rubrum on Nutrient Dynamics and Bacterial Composition in a Microcosm
Previous Article in Journal
Combining E-Nose and Lateral Flow Immunoassays (LFIAs) for Rapid Occurrence/Co-Occurrence Aflatoxin and Fumonisin Detection in Maize
Previous Article in Special Issue
Dinophysis acuta in Scottish Coastal Waters and Its Influence on Diarrhetic Shellfish Toxin Profiles
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

RNA-Seq Analysis for Assessing the Early Response to DSP Toxins in Mytilus galloprovincialis Digestive Gland and Gill

Grupo Xenomar, Departamento de Bioloxía, Facultade de Ciencias and CICA (Centro de Investigacións Científicas Avanzadas), Universidade da Coruña, Campus de A Zapateira, 15071 A Coruña, Spain
*
Author to whom correspondence should be addressed.
Toxins 2018, 10(10), 417; https://doi.org/10.3390/toxins10100417
Submission received: 5 September 2018 / Revised: 11 October 2018 / Accepted: 13 October 2018 / Published: 16 October 2018
(This article belongs to the Special Issue Dinophysis Toxins: Distribution, Fate in Shellfish and Impacts)

Abstract

:
The harmful effects of diarrhetic shellfish poisoning (DSP) toxins on mammalian cell lines have been widely assessed. Studies in bivalves suggest that mussels display a resistance to the cytogenotoxic effects of DSP toxins. Further, it seems that the bigger the exposure, the more resistant mussels become. To elucidate the early genetic response of mussels against these toxins, the digestive gland and the gill transcriptomes of Mytilus galloprovincialis after Prorocentrum lima exposure (100,000 cells/L, 48 h) were de novo assembled based on the sequencing of 8 cDNA libraries obtained using an Illumina HiSeq 2000 platform. The assembly provided 95,702 contigs. A total of 2286 and 4523 differentially expressed transcripts were obtained in the digestive gland and the gill, respectively, indicating tissue-specific transcriptome responses. These transcripts were annotated and functionally enriched, showing 44 and 60 significant Pfam families in the digestive gland and the gill, respectively. Quantitative PCR (qPCR) was performed to validate the differential expression patterns of several genes related to lipid and carbohydrate metabolism, energy production, genome integrity and defense, suggesting their participation in the protective mechanism. This work provides knowledge of the early response against DSP toxins in the mussel M. galloprovincialis and useful information for further research on the molecular mechanisms of the bivalve resistance to these toxins.
Key Contribution: This work describes the transcriptome and gene expression profiles of M. galloprovincialis digestive gland and gill after early exposure to DSP toxins. Results showed that differentially expressed genes (DEGs) include genes involved in defense, immunity and metabolism, although some of them have been described as DEGs in response to other stimuli. This indicates that the mussel defense reaction is to some extent unspecific. This study also indicated that the expression of rpS4 and TPM genes in the digestive gland under these experimental conditions is stable and, therefore, these genes can be employed as reference genes to normalize gene expression in qPCR experiments carried out in mussels exposed to low concentrations of DSP toxins for short time periods.

Graphical Abstract

1. Introduction

Nowadays, harmful algal blooms (HABs) constitute one of the most important sources of natural contamination in the marine environment. This term refers not only to the phenomena originated by the proliferation of harmful algae, but also the phenomena caused by proliferation of toxic algae [1]. Although there is still a considerable absence of high quality time-series data in most regions affected by HABs [2], the blooms caused by the outbreaks of diarrhetic shellfish poisoning (DSP) toxin producing species seem to be associated with most of the HABs detected in European coasts [3]. These toxins are produced by dinoflagellates of the Dinophysis and Prorocentrum genera and constitute a heterogenous group of polyethers, including okadaic acid (OA) and its analogs, the dinophysis toxins (DTXs) [3,4,5,6,7,8]. In terms of abundance and consequent toxicity, OA is considered the main DSP toxin followed by DTX1, while DTX3—a less abundant DSP toxin—has become important because of its production through metabolic transformations that occur in some bivalves [7]. DTX1 seems to have similar toxicity levels to that of OA, while DTX2, DTX3 and DTX4 are less acutely toxic. On the other hand, the acylation of the 7-hydroxyl group with a saturated fatty acid forms compounds which are approximately 20 times less toxic than OA [9]. DSP toxins have a high lipophilic character, which allows for them to be accumulated in the fatty tissues of filter-feeding organisms—mainly in bivalve mollusks—and be transferred across the food chain, causing several gastrointestinal disorders [6]. Currently, efficient monitoring programs have been established by many countries to ban the harvesting of contaminated seafood and therefore, avert human intoxications [3]. However, seafood with small quantities of DSP toxins is still commercialized.
Since the ability of OA to inhibit several types of serine/threonine protein phosphatases was discovered by Bialojan and Takai [4], numerous works have studied the harmful effects of this toxic compound on different model systems, including different mammalian cell lines [8]. However, studies that assess the effects of these toxins in their main vectors—bivalve mollusks—are scarce. Recent studies carried out by our research group showed that DSP toxins cause more severe genotoxic and cytotoxic effects in bivalve cells at low concentrations and short exposition times, while these effects decrease or disappear as exposure increases in concentration and time [5,10,11,12]. This suggests that these organisms may have developed a quick protection mechanism against these toxic compounds. This may be associated with the accumulation, transformation and elimination of DSP toxins. This still unknown mechanism is of great interest for predicting the time course of toxic episodes and for reducing their negative consequences. With the aim of obtaining knowledge about this early genetic response, our research group has assessed the immediate effects caused by DSP toxins in the mussel Mytilus galloprovincialis using different stress indicators: DNA breaks, number of apoptotic cells [12], lipid peroxidation and antioxidant enzyme activities [10]. Although these indicators constitute a good approach to assess the first harmful effects produced by these toxins, they offer just a partial view on mussel response to toxic compounds. Taking this into account, it seems necessary to carry out analyses on the transcriptome response of mussels to DSP toxins to obtain a global perspective on their defense mechanisms against these toxins. Previous works used transcriptomic techniques to determine M. galloprovincialis transcriptome response to several stimuli, including marine toxins and pathogens [13,14,15,16,17,18,19]. Transcriptomic techniques such as RNA-Seq provide a valuable contribution to determining which gene pool expression is induced or suppressed depending on its physiological role in response to different treatments [20].
Some works have determined that the accumulation and distribution of DSP toxins in mussels is tissue specific [21,22]. The digestive gland is the mussel tissue that accumulates the most DSP toxins and is considered the main site of toxin bioconversion [23]. Furthermore, gills have numerous functions related to feeding, digestion and elimination of wastes and contaminants. The large surface and thin epithelium of the mussel gill make it an efficient site for direct interaction with the environment. Thus, gills efficiently capture suspended food particles—thanks to the mucus produced by them—and mediate their transport through the mussel mouth and digestive system [24].
In this work the whole transcriptome of the mussel M. galloprovincialis was de novo assembled and differentially expressed genes (DEGs) in digestive gland and gill after early exposure to DSP toxin-producer Prorocentrum lima were identified in order to determine the first response of these bivalve mollusks to these toxins and identify transcripts which could participate in the resistance mechanisms of mussels against the harmful effects of DSP toxins. Previous studies have characterized gene expression changes related to exposition to OA in bivalve mollusks [17,18,25,26] but to our knowledge, this is the first work that uses RNA-Seq to study the early transcriptional response of the mussel M. galloprovincialis to DSP toxins under short exposure to low concentrations of P. lima.

2. Results

2.1. Toxin Accumulation

According to the High Performance Liquid Chromatography/Mass Spectrometry (HPLC/MS) analyses, the P. lima strain AND-A0605 had an average toxin content of 0.4 pg OA/cell. Control mussels, fed with a mixture of Isochrysis galbana and Tetraselmis suecica, did not accumulate OA (<0.1 ng/g dry weight), while OA accumulated in treated mussels—fed also with P. lima—was 112.12 ng/g dry weight. Based on these results, and since these levels are well below the limit allowed by the European Commission Regulation for harvesting and sale (160 µg of OA equivalent/kg dry weight), we could consider that the mussels were exposed to low microalga cell densities, similar to those at the early stages of a HAB [27].

2.2. Transcriptome Sequencing and De Novo Assembly

In order to investigate the defense mechanisms of mussels exposed to DSP toxins, eight libraries derived from the digestive gland and the gill of the mussel M. galloprovincialis, in the absence of and under low densities of P. lima exposure, were constructed and sequenced using an Illumina sequencing platform. After de novo assembly with Trinity and Oases and their subsequent clustering by homology, 95,702 transcripts were obtained. Mean transcript size was 748 bp, with lengths ranging from as small as 100 bp to as a large as 16,082 bp. About 78% of the final assemblies were >200 bp and a N50 length of 1062 bp was obtained (Table 1).

2.3. DEGs Among Samples

Transcriptomic analyses were performed with the aim of identifying the main molecular mechanisms involved in the response of mussels to early contamination by DSP toxins. Using a RNA-Seq experiment, we generated transcriptome profiles for the digestive gland and the gill of the mussel M. galloprovincialis exposed to low densities of P. lima (100,000 cells/L) for a short period of time (48 h) and compared these data with profiles obtained from the digestive glands and the gills of control mussels. Sequences of all DEGs obtained are listed in File S1. A Venn diagram was used to depict the overlapping of DEGs when libraries were compared (Figure 1). Regarding the digestive gland, there were a total of 2286 DEGs between treatment and control groups, from which 1198 and 1088 transcripts were up- and down-regulated, respectively. Regarding the gills, there were a total of 4523 DEGs between both groups (treatment and control), from which 2579 and 1944 transcripts were up- and down-regulated, respectively. As a complementary analysis, the comparison of treated digestive glands and gills showed a total of 27,174 DEGs; 14,985 of them were up-regulated transcripts, while 12,189 were down-regulated (File S2). Only 26 transcripts out of all DEGs obtained were detected in all comparisons, with 17 and 9 of them being up- and down-regulated, respectively. The comparison of digestive glands and gills showed a total of 253 DEGs, from which 110 and 143 transcripts were up- and down-regulated, respectively. These DEGs could be useful for discovering genes involved in the early response to DSP toxins and, thereby, for identifying putative biomarkers for monitoring in advance of contamination episodes in the marine environment.

2.4. Gene Functional Annotations

Only 6% of the contigs included in the reference transcriptome showed BLAST similarity to proteins. About 20% of transcripts showed similarity to protein sequences deposited in the UniProt database and approximately 50% showed Pfam annotations. Thus, a relevant fraction of the contigs included in the reference transcriptome obtained in this work did not display any BLAST similarity or annotation.
Table 2, Table 3, Table 4 and Table 5 show the 25 most significantly up- and down-regulated genes in the digestive gland and the gill after exposure to low concentrations of DSP toxins (100,000 cells/L) for a short time period (48 h). Among the top over-represented DEGs in the digestive gland are genes that encode enzymes involved in the electron transport chain or mitochondrial oxidative phosphorilation (cytochrome c oxidase), as well as genes that encode ribosomal proteins or proteolytic enzymes (ribosomal protein L23a) (Table 2). Among the infra-represented genes in this tissue are also genes that encode enzymes of the electron transport chain (NADH dehydrogenase subunit 5) and ribosomal proteins (40S ribosomal protein S10-like). On the other hand, there are genes related to apoptosis (GTPase IMAP family member 7) and genes that encode proteins involved in the formation of nacre, promoting the crystallization of calcium carbonate (Perlucin) (Table 3). Similar to the digestive gland, among the over-represented genes in the gill (Table 4) are genes that encode enzymes of the electron transport chain (NADH dehydrogenase subunit 6) and proteins that play a role in the regulation of ion transport (calcyphosin-like protein). In contrast to the results obtained in the digestive gland, a gene encoding the cytochrome c oxidase subunit I is significantly down-regulated (Table 5). Also, a gene that encodes a protein involved in lipid metabolic processes and endocytosis is down-regulated in this tissue in the early response to DSP toxins (Table 5).
Functional enrichment studies performed using Pfam annotations obtained from the DEGs, showed 44 and 60 Pfam families significantly enriched in the digestive gland and the gill, respectively (File S3). Among these enriched domains, we found genes coding for proteins involved in GTP and calcium ion binding, transport, antibacterial activity and immune system in the digestive gland (Table 6). On the other hand, domains related to cell adhesion, cell-cell recognition, protein binding, immune system and correct folding of proteins were found in the gill (Table 7).
All DEGs from each tissue were classified according to the three main Gene Ontology (GO) aspects (biological processes, molecular functions and cellular components) and subcategories within (Figure 2 and Figure 3). Among the biological processes obtained for the digestive gland, proteolysis involved in the cellular protein catabolic process deserved special recognition for its down-regulation, while protein folding and translation are two of the most up-regulated processes. Regarding molecular functions, zinc and metal ion binding, as well as NADH dehydrogenase activity, showed considerable down-regulation in the digestive gland exposed to DSP toxins, while protein, GTP and RNA binding were up-regulated when the digestive gland responded to these toxins. The cellular components most involved in the response against DSP toxins seem to be the cytosol and the mitochondrion (cellular components up-regulated), while numerous sequences related to the extracellular exosome are down-regulated.
Regarding gills, the main down-regulated biological process when this tissue is exposed to DSP toxins is apoptosis. On the contrary, processes such as translational initiation or ATP synthesis coupled proton transport are over-represented after exposure to DSP toxins. When molecular functions are considered, RNA binding and NADH dehydrogenase activity are mostly up-regulated, while iron ion binding, sequence-specific DNA binding or cytochrome c oxidase activity are mainly down-regulated in the presence of DSP toxins. In this tissue, those cellular components most involved in the response against DSP toxins seem to be the nucleolus and the mitochondrion.

2.5. Real-Time Quantitative PCR (qPCR) Validation

We selected 10 DEGs for real-time qPCR confirmation based on their functions (lipid metabolism and immunity): seven up-regulated, two down-regulated and one with no differential expression. Regarding the digestive gland, big defensin 2 (BD2), NADH dehidrogenase subunit 5 (NADH5) and KAZAL domain containing protein (KAZAL DC) were up-regulated, GIY-YIG domain containing protein (GIY-YIG DC) was down-regulated and Dynactin-subunit-6-like (DYNA) showed no expression changes. Regarding the gills, Cytosolic phospholipase A-2 like (CPLA2), Arachidonate 15-lipoxygenase B-like (ALOX15B), Alpha-L-fucosidase-like (FUCA) and H_Lectin domain containing protein (H_Lectin DC) were up-regulated, while Fibrinogen_C domain containing protein (Fibrinogen_C DC) was down-regulated.
The heatmap provided in Figure 4 illustrates the expression levels of these genes in each library.
To confirm these patterns of expression by means of real-time qPCR, specific primers were designed. Sequences of these primers are shown in Table 8.
NormFinder software showed that rpS4 and TPM genes were the most stable genes and identified them as the best two-gene combination among all potential reference genes. These two genes also showed the lowest SD values when analyzed with BestKeeper. Moreover, their suitability as reference genes was supported by RefFinder results. Therefore, taking into account the combination of all results from the different analysis methods used (Table 9), TPM and rpS4 were identified as the most stable pair of reference genes in the digestive gland. These two genes were used for the normalization of gene expression in real-time qPCR.
The results of normalized expression (Figure 5 and Figure 6) validated the previous observations obtained using RNA-Seq. Fibrinogen_C DC, FUCA and NADH5 qPCR analyses were carried out using three biological replicates.

3. Discussion

Given the scarce knowledge of the resistance mechanisms involved in the early response of bivalve mollusks to marine toxins, the data presented in this work represent an important resource. Compared to other transcriptional works carried out in the digestive gland of the mussel M. galloprovincialis [13,15], a great number of DEGs were identified in the present study. This suggests a major impact of DSP toxins on gene expression regulation in the digestive gland and the gill of this species.
This study also revealed numerous transcripts assigned to Pfam families related to transport, cell adhesion, protein binding, calcium-binding proteins or immune system, among others. Many of these domains were also identified when haemolymph and digestive gland transcriptomes of mussel were analyzed in response to Vibrio alginnolyticus infection and domoic acid exposure [15,28,29]. Previous works carried out in bivalves exposed to marine toxins have shown significant changes in the expression levels of genes and proteins related to detoxification processes, such as cytochromes p450, ATP-binding cassette (ABC) transporters or glutathione S-transferases (GST) [12,15,30,31]. Surprisingly, although some of these genes are included among the DEGs in our results, they were not found among the most significant ones. Guo et al. [32] suggested the possible implication of p450 genes in OA metabolism in humans, generating new metabolites with less capacity to inhibit PP2A in comparison to OA. However, these transformations would not be completely effective to OA detoxification, which could explain our results.
Regarding GO, cellular organization and biogenesis, protein metabolism and modification, catabolism, response to stress and death and cell death are some of the biological processes most involved in the mussel response against toxins [33]. This is partially consistent with the main biological processes assigned in the present work when the digestive gland and the gill of mussels were exposed to DSP toxins. However, our data showed an important down-regulation of genes related to metabolic and apoptotic processes in the digestive gland and the gill, respectively, which may lay behind the first harmful effects of DSP toxins in these tissues. This result is not in agreement with the apoptosis induction observed in digestive glands when Mediterranean mussels were fed OA-contaminated nutrients [19]. Among the molecular functions involved in mussel response to toxins are protein binding, catalytic activity and transporter activity [33]. A similar result was obtained in the present work, although with important cytochrome-c oxidase and NADH dehydrogenase activities. On the other hand, the main cellular components shown in comparative transcriptomic studies of bivalves exposed to toxins were cytoplasm, nucleus, extracellular region and mitochondrion [33]. This is in agreement with some of the cellular components identified in the present work. However, our results also seem to show a key role of the extracellular exosome and respiratory chain in both mussel tissues—the digestive gland and the gill—in the early response to DSP toxins. Yamashita et al. [34] had already determined that exosomal secretion mechanisms are essential for methylmercury detoxification in the zebrafish embryo. Also, our work revealed an important participation of membrane integral components in the response to DSP toxins. This may be related to the known inhibitory effect of OA on intercellular channels in mammalian cells [35].
A large amount of contigs included in the reference transcriptome obtained in this study did not display any BLAST similarity or annotation, even with the recently sequenced M. galloprovincialis genome [36] or with Crassostrea gigas genome [37]. That was also the case for many of the top DEGs identified in this work that, despite their implication in the early response of mussel to DSP toxins, could not be identified. Similar results were obtained in a previous RNA-Seq study when digestive gland transcriptome of M. galloprovincialis was analyzed after exposure to the dinoflagellate Alexandrium minutum, a paralytic toxin producer [13]. Taking into account the length of some of these contigs as well as previous suggestions made by some authors, these sequences could be candidates to long non-coding RNA (lncRNA). lncRNA can regulate the activity of other genes by interacting with protein-coding mRNAs [38]. Milan et al. [39] observed that approximately 10% of the contigs obtained from the transcriptome of the clam Ruditapes philippinarum were originated by natural antisense transcription (NAT), a process that seems to be highly prevalent in bivalves.
When the data represented in the heatmap and the results obtained by qPCR were compared, a high correlation was observed between them, clear evidence that the RNA-Seq analysis conducted in this work was robust. Analyses in the digestive gland showed that the two most suitable genes for qPCR gene expression normalization were rpS4 and TPM. This result is in agreement with previous reports in which rpS4 was proposed as an optimal housekeeping gene to use under similar conditions [10,40]. To our knowledge, this is the first time that TPM is proposed and used in mussels to normalize qPCR data.
Our digestive gland data showed the up-regulation of different genes related to immune defense, including BD2, NADH5 and a KAZAL DC protein. Big defensins belong to a diverse family of peptides not only in terms of sequence, but also in terms of genomic organization and regulation of gene expression [41]. High gene expression levels of big defensin were identified in gills of the mussel Bathymodiolus azoricus [42]. Also, their up-regulation in haemocytes of the oyster C. gigas exposed to A. minutum—a paralytic toxin producer—has been associated with alterations in the immune system [43]. Similarly, big defensin gene expression was significantly up-regulated in haemolymph of the scallop Argopecten irradians when it was exposed to OA [44]. However, in a work in which the mussel M. galloprovincialis was exposed to the marine contaminant tris (1-choro-2-propyl) phosphate (TCPP) big defensin was down-regulated, affecting immunocompetence. Taking into account the high diversity of these genes [41] the big defensin identified in our work may correspond to a new variant of M. galloprovincialis related to the response to DSP toxins. On the other hand, the NADH protein family participates in transport and energy production. NADH is the third most frequently detected protein in comparative transcriptional studies that are carried out in bivalves exposed to different toxins [33]. Our results showed a significant increase in NADH5 gene expression. This is in line with an important up-regulation of NADH observed in a microarray designed based on data from normalized and suppression hybridization (SSH) libraries obtained from digestive gland and gill of the mussel M. galloprovincialis after exposition to sublethal concentrations of OA [17]. Our results also showed high expression levels of a putative KAZAL DC protein. Gerdol and Venier [45] have suggested that some bivalves can express Kazal-like protease inhibitors to counteract protease variants produced by invading microbes.
On the contrary, our digestive gland data showed the down-regulation of a putative GIY-YIG DC protein. This domain is present in many endonucleases involved in cellular processes such as DNA repair, the restriction of incoming foreing DNA, the movement of non-LTR retrotransposons or the maintenance of genome stability [46]. Indeed, Biscotti et al. [47] suggested that the expansion of this family in lungfish might be a genomic defense mechanism against the threat of spreading mobilome. Furthermore, Dittrich et al. [48] reported a gene which contains a GIY-YIG nuclease domain as an essential gene for proper DNA damage response in Caenorhabditis elegans embryos. However, mutants for this gene seem to have normal cell cycle arrest and apoptosis, which means this gene is not involved in the initial signalling process following DNA damage. This fact might partially explain the down-regulation of this transcript in the digestive gland of mussels during the early stages of DSP exposure, the situation simulated in the present study.
Our gill data showed an up-regulation of different genes related to lipid and carbohydrate metabolism, inflammatory response or immune defense, including CPLA2, ALOX15B, FUCA and a H_Lectin DC protein. CPLA2 is an enzyme that plays an important role as the primary generator of free arachidonic acid (AA)—a common precursor of a family of compounds with roles in inflammation [49]—released from membrane phospholipids. CPLA2 expression and activity are increased by reactive oxygen species (ROS) [50]. However, in a previous work, a decrease in lipid peroxidation levels was observed when mussel gills were exposed to the same DSP treatment [10]. This suggests the existence of an alternative defense mechanism. On the other hand, lipoxygenases (LOX) catalyze the generation of leukotrienes from AA producing byproducts that can function as ROS [51]. Some mussel extracts contain fatty acids with the ability to inhibit AA oxygenation by the cycloxigenase and LOX pathways, thus preventing inflammation [52]. In mammals, CPLA2 can cause membrane degradation, changes in plasma and mitochondrial membrane bioenergetics and permeability [53] and lysosomal membrane destabilization [54]. Indeed, CPLA2 is used as a stress indicator in biomonitoring programs. Some authors have also suggested that the up-regulation of genes involved in the inflammatory process, which was observed when digestive glands of the oyster C. gigas were exposed to P. lima, might represent a risk to this bivalve’s integrity [55]. Heavy metals functionally alter lysosomal membranes in haemocytes of mussels [56]. Ca2+ dependent CPLA2 enzymes play an important role in the lysosomal membrane destabilization induced by mercury and copper in the haemolymph cells of mussels [57]. Mussel gill exposed to low DSP toxin concentration produces an inflammatory response associated with the up-regulation of CPLA2 and ALOX15B that may be partially compensated by the up-regulation of antioxidant enzymes shown in many studies [10,58].
FUCA is an enzyme located in lysosomes and involved in carbohydrate metabolism. Based on our results, this gene seems to take part in the early response of mussel gills to DSP toxins. However, FUCA did not show gene deregulation when the gill of the scallop Nodipecten subnodosus was exposed to Gymnodinium catenatum, while an up-regulation was observed in the adductor muscle [59]. A down-regulation of FUCA protein was observed when the scallop Pecten maximus was exposed to hypoxia at different temperatures, suggesting an energy saving strategy by reducing protein turnover [60]. Nevertheless, the restriction of carbohydrate metabolism does not seem to be an important part of the early response of mussel gill to DSP toxins. Our gill data also showed up-regulation of a putative H_Lectin DC protein. This is a common finding in this type of studies, since type C lectins are usually overrepresented in bivalve transcriptomes exposed to marine toxins [17,61]. However, there is still relatively little information available about this domain related to cell adhesion and carbohydrate binding.
On the other hand, our gill data showed the down-regulation of a putative Fibrinogen_C DC protein. A study about the immune system of the mussel M. galloprovincialis identified fibrinogen as one of the most abundant transcripts in the Mytibase collection [62]. More specifically, C-terminal fibrinogen-like domain has a structure that binds to the carbohydrate residues of foreing and apoptotic cells. Indeed, some fibrinogen-like domains are included in many lectins [63] and, consequently, are involved in microorganism recognition by the activation of the lectin pathway, constituting a first line of immune defense. Although fibrinogen was first associated with haemolymph, the gill together with the digestive gland were the following tissues with the highest gene expression levels when three fibrinogen-related proteins were evaluated in the mussel M. galloprovincialis [64]. Down-regulation of fibrinogen was also observed when haemolymph of the scallop A. irradians was exposed to low concentrations of OA (50 nM) for short exposure times (48 h), suggesting the potential of this toxin to inhibit the ability of scallops to recognize and remove non-self particles [65]. Gene expression levels of Fibrinogen C also decreased when bay scallop gill tissue was exposed to 500 nM of OA for 48 h [30]. Differences in gene expression of fibrinogen C were also detected in the digestive gland of the mussel M. galloprovincialis after exposure to domoic acid-producing Pseudo-nitzschia [15]. However, fibrinogen gene expression was significantly up-regulated when the haemolymph of the scallop A. irradians was challenged with Listonella anguillarum [66] or when the haemolymph of the mussel Mytilus chilensis was exposed to saxitoxins [58]. It is important to note that, as in the case of big defensins, proteins that contain this domain present high individual variability. Thus, different mussels usually have different gene sequences, which demonstrates the extraordinary complexity of the immune system in these organisms [62].

4. Conclusions

This work represents the first RNA-Seq approach used in the mussel M. galloprovincialis to analyze tissue-specific mussel transcriptome after early exposure to DSP toxins. It describes the transcriptome and gene expression profiles of M. galloprovincialis digestive gland and gill, therefore increasing available genomic resources for this organism.
Furthermore, results showed that DEGs in early response to DSP toxins include genes involved in defense, immunity and metabolism, sheding some light into the resistance mechanisms that these organisms have against harmful effects of DSP toxins. In the digestive gland, BD2, KAZAL DC and NADH5 genes were up-regulated while GIY-YIG DC was down-regulated and DYNA showed no expression changes. On the other hand, ALOX15B, H_Lectin DC, CPLA2 and FUCA genes were up-regulated and Fibrinogen_C DC was down-regulated in gill. Nevertheless, many of the genes that responded to these toxins have been described as DEGs in response to other stimuli, indicating that the mussel defense reaction is to some extent unspecific, which may be beneficial when faced with other potentially harmful compounds.
This study also indicated that the expression of rpS4 and TPM genes in the digestive gland under these experimental conditions is stable and, therefore, these genes can be employed as reference genes to normalize gene expression in qPCR experiments carried out in mussels exposed to low concentrations of DSP toxins for short time periods.

5. Materials and Methods

5.1. Sample Collection and Experimental Design

Adult individuals of the mussel M. galloprovincialis (34 ± 0.5 mm anterior-posterior shell length) were collected from a natural population in the rocky shores of O Rañal beach (43°19′40.1″ N, 8°30′45.1″ W, A Coruña, NW Spain) in April 2015. This location (used by our research group in other studies [10]) was chosen as our sampling site based on its low density of DSP toxin-producing dinoflagellates [67]. The invertebrate animal experiment was assessed by the Spanish Ministry of Economy and Competitivity (project AGL2012-30897 approved on 28 December 2012). In the laboratory, specimens were acclimated for seven days at 17 °C with constant aeration in a photoperiod chamber with a 12 h light-dark cycle and fed twice a day with a 1:1 mixture of two cultures of nontoxic microalga species, I. galbana (3 × 106 cells/L) and T. suecica (12 × 106 cells/L). After acclimatization, mussels were randomly divided into two groups (n = 30 per experimental group) (Figure 7): a control group fed only with the microalga mixture used during acclimation period, and a treatment group additionally fed with 100,000 cells/L of the DSP toxin-producing alga P. lima. The culture of P. lima (strain AND-A0605) was obtained from the Quality Control Laboratory of Fishery Resources (Huelva, Spain). The treatment group was fed, four times a day, with 100,000 cells/L of P. lima during 48 h. These exposure characteristics were selected based on the results obtained in previous works by our research group in which these conditions showed the most interesting response at both the cytogenotoxic and the transcriptional level [10,12]. Cell concentrations of the nontoxic microalga cultures were determined by means of a Thoma cell counting chamber (Marienfeld, Lauda-Köningshofen, Germany), while that of the P. lima culture was estimated using the Sedgwich-Refter counting slide (Pyser-Sgi, Edenbridge, UK) after fixation with Lugol’s solution. After exposure, 12 individuals from each group—control and treatment—were dissected for digestive gland and gill tissues. These tissues were frozen in liquid nitrogen and stored at −80 °C until their use for RNA extraction, while the remaining individuals were used to estimate OA—the main DSP toxin—accumulation in the mussels by means of High Performance Liquid Chromatography/Mass Spectrometry (HPLC/MS). HPLC/MS analyses were carried out by the chromatography unit at Servizos de Apoio á Investigación (SAI)-University of A Coruña, following the protocol of the European Union Reference Laboratory for Marine Biotoxins [68].

5.2. RNA Extraction

Total RNA of digestive gland and gill from six control and six treated mussels was individually extracted using TRIzol (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s instructions (Figure 7). Isolated RNA was initially quantified using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). With the aim of reducing inter-individual variability, these RNAs were pooled (in equal quantities) in groups of three to provide a template for Illumina libraries (Figure 7). Additionally, quantity and integrity of RNA pools were checked using a Qubit 2.0 fluorometer (Life Technologies, Saint-Aubin, France) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), respectively.

5.3. Library Preparation and Sequencing

cDNA libraries were prepared and sequenced by Sistemas Genómicos (Valencia, Spain). Eight cDNA libraries were obtained from the digestive gland and the gill of mussels (two from control mussels and two from mussels exposed to P. lima, for each tissue, Figure 7). Poly(A)+mRNA fraction was isolated from total RNA and cDNA libraries were constructed following Illumina’s recommendations. cDNA libraries were sequenced using an Illumina HiSeq 2000 sequencer (Illumina, San Diego, CA, USA) and a paired-end sequencing strategy (100 × 2 bp). Raw data are accessible from the NCBI Short Read Archive (SRA accession: SRP158485).

5.4. De Novo Assembly

A preliminary bioinformatic analysis was performed by Sistemas Genómicos (Valencia, Spain). Initially, short sequence reads were quality checked using FastQC [69] and the TrueSeq adapters were trimmed using Trim Galore software version 0.3.3 (Babraham Bioinformatics, Cambridge, UK), keeping those reads with a mean phred score >30. With the aim of obtaining a reference transcriptome, all generated results were combined in a single data set. Then, low quality reads were re-identified and removed using PrinSeq-lite software version 0.20.4 [70], while duplicate reads were then removed using FastX-Toolkit (fastx_collapser option) [71]. Subsequently, de novo transcriptome assembly was conducted with the software Oases (version 2.0.9) and Trinity (version 2.1.1). Both assemblies were correlated by combining contigs with sequence similarity (>90% homology) using cd-hit (version 4.6). Potential ORFs were predicted using TransDecoder (version 2.0) with default settings. Then, each library was mapped against the reference transcriptome obtained in the previous step using Bowtie2 (version 2.2.6) and high quality reads were selected—high mapping quality with a 1 × 10−4 error probability—to increase count expression resolution. Finally, expression inference was carried out using the counts of properly paired reads by transcript.

5.5. Differential Expression, Functional Annotation and Functional Enrichment Analysis of DEGs

The expression of each sample was normalized by library size (initial number of reads) using the R package DESeq2 version 1.8.2 [72] (R software version 3.2.3 [73]) based on a negative binomial distribution, with the aim of analyzing differential expression. Those genes with a fold change lower than −2 or higher than 2, and an adjusted p-value < 0.05 were considered differentially expressed. Additionally, the method for controlling FDR was used to calculate the adjusted p-values [74].
DEGs were initially annotated using blastx against UniProt database and blastn against the NCBI nucleotide database, using an E-value threshold of 0.01. Subsequently, sequences annotated with RNAs were identified, while sequences associated with P. lima were removed from further analysis. Additionally, DEGs were re-annotated by a blastx analysis (ncbi-blast/2.3.0+)—using an e-value of 1 × 10−6 as cut-off—performed through the Supercomputing Centre of Galicia (CESGA). Subsequently, to know the biological processes, molecular functions and cellular components related to DEGs, annotated sequences were analyzed using GO implemented in Blast2GO software [75,76]. A functional enrichment analysis was performed using the Pfam [77] functional information, with the aim of annotating protein domains. Additionally, a subset of annotated DEGs was selected based on their biological function and their gene expression levels were represented in a heat map using CIMminer [78].

5.6. Real-Time Quantitative PCR Validation

A subset of annotated DEGs was selected based on their biological function to validate their gene expression using real-time qPCR. Reference genes for expression quantification were selected among six potential candidate housekeeping genes, including two primers for 18S ribosomal RNA (18S) [79], ribosomal protein S4 (rpS4), glyceraldehyde 3-phosphate-dehydrogenase (GAPDH) [40], elongation factor 1 (EF1) [10] and tropomyosin (TPM). TPM primers were designed as part of this work from an annotated gene with very stable expression levels. These primers and the specific primers to amplify the selected DEGs were designed using the Universal Probe Library software [80] (Roche Diagnostics, Mannheim, Germany). Primer specificities were verified using agarose gel electrophoresis, showing one single DNA product of the expected length. Two different algorithms, Normfinder and BestKeeper, were initially used to rank candidate reference genes according to their stability in the digestive gland and to decide on the optimal number of reference genes required for accurate normalization. Normfinder was used with R version 3.0.1 [73] and BestKeeper is an Excel-based tool that uses pairwise correlations [81]. Whenever BestKeeper analysis showed genes with SD values > 1, those genes were excluded from correlation coefficient calculations. Subsequently, results were checked using RefFinder [82], a web-based tool that integrates four different algorithms (Normfinder, BesKeeper, GeNorm and Delta Ct).
RNA samples from those individuals previously used for library preparation were used for the real-time qPCR validation. Four independent biological replicates and two technical replicates were analyzed together using the sample maximization approach [83]. cDNA was synthesized using 1 µg of RNA using the First Strand cDNA Synthesis kit according to the manufacturer’s instructions (Roche Diagnostics, Mannheim, Germany). qPCR amplifications were carried out using the FastStart Essential DNA Green Master kit (Roche Diagnostics, Mannheim, Germany) following the manufacturer’s instructions with the following modifications. All reactions were performed in a final volume of 20 µL of master mix containing 6.4 µL H2O, 0.8 µL of each primer (10 µM), 10 µL of the SYBR Green Mix (Roche Diagnostics, Mannheim, Germany) and 2 µL of each reverse transcribed RNA (cDNA). Reactions consisted of an initial denaturation step of 10 min at 95 °C, followed by an amplification of the target cDNA for 40 cycles (denaturation at 95 °C for 10 s, annealing at 60 °C for 10 s, elongation at 72 °C for 10 s), melting curve analysis (1 cycle at 95 °C for 5 s, 65 °C for 60 s and 95 °C for 1 s), and cooling at 40 °C for 20 s. Specificity of the qPCR product was analyzed by melting curve analysis.
Efficiency of the reaction for each mRNA was determined using LinRegPCR 2014.x software [84]. Gene relative expression levels were normalized using rpS4 and TPM as reference genes. For data analyses, Cq values were extracted with the qPCR instrument software LightCycler Software 1.5.0 (Roche Diagnostics, Mannheim, Germany). Cq values were then exported to Excel (Microsoft, Redmond, WA, USA), and differences in expression were calculated using the Pfaffl method with two reference genes [85]. Whenever a single individual sample showed a Cq value with an over five point difference to the mean Cq for the condition, that value was considered an amplification error, therefore, that sample was removed and analyses were carried out using three biological replicates instead of four. Normalized relative quantities (NRQ) for each gene were represented in bar plots (control vs. treatment) using GraphPad Prism version 6 (GraphPad Prism Software Inc., La Jolla, CA, USA). For better visualization of results some data were log transformed for graphic representation. Differences in gene expression between control and treatment samples were determined by Mann-Whitney non-parametric U test using the SPSS IBM software package version 22 (IBM, Armon, NY, USA). An additional analysis to confirm the obtained gene expression differences was conducted in REST 2009 (Qiagen, Hilden, Germany) [86].

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2072-6651/10/10/417/s1, File S1: Nucleotide sequences—in fasta format—of all differentially expressed genes (DEGs), File S2: List of DEGs. Each spreadsheet shows DEGs from the comparison of either the same tissue under different conditions—MBT_vs_MBC_DEGs for gills and MGT_vs_MGC_DEGs for digestive gland—or two tissues under the same condition—MGT_vs_MBT_DEGs for treated digestive gland and gill. For each DEG, sequence ID, baseMean, length, Log2 Fold Change (FC), FC, p-value and adjusted p-value are given. Also, Blast_nucleotide, Blast_UniProt and Pfam columns show the best hit against Nucleotide, UniProt and Pfam databases, respectively, File S3: List of Pfam families functionally enriched. Each spreadsheet shows DEGs from the comparison of the same tissue under different conditions—MGT_vs_MGC for digestive gland and MBT_vs_MBC for gill. For each Pfam category, number of genes, p-value, expression patterns and gene IDs are given.

Author Contributions

M.V.P.-F. and J.M. conceived and designed the experiments; M.V.P.-F. and L.M. performed the experiments and analyzed the data; J.M. contributed reagents/materials/analysis tools; M.V.P.-F., L.M. and J.M. wrote and revised the manuscript.

Funding

This study was supported by grants from the Spanish Ministry of Economy and Competitivity (AGL2012-30897, Josefina Mendez) and M.V.P.-F. work was funded through a fellowship by Deputación da Coruña (BINV-CC/2017).

Acknowledgments

The authors would like to thank the Consello Regulador do Mexillón de Galicia for its support, Juan C. Triviño for help with bioinformatic analyses and CESGA (www.cesga.es) in Santiago de Compostela, Spain for access to computing facilities.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Anderson, D.M.; Andersen, P.; Bricelj, V.M.; Cullen, J.J.; Rensel, J.E.J. Monitoring and Management Strategies for Harmful Algal Blooms in Coastal Waters; Unesco: Paris, France, 2001. [Google Scholar]
  2. Wells, M.L.; Trainer, V.L.; Smayda, T.J.; Karlson, B.S.O.; Trick, C.G.; Kudela, R.M.; Ishikawa, A.; Bernard, S.; Wulff, A.; Anderson, D.M.; et al. Harmful algal blooms and climate change: Learning from the past and present to forecast the future. Harmful Algae 2015, 49, 68–93. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Visciano, P.; Schirone, M.; Berti, M.; Milandri, A.; Tofalo, R.; Suzzi, G. Marine biotoxins: Occurrence, toxicity, regulatory limits and reference methods. Front. Microbiol. 2016, 7, 1051. [Google Scholar] [CrossRef] [PubMed]
  4. Bialojan, C.; Takai, A. Inhibitory effect of a marine-sponge toxin, okadaic acid, on protein phosphatases. Specificity and kinetics. Biochem. J. 1988, 256, 283–290. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Prado-Alvarez, M.; Flórez-Barrós, F.; Sexto-Iglesias, A.; Méndez, J.; Fernandez-Tajes, J. Effects of okadaic acid on haemocytes from Mytilus galloprovincialis: A comparison between field and laboratory studies. Mar. Environ. Res. 2012, 81, 90–93. [Google Scholar] [CrossRef] [PubMed]
  6. Prego-Faraldo, M.V.; Valdiglesias, V.; Méndez, J.; Eirín-López, J.M. Okadaic acid meet and greet: An insight into detection methods, response strategies and genotoxic effects in marine invertebrates. Mar. Drugs 2013, 11, 2829–2845. [Google Scholar] [CrossRef] [PubMed]
  7. Reguera, B.; Velo-Suárez, L.; Raine, R.; Park, M.G. Harmful dinophysis species: A review. Harmful Algae 2012, 14, 87–106. [Google Scholar] [CrossRef]
  8. Valdiglesias, V.; Prego-Faraldo, M.V.; Pásaro, E.; Méndez, J.; Laffon, B. Okadaic acid: More than a diarrheic toxin. Mar. Drugs 2013, 11, 4328–4349. [Google Scholar] [CrossRef] [PubMed]
  9. Munday, R. Is protein phosphatase inhibition responsible for the toxic effects of okadaic acid in animals? Toxins 2013, 5, 267–285. [Google Scholar] [CrossRef] [PubMed]
  10. Prego-Faraldo, M.; Vieira, L.; Eirin-Lopez, J.; Méndez, J.; Guilhermino, L. Transcriptional and biochemical analysis of antioxidant enzymes in the mussel Mytilus galloprovincialis during experimental exposures to the toxic dinoflagellate Prorocentrum lima. Mar. Environ. Res. 2017, 129, 304–315. [Google Scholar] [CrossRef] [PubMed]
  11. Prego-Faraldo, M.V.; Valdiglesias, V.; Laffon, B.; Eirín-López, J.M.; Méndez, J. In vitro analysis of early genotoxic and cytotoxic effects of okadaic acid in different cell types of the mussel Mytilus galloprovincialis. J. Toxicol. Environ. Health A 2015, 78, 814–824. [Google Scholar] [CrossRef] [PubMed]
  12. Prego-Faraldo, M.V.; Valdiglesias, V.; Laffon, B.; Mendez, J.; Eirin-Lopez, J.M. Early genotoxic and cytotoxic effects of the toxic dinoflagellate Prorocentrum lima in the mussel Mytilus galloprovincialis. Toxins 2016, 8, 159. [Google Scholar] [CrossRef] [PubMed]
  13. Gerdol, M.; De Moro, G.; Manfrin, C.; Milandri, A.; Riccardi, E.; Beran, A.; Venier, P.; Pallavicini, A. RNA sequencing and de novo assembly of the digestive gland transcriptome in Mytilus galloprovincialis fed with toxinogenic and non-toxic strains of Alexandrium minutum. BMC Res. Notes 2014, 7, 722. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Moreira, R.; Pereiro, P.; Canchaya, C.; Posada, D.; Figueras, A.; Novoa, B. RNA-Seq in Mytilus galloprovincialis: Comparative transcriptomics and expression profiles among different tissues. BMC Genom. 2015, 16, 728. [Google Scholar] [CrossRef] [PubMed]
  15. Pazos, A.J.; Ventoso, P.; Martínez-Escauriaza, R.; Pérez-Parallé, M.L.; Blanco, J.; Triviño, J.C.; Sánchez, J.L. Transcriptional response after exposure to domoic acid-producing Pseudo-nitzschia in the digestive gland of the mussel Mytilus galloprovincialis. Toxicon 2017, 140, 60–71. [Google Scholar] [CrossRef] [PubMed]
  16. Rosani, U.; Varotto, L.; Rossi, A.; Roch, P.; Novoa, B.; Figueras, A.; Pallavicini, A.; Venier, P. Massively parallel amplicon sequencing reveals isotype-specific variability of antimicrobial peptide transcripts in Mytilus galloprovincialis. PLoS ONE 2011, 6, e26680. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Suarez-Ulloa, V.; Fernandez-Tajes, J.; Aguiar-Pulido, V.; Prego-Faraldo, M.V.; Florez-Barros, F.; Sexto-Iglesias, A.; Mendez, J.; Eirin-Lopez, J.M. Unbiased high-throughput characterization of mussel transcriptomic responses to sublethal concentrations of the biotoxin okadaic acid. PeerJ 2015, 3, e1429. [Google Scholar] [CrossRef] [PubMed]
  18. Suárez-Ulloa, V.; Fernández-Tajes, J.; Aguiar-Pulido, V.; Rivera-Casas, C.; González-Romero, R.; Ausio, J.; Méndez, J.; Dorado, J.; Eirín-López, J.M. The CHROMEVALOA database: A resource for the evaluation of okadaic acid contamination in the marine environment based on the chromatin-associated transcriptome of the mussel Mytilus galloprovincialis. Mar. Drugs 2013, 11, 830–841. [Google Scholar] [CrossRef] [PubMed]
  19. Manfrin, C.; Dreos, R.; Battistella, S.; Beran, A.; Gerdol, M.; Varotto, L.; Lanfranchi, G.; Venier, P.; Pallavicini, A. Mediterranean mussel gene expression profile induced by okadaic acid exposure. Environ. Sci. Technol. 2010, 44, 8276–8283. [Google Scholar] [CrossRef] [PubMed]
  20. Wang, Z.; Gerstein, M.; Snyder, M. RNA-Seq: A revolutionary tool for transcriptomics. Nat. Rev. Genet. 2009, 10, 57. [Google Scholar] [CrossRef] [PubMed]
  21. Blanco, J.; Mariño, C.; Martín, H.; Acosta, C.P. Anatomical distribution of diarrhetic shellfish poisoning (DSP) toxins in the mussel Mytilus galloprovincialis. Toxicon 2007, 50, 1011–1018. [Google Scholar] [CrossRef] [PubMed]
  22. Moroño, A.; Arévalo, F.; Fernández, M.; Maneiro, J.; Pazos, Y.; Salgado, C.; Blanco, J. Accumulation and transformation of DSP toxins in mussels Mytilus galloprovincialis during a toxic episode caused by Dinophysis acuminata. Aquat. Toxicol. 2003, 62, 269–280. [Google Scholar] [CrossRef]
  23. Manfrin, C.; De Moro, G.; Torboli, V.; Venier, P.; Pallavicini, A.; Gerdol, M. Physiological and molecular responses of bivalves to toxic dinoflagellates. Invertebr. Surv. J. 2012, 9, 184–199. [Google Scholar]
  24. Beyer, J.; Green, N.W.; Brooks, S.; Allan, I.J.; Ruus, A.; Gomes, T.; Bråte, I.L.N.; Schøyen, M. Blue mussels (Mytilus edulis spp.) as sentinel organisms in coastal pollution monitoring: A review. Mar. Environ. Res. 2017. [Google Scholar] [CrossRef] [PubMed]
  25. Romero-Geraldo, R.d.J.; García-Lagunas, N.; Hernandez-Saavedra, N.Y. Effects of in vitro exposure to diarrheic toxin producer Prorocentrum lima on gene expressions related to cell cycle regulation and immune response in Crassostrea gigas. PLoS ONE 2014, 9, e97181. [Google Scholar] [CrossRef]
  26. Romero-Geraldo, R.d.J.; Hernández-Saavedra, N.Y. Stress gene expression in Crassostrea gigas (Thunberg, 1793) in response to experimental exposure to the toxic dinoflagellate Prorocentrum lima (Ehrenberg) Dodge, 1975. Aquac. Res. 2014, 45, 1512–1522. [Google Scholar] [CrossRef]
  27. Díaz, P.A.; Reguera, B.; Ruiz-Villarreal, M.; Pazos, Y.; Velo-Suárez, L.; Berger, H.; Sourisseau, M. Climate variability and oceanographic settings associated with interannual variability in the initiation of Dinophysis acuminata blooms. Mar. Drugs 2013, 11, 2964–2981. [Google Scholar] [CrossRef] [PubMed]
  28. Dong, W.; Chen, Y.; Lu, W.; Wu, B.; Qi, P. Transcriptome analysis of Mytilus coruscus hemocytes in response to Vibrio alginnolyficus infection. Fish Shellfish Immunol. 2017, 70, 560–567. [Google Scholar] [CrossRef] [PubMed]
  29. Ventoso, P.; Martínez-Escauriaza, R.; Sánchez, J.; Pérez-Parallé, M.; Blanco, J.; Triviño, J.; Pazos, A. In Sequencing and de novo assembly of the digestive gland transcriptome in Mytilus galloprovincialis and analysis of differentially expressed genes in response to domoic acid. In Proceedings of the International Symposium on Genetics in Aquaculture XII, Santiago de Compostela, Spain, 21–27 June 2015; p. 93, 229. [Google Scholar]
  30. Chi, C.; Giri, S.; Jun, J.; Kim, S.; Kim, H.; Kang, J.; Park, S. Detoxification- and immune-related transcriptomic analysis of gills from bay scallops (Argopecten irradians) in response to algal toxin okadaic acid. Toxins 2018, 10, 308. [Google Scholar] [CrossRef] [PubMed]
  31. Huang, L.; Zou, Y.; Weng, H.-W.; Li, H.-Y.; Liu, J.-S.; Yang, W.-D. Proteomic profile in Perna viridis after exposed to Prorocentrum lima, a dinoflagellate producing DSP toxins. Environ. Pollut. 2015, 196, 350–357. [Google Scholar] [CrossRef] [PubMed]
  32. Guo, F.; An, T.; Rein, K.S. The algal hepatoxoxin okadaic acid is a substrate for human cytochromes CYP3A4 and CYP3A5. Toxicon 2010, 55, 325–332. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Miao, J.; Chi, L.; Pan, L.; Song, Y. Generally detected genes in comparative transcriptomics in bivalves: Toward the identification of molecular markers of cellular stress response. Environ. Toxicol. Pharmacol. 2015, 39, 475–481. [Google Scholar] [CrossRef] [PubMed]
  34. Yamashita, M.; Yamashita, Y.; Suzuki, T.; Kani, Y.; Mizusawa, N.; Imamura, S.; Takemoto, K.; Hara, T.; Hossain, M.A.; Yabu, T. Selenoneine, a novel selenium-containing compound, mediates detoxification mechanisms against methylmercury accumulation and toxicity in zebrafish embryo. Mar. Biotechnol. 2013, 15, 559–570. [Google Scholar] [CrossRef] [PubMed]
  35. Creppy, E.E.; Traoré, A.; Baudrimont, I.; Cascante, M.; Carratú, M.-R. Recent advances in the study of epigenetic effects induced by the phycotoxin okadaic acid. Toxicology 2002, 181, 433–439. [Google Scholar] [CrossRef]
  36. Murgarella, M.; Puiu, D.; Novoa, B.; Figueras, A.; Posada, D.; Canchaya, C. A first insight into the genome of the filter-feeder mussel Mytilus galloprovincialis. PLoS ONE 2016, 11, e0151561. [Google Scholar] [CrossRef]
  37. Zhang, G.; Fang, X.; Guo, X.; Li, L.; Luo, R.; Xu, F.; Yang, P.; Zhang, L.; Wang, X.; Qi, H. The oyster genome reveals stress adaptation and complexity of shell formation. Nature 2012, 490, 49–54. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Ilott, N.E.; Ponting, C.P. Predicting long non-coding RNAs using RNA sequencing. Methods 2013, 63, 50–59. [Google Scholar] [CrossRef] [PubMed]
  39. Milan, M.; Coppe, A.; Reinhardt, R.; Cancela, L.M.; Leite, R.B.; Saavedra, C.; Ciofi, C.; Chelazzi, G.; Patarnello, T.; Bortoluzzi, S. Transcriptome sequencing and microarray development for the manila clam, Ruditapes philippinarum: Genomic tools for environmental monitoring. BMC Genom. 2011, 12, 234. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. Lozano, V.; Martínez-Escauriaza, R.; Pérez-Parallé, M.; Pazos, A.; Sánchez, J. Two novel multidrug resistance associated protein (MRP/ABCC) from the mediterranean mussel (Mytilus galloprovincialis): Characterization and expression patterns in detoxifying tissues. Can. J. Zool. 2015, 93, 567–578. [Google Scholar] [CrossRef]
  41. Rosa, R.D.; Santini, A.; Fievet, J.; Bulet, P.; Destoumieux-Garzón, D.; Bachère, E. Big defensins, a diverse family of antimicrobial peptides that follows different patterns of expression in hemocytes of the oyster Crassostrea gigas. PLoS ONE 2011, 6, e25594. [Google Scholar] [CrossRef] [PubMed]
  42. Bettencourt, R.; Pinheiro, M.; Egas, C.; Gomes, P.; Afonso, M.; Shank, T.; Santos, R.S. High-throughput sequencing and analysis of the gill tissue transcriptome from the deep-sea hydrothermal vent mussel Bathymodiolus azoricus. BMC Genom. 2010, 11, 559. [Google Scholar] [CrossRef] [PubMed]
  43. Mello, D.F.; da Silva, P.M.; Barracco, M.A.; Soudant, P.; Hégaret, H. Effects of the dinoflagellate Alexandrium minutum and its toxin (saxitoxin) on the functional activity and gene expression of Crassostrea gigas hemocytes. Harmful Algae 2013, 26, 45–51. [Google Scholar] [CrossRef]
  44. Chi, C.; Giri, S.S.; Jun, J.W.; Kim, H.J.; Kim, S.W.; Yun, S.; Park, S.C. Effects of algal toxin okadaic acid on the non-specific immune and antioxidant response of bay scallop (Argopecten irradians). Fish Shellfish Immunol. 2017, 65, 111–117. [Google Scholar] [CrossRef] [PubMed]
  45. Gerdol, M.; Venier, P. An updated molecular basis for mussel immunity. Fish Shellfish Immunol. 2015, 46, 17–38. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Dunin-Horkawicz, S.; Feder, M.; Bujnicki, J.M. Phylogenomic analysis of the GIY−YIG nuclease superfamily. BMC Genom. 2006, 7, 98. [Google Scholar] [CrossRef] [PubMed]
  47. Biscotti, M.A.; Gerdol, M.; Canapa, A.; Forconi, M.; Olmo, E.; Pallavicini, A.; Barucca, M.; Schartl, M. The lungfish transcriptome: A glimpse into molecular evolution events at the transition from water to land. Sci. Rep. 2016, 6, 21571. [Google Scholar] [CrossRef] [PubMed]
  48. Dittrich, C.M.; Kratz, K.; Sendoel, A.; Gruenbaum, Y.; Jiricny, J.; Hengartner, M.O. LEM−3–A LEM domain containing nuclease involved in the DNA damage response in C. elegans. PLoS ONE 2012, 7, e24555. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Balsinde, J. Phospholipase A2; Cellular Regulation, Function, and Inhibition 2016. Available online: http://www.balsinde.org/publists/engplasic.pdf (accessed on 24 May 2018).
  50. Korbecki, J.; Baranowska-Bosiacka, I.; Gutowska, I.; Chlubek, D. The effect of reactive oxygen species on the synthesis of prostanoids from arachidonic acid. J. Physiol. Pharmacol. 2013, 64, 409–421. [Google Scholar] [PubMed]
  51. Kim, C.; Kim, J.-Y.; Kim, J.-H. Cytosolic phospholipase A2, lipoxygenase metabolites, and reactive oxygen species. BMB Rep. 2008, 41, 555–559. [Google Scholar] [CrossRef] [PubMed]
  52. Bierer, T.L.; Bui, L.M. Improvement of arthritic signs in dogs fed green-lipped mussel (Perna canaliculus). J. Nutr. 2002, 132, 1634S–1636S. [Google Scholar] [CrossRef] [PubMed]
  53. Zhao, M.; Brunk, U.T.; Eaton, J.W. Delayed oxidant-induced cell death involves activation of phospholipase A2. FEBS Lett. 2001, 509, 399–404. [Google Scholar] [CrossRef] [Green Version]
  54. Mukherjee, A.; Ghosal, S.; Maity, C. Lysosomal membrane stabilization by α-tocopherol against the damaging action of Vipera russelli venom phospholipase A2. Cell. Mol. Life Sci. 1997, 53, 152–155. [Google Scholar] [CrossRef] [PubMed]
  55. Romero-Geraldo, R.d.J.; García-Lagunas, N.; Hernández-Saavedra, N.Y. Crassostrea gigas exposure to the dinoflagellate Prorocentrum lima: Histological and gene expression effects on the digestive gland. Mar. Environ. Res. 2016, 120, 93–102. [Google Scholar] [CrossRef] [PubMed]
  56. Viarengo, A.; Marro, A.; Marchi, B.; Burlando, B. Single and combined effects of heavy metals and hormones on lysosomes of haemolymph cells from the mussel Mytilus galloprovincialis. Mar. Biol. 2000, 137, 907–912. [Google Scholar] [CrossRef]
  57. Marchi, B.; Burlando, B.; Moore, M.; Viarengo, A. Mercury- and copper-induced lysosomal membrane destabilisation depends on [Ca2+]i dependent phospholipase A2 activation. Aquat. Toxicol. 2004, 66, 197–204. [Google Scholar] [CrossRef] [PubMed]
  58. Núñez-Acuña, G.; Aballay, A.E.; Hégaret, H.; Astuya, A.P.; Gallardo-Escárate, C. Transcriptional responses of Mytilus chilensis exposed in vivo to saxitoxin (STX). J. Mollus. Stud. 2013, 79, 323–331. [Google Scholar] [CrossRef]
  59. Estrada, N.; de Jesús Romero, M.; Campa-Córdova, A.; Luna, A.; Ascencio, F. Effects of the toxic dinoflagellate, Gymnodinium catenatum on hydrolytic and antioxidant enzymes, in tissues of the giant lions-paw scallop Nodipecten subnodosus. Comp. Biochem. Phys. C Toxicol. Pharmacol. 2007, 146, 502–510. [Google Scholar] [CrossRef] [PubMed]
  60. Artigaud, S.; Lacroix, C.; Richard, J.; Flye-Sainte-Marie, J.; Bargelloni, L.; Pichereau, V. Proteomic responses to hypoxia at different temperatures in the great scallop (Pecten maximus). PeerJ 2015, 3, e871. [Google Scholar] [CrossRef] [PubMed]
  61. Detree, C.; Núñez-Acuña, G.; Roberts, S.; Gallardo-Escárate, C. Uncovering the complex transcriptome response of Mytilus chilensis against saxitoxin: Implications of harmful algal blooms on mussel populations. PLoS ONE 2016, 11, e0165231. [Google Scholar] [CrossRef] [PubMed]
  62. Venier, P.; Varotto, L.; Rosani, U.; Millino, C.; Celegato, B.; Bernante, F.; Lanfranchi, G.; Novoa, B.; Roch, P.; Figueras, A. Insights into the innate immunity of the mediterranean mussel Mytilus galloprovincialis. BMC Genom. 2011, 12, 69. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Domeneghetti, S.; Manfrin, C.; Varotto, L.; Rosani, U.; Gerdol, M.; De Moro, G.; Pallavicini, A.; Venier, P. How gene expression profiles disclose vital processes and immune responses in Mytilus spp. Invertebr. Surv. J. 2011, 8, 179–189. [Google Scholar]
  64. Romero, A.; Dios, S.; Poisa-Beiro, L.; Costa, M.M.; Posada, D.; Figueras, A.; Novoa, B. Individual sequence variability and functional activities of fibrinogen-related proteins (FREPs) in the mediterranean mussel (Mytilus galloprovincialis) suggest ancient and complex immune recognition models in invertebrates. Dev. Comp. Immunol. 2011, 35, 334–344. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Cheng, C. Physico-Immunological Characterizations of Exogenous Substances (Palmitoleic Acid and Okadaic Acid) in Bivalves. Ph.D. Thesis, The Graduate School of Seoul National University, Seoul, Korea, 2017. [Google Scholar]
  66. Zhang, X.-J.; Qin, G.-M.; Yan, B.-L.; Xu, J.; Bi, K.-R.; Qin, L. Phenotypic and molecular characterization of pathogenic Listonella anguillarum isolated from half-smooth tongue sole Cynoglossus semilaevis. Acta Oceanol. Sin. 2009, 5, 012. [Google Scholar]
  67. Intecmar, Xunta de Galicia. Available online: http://www.intecmar.gal/ (accessed on 4 January 2015).
  68. EU-Harmonised Standard Operating Procedure for Determination of Lipophilic Marine Biotoxins in Molluscs by LC-MS/MS. Available online: http://www.aecosan.msssi.gob.es/CRLMB/docs/docs/metodos_analiticos_de_desarrollo/EU-Harmonised-SOP-LIPO-LCMSMS_Version5.pdf (accessed on 21 September 2018).
  69. Andrews, S. FastQC: A Quality Control Tool for High throughput Sequence Data. 2010, unpublished. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 7 October 2015).
  70. Schmieder, R.; Edwards, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics 2011, 27, 863–864. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  71. Gordon, A.; Hannon, G. Fastx-Toolkit. FASTQ/A Short-Reads Preprocessing Tools. 2010, unpublished. Available online: http://hannonlab.cshl.edu/fastx_toolkit/ (accessed on 16 October 2015).
  72. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed]
  73. Team, R. RStudio: Integrated Development for R; RStudio, Inc.: Boston, MA, USA, 2015; Available online: http://www.rstudio.com (accessed on 27 October 2015).
  74. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. B Methodol. 1995, 57, 289–300. [Google Scholar]
  75. Conesa, A.; Götz, S.; García-Gómez, J.M.; Terol, J.; Talón, M.; Robles, M. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21, 3674–3676. [Google Scholar] [CrossRef] [PubMed]
  76. Götz, S.; García-Gómez, J.M.; Terol, J.; Williams, T.D.; Nagaraj, S.H.; Nueda, M.J.; Robles, M.; Talón, M.; Dopazo, J.; Conesa, A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36, 3420–3435. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Finn, R.D.; Mistry, J.; Schuster-Böckler, B.; Griffiths-Jones, S.; Hollich, V.; Lassmann, T.; Moxon, S.; Marshall, M.; Khanna, A.; Durbin, R. Pfam: Clans, web tools and services. Nucleic Acids Res. 2006, 34, D247–D251. [Google Scholar] [CrossRef] [PubMed]
  78. CIMminer. Available online: http://discover.nci.nih.gov/cimminer (accessed on 19 June 2018).
  79. Cubero-Leon, E.; Ciocan, C.M.; Minier, C.; Rotchell, J.M. Reference gene selection for qPCR in mussel, Mytilus edulis, during gametogenesis and exogenous estrogen exposure. Environ. Sci. Pollut. Res. Int. 2012, 19, 2728–2733. [Google Scholar] [CrossRef] [PubMed]
  80. Universal ProbeLibrary. Available online: https://lifescience.roche.com/en_es/brands/universal-probe-library.html#assay-design-centre (accessed on 20 July 2018).
  81. Pfaffl, M.W.; Tichopad, A.; Prgomet, C.; Neuvians, T.P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: Bestkeeper–Excel-based tool using pair-wise correlations. Biotechnol. Lett. 2004, 26, 509–515. [Google Scholar] [CrossRef] [PubMed]
  82. Xie, F.; Xiao, P.; Chen, D.; Xu, L.; Zhang, B. miRDeepFinder: A miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol. Biol. 2012, 80, 75–84. [Google Scholar] [CrossRef] [PubMed]
  83. Hellemans, J.; Mortier, G.; De Paepe, A.; Speleman, F.; Vandesompele, J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007, 8, R19. [Google Scholar] [CrossRef] [PubMed]
  84. Ruijter, J.; Ramakers, C.; Hoogaars, W.; Karlen, Y.; Bakker, O.; Van den Hoff, M.; Moorman, A. Amplification efficiency: Linking baseline and bias in the analysis of quantitative PRC data. Nucleic Acids Res. 2009, 37, e45. [Google Scholar] [CrossRef] [PubMed]
  85. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Res. 2001, 29, e45. [Google Scholar] [CrossRef] [PubMed]
  86. Pfaffl, M.W.; Horgan, G.W.; Dempfle, L. Relative expression software tool (REST©) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002, 30, e36. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Venn diagram indicating the overlaping of genes significantly up-regulated (green arrows) and down-regulated (red arrows) when DEGs from different libraries were compared.
Figure 1. Venn diagram indicating the overlaping of genes significantly up-regulated (green arrows) and down-regulated (red arrows) when DEGs from different libraries were compared.
Toxins 10 00417 g001
Figure 2. GO classification of DEGs from the digestive gland of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Overrepresented and infrarrepresented biological processes, molecular functions and cellular components are shown. Red and green bars represent the number of down- and up-regulated genes in each category, respectively. The length of the bars is determined by the number of genes identified within each subcategory.
Figure 2. GO classification of DEGs from the digestive gland of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Overrepresented and infrarrepresented biological processes, molecular functions and cellular components are shown. Red and green bars represent the number of down- and up-regulated genes in each category, respectively. The length of the bars is determined by the number of genes identified within each subcategory.
Toxins 10 00417 g002
Figure 3. GO classification of DEGs from the gill of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. The overrepresented and infrarrepresented biogical processes, molecular functions and cellular components are shown. Red and green bars represent the number of down- and up-regulated genes in each category, respectively. The length of the bars is determined by the number of genes identified within each subcategory.
Figure 3. GO classification of DEGs from the gill of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. The overrepresented and infrarrepresented biogical processes, molecular functions and cellular components are shown. Red and green bars represent the number of down- and up-regulated genes in each category, respectively. The length of the bars is determined by the number of genes identified within each subcategory.
Toxins 10 00417 g003
Figure 4. Heatmap showing expression levels of a set of annotated genes involved in the early response to DSP toxins in mussels and selected for qPCR validation. Columns represent one library each and cells depict gene expression levels based on the number of reads. MGC: library obtained from digestive glands of control mussels. MGT: library obtained from digestive glands of treated mussels. MBC: library obtained from gills of control mussels. MBT: library obtained from gills of treated mussels.
Figure 4. Heatmap showing expression levels of a set of annotated genes involved in the early response to DSP toxins in mussels and selected for qPCR validation. Columns represent one library each and cells depict gene expression levels based on the number of reads. MGC: library obtained from digestive glands of control mussels. MGT: library obtained from digestive glands of treated mussels. MBC: library obtained from gills of control mussels. MBT: library obtained from gills of treated mussels.
Toxins 10 00417 g004
Figure 5. Relative transcript levels for each selected gene of digestive gland of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Blue bars: control samples. Green bars: samples treated with 100,000 cells/L for 48 h (mean ± SE). NRQ: Normalized Relative Quantification. n = 4. * indicates significant differences to control according to Mann-Whitney’s U-test (p-value < 0.05).
Figure 5. Relative transcript levels for each selected gene of digestive gland of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Blue bars: control samples. Green bars: samples treated with 100,000 cells/L for 48 h (mean ± SE). NRQ: Normalized Relative Quantification. n = 4. * indicates significant differences to control according to Mann-Whitney’s U-test (p-value < 0.05).
Toxins 10 00417 g005
Figure 6. Relative transcript levels for each validated candidate gene of gill of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Blue bars: control samples. Green bars: samples treated with 100,000 cells/L for 48 h (mean ± SE). NRQ: Normalized Relative Quantification. n = 4. * indicates significant differences to control in Mann-Whitney’s U-test (p-value < 0.05).
Figure 6. Relative transcript levels for each validated candidate gene of gill of the mussel M. galloprovincialis exposed to the DSP toxin-producing dinoflagellate P. lima. Blue bars: control samples. Green bars: samples treated with 100,000 cells/L for 48 h (mean ± SE). NRQ: Normalized Relative Quantification. n = 4. * indicates significant differences to control in Mann-Whitney’s U-test (p-value < 0.05).
Toxins 10 00417 g006
Figure 7. Experimental design diagram. Mussels from rocky shores were acclimated to laboratory conditions and subsequently exposed to 100,000 cells/L of P. lima for 48 h. Afterwards, gills and digestive gland were used for RNA extraction. RNA from 3 individuals was pooled for library construction and sequencing. MGC: RNA pool obtained from digestive glands of control mussels. MGT: RNA pool obtained from digestive glands of treated mussels. MBC: RNA pool obtained from gills of control mussels. MBT: RNA pool obtained from gills of treated mussels.
Figure 7. Experimental design diagram. Mussels from rocky shores were acclimated to laboratory conditions and subsequently exposed to 100,000 cells/L of P. lima for 48 h. Afterwards, gills and digestive gland were used for RNA extraction. RNA from 3 individuals was pooled for library construction and sequencing. MGC: RNA pool obtained from digestive glands of control mussels. MGT: RNA pool obtained from digestive glands of treated mussels. MBC: RNA pool obtained from gills of control mussels. MBT: RNA pool obtained from gills of treated mussels.
Toxins 10 00417 g007
Table 1. Summary of reference transcriptome assembly for M. galloprovincialis.
Table 1. Summary of reference transcriptome assembly for M. galloprovincialis.
Total number of contigs95,702L251682 bp
Total length71,623.079 KbN5021,152
Maximum contig length16,082 KbL501062 bp
Minimum contig length102 pbN7542,376
Average contig length748 bpL75668 bp
N257537%GC33.20%
Table 2. List of the 25 putative top up-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the digestive gland of M. galloprovincialis.
Table 2. List of the 25 putative top up-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the digestive gland of M. galloprovincialis.
Sequence IDDescriptionLength (bp)baseMeanLog2FCFCp-ValueAdjusted p-Value
ci|000006456|Bact|Sample_MGT2|2cytochrome c oxidase subunit 1, partial91023,389.227.09136.291.96 × 101121.70 × 10−107
ci|000001182|Bact|Sample_MBT2|2* ATP-synt_A5783975.187.21147.766.65 × 10 −661.92 × 10−61
Contig39610ribosomal protein L23a, partial11662205.277.92243.023.89 × 10−436.71 × 10−39
ci|000005084|Bact|Sample_MGT2|2cytochrome c oxidase subunit I184812,611.716.2375.081.91 × 10−372.75 × 10−33
Contig34888NA529608.038.47355.724.04 × 10−354.98 × 10−31
ci|000015505|Bact|Sample_MGT1|2NA5881046.845.0533.229.05 × 10−359.78 × 10−31
ci|000014133|Bact|Sample_MGT2|2* Glyco_hydro_16949411.926.95124.043.09 × 10−342.97 × 10−30
Contig22742NA1165944.386.6197.634.80 × 10−334.14 × 10−29
Contig33832Kazal-like serine protease inhibitor domain-containing protein507514.556.70103.846.63 × 10−335.20 × 10−29
ci|000004031|Bact|Sample_MGT1|2* Porin_31024888.455.1936.443.14 × 10−312.09 × 10−27
ci|000016700|Bact|Sample_MGT1|2† COX1_MYTED75012,157.036.86116.522.32 × 10−301.43 × 10−26
ci|000022316|Bact|Sample_MBT1|2* Ribosomal_L23340280.826.68102.352.26 × 10−281.15 × 10−24
Contig17884PREDICTED: 60 kDa SS-A/Ro ribonucleoprotein1726227.655.5145.576.21 × 10−252.68 × 10−21
ci|000012420|Bact|Sample_MGT2|2NA560761.819.971004.612.04 × 10−248.01 × 10−21
ci|000010593|Bact|Sample_MBC1|2* Ribosomal_L7Ae458403.478.29313.916.26 × 10−232.16 × 10−19
ci|000001186|Bact|Sample_MGT2|2NA316209.876.2676.557.11 × 10−232.36 × 10−19
ci|000001089|Bact|Sample_MGT2|2NA1106224.575.3540.771.03 × 10−203.08 × 10−17
Contig35276NA420188.917.45174.326.35 × 10−201.77 × 10−16
Contig38903* Myticin-prepro506510.919.28623.481.02 × 10−192.68 × 10−16
ci|000022507|Bact|Sample_MBT2|2* Ribosomal_S91568918.245.2036.823.34 × 10−197.79 × 10−16
ci|000000480|Bact|Sample_MGT2|2* Astacin858563.336.67101.714.63 × 10−191.00 × 10−15
ci|000018470|Bact|Sample_MGT1|2* Lectin_C633186.737.90238.405.47 × 10−191.13 × 10−15
ci|000004710|Bact|Sample_MBC2|2NA1395219.146.0666.621.35 × 10−182.66 × 10−15
ci|000008308|Bact|Sample_MGT1|2NA464637.375.0432.921.82 × 10−183.49 × 10−15
ci|000004147|Bact|Sample_MGT2|2NA622507.896.81112.022.47 × 10−184.65 × 10−15
FC: Fold Change. NA: No gene annotation for the transcript. * Pfam result: protein containing the specified domain. † BlastUniProt result.
Table 3. List of the 25 putative top down-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the digestive gland of M. galloprovincialis.
Table 3. List of the 25 putative top down-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the digestive gland of M. galloprovincialis.
Sequence IDDescriptionLength (bp)baseMeanLog2FCFCp-ValueAdjusted p-Value
ci|000007816|Bact|Sample_MGC1|2NA5397514.26−7.48−178.701.66 × 10−827.17 × 10−78
Contig22552NA622131,251.75−6.93−121.525.49 × 10−601.18 × 10−55
Contig26868NADH dehydrogenase subunit 5, partial7191916.17−6.96−124.622.16 × 10−321.56 × 10−28
Contig2813540S ribosomal protein S10-like559406.63−5.12−34.815.46 × 10−293.14 × 10−25
Contig30578* DUF10825293132.77−9.97−1005.572.05 × 10−281.11 × 10−24
Contig28105* SRCR1419329.84−7.74−213.172.60 × 10−261.25 × 10−22
ci|000000372|Bact|Sample_MGC1|2NA 723851.06−6.15−70.817.19 × 10−263.27 × 10−22
ci|000009048|Bact|Sample_MBC1|2NA 703359.37−8.39−334.641.14 × 10−244.67 × 10−21
Contig26906NA530286.75−8.12−279.101.47 × 10−235.52 × 10−20
ci|000018684|Bact|Sample_MBC2|2* Cytochrom_B_N_26432708.10−7.79−221.905.93 × 10−232.13 × 10−19
ci|000000728|Bact|Sample_MGC2|2NA768510.67−9.66−810.788.67 × 10−232.77 × 10−19
Contig29976uncharacterized protein LOC567525 isoform X1/* Fibrinogen_C1089167.09−6.83−114.106.31 × 10−211.95 × 10−17
ci|000000734|Bact|Sample_MGC2|2* Zona_pellucida1185255.85−9.02−518.064.16 × 10−201.20 × 10−16
Contig26843NA984230.17−4.23−18.767.68 × 10−202.07 × 10−16
ci|000002253|Bact|Sample_MGC2|2PREDICTED: GTPase IMAP family member 7/* AIG11188255.05−7.83−228.111.33 × 10−193.37 × 10−16
ci|000003979|Bact|Sample_MGC1|2NA1136233.95−8.88−471.602.94 ×10−197.25 × 10−16
ci|000021317|Bact|Sample_MBC1|2NA 8341818.89−6.97−125.283.04 × 10−197.29 × 10−16
ci|000008655|Bact|Sample_MGC2|2NA 3741501.18−8.07−267.934.48 × 10−191.00 × 10−15
ci|000004674|Bact|Sample_MGC2|2Perlucin660174.13−5.61−48.704.57 × 10−191.00 × 10−15
ci|000023153|Bact|Sample_MBT1|2* COX16054253.61−6.85−115.035.47 × 10−191.13 × 10−15
ci|000001983|Bact|Sample_MBC2|2* KOW6074924.53−2.78−6.888.86 × 10−191.78 × 10−15
ci|000005149|Bact|Sample_MGC1|2* TIG612148.86−6.00−64.192.90 × 10−175.11 × 10−14
ci|000009215|Bact|Sample_MGC2|2* Glyco_hydro_10946136.75−7.11−138.624.08 × 10−177.06 × 10−14
Contig28020NA570166.31−8.45−350.564.96 × 10−178.24 × 10−14
ci|000015516|Bact|Sample_MBT2|2* Ribosomal_L2213381356.05−4.92−30.371.19 × 10−161.87 × 10−13
FC: Fold Change. NA: No gene annotation for the transcript. * Pfam result: protein containing the specified domain.
Table 4. List of the 25 putative top up-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the gill of M. galloprovincialis.
Table 4. List of the 25 putative top up-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the gill of M. galloprovincialis.
Sequence IDDescriptionLength (bp)baseMeanLog2FCFCp-ValueAdjusted p-Value
ci|000029194|Bact|Sample_MBT1|2* EF-hand_1 and 75083570.449.76868.307.36 × 10−986.48 × 10−93
ci|000006043|Bact|Sample_MBT1|2NA 4714432.186.0566.474.32 × 10−601.90 × 10−55
ci|000001929|Bact|Sample_MBT2|2NA690803.537.61195.154.55 × 10−491.00 × 10−44
ci|000002899|Bact|Sample_MBT1|2NA7791066.249.48715.099.38 × 10−401.65 × 10−35
Contig35833NA944520.968.45350.516.42 × 10−296.28 × 10−25
ci|000022507|Bact|Sample_MBT2|2NADH dehydrogenase subunit 61568475.013.9014.881.10 × 10−289.70 × 10−25
ci|000017597|Bact|Sample_MBT1|2* Antistasin795266.997.13140.361.22 × 10−289.79 × 10−25
ci|000007496|Bact|Sample_MBT2|2NA745425.526.3984.074.01 × 10−282.94 × 10−24
ci|000020755|Bact|Sample_MBT2|2† NU4M_MYTED1483849.565.1535.441.44 × 10−269.05 × 10−23
ci|000025759|Bact|Sample_MBT2|2NA2007257.037.56188.101.81 × 10−261.06 × 10−22
Contig39610* Ribosomal_L231166453.856.4687.803.94 × 10−262.16 × 10−22
Contig15942NA482916.445.8457.296.24 × 10−263.23 × 10−22
ci|000001411|Bact|Sample_MGC2|2* HSBP1580421.114.5523.411.82 × 10−258.92 × 10−22
ci|000005084|Bact|Sample_MGT2|2* COX118483923.983.9815.753.81 × 10−251.77 × 10−21
ci|000003417|Bact|Sample_MBT2|2* Phospholip_A2_1562203.186.0867.551.10 × 10−234.05 × 10−20
Contig20144NA2258205.663.8114.014.48 × 10−231.58 × 10−19
ci|000019916|Bact|Sample_MBT1|2NA768241.647.15141.881.42 × 10−224.45 × 10−19
ci|000001302|Bact|Sample_MBT1|2NA675281.746.6197.651.36 × 10−214.00 × 10−18
ci|000018492|Bact|Sample_MBT1|2NA1552159.565.5045.191.59 × 10−214.36 × 10−18
Contig13066Calcyphosin-like protein2325853.763.6812.852.01 × 10−215.35 × 10−18
ci|000018122|Bact|Sample_MBT2|2* HYR and TMEM15433212042.074.3420.272.30 × 10−215.95 ×10−18
Contig12937† RS27L_HUMAN2183321.157.71210.043.66 × 10−219.21 × 10−18
ci|000000451|Bact|Sample_MGT2|2NA655587.432.927.571.42 × 10−203.38 × 10−17
ci|000003122|Bact|Sample_MBT2|2NA1513215.457.54185.874.85 × 10−201.12 × 10−16
Contig40138NA584136.636.1772.157.64 × 10−201.72 × 10−16
FC: Fold Change. NA: No gene annotation for the transcript. * Pfam result: protein containing the specified domain. † Blast UniProt result.
Table 5. List of the 25 putative top down-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the gill of M. galloprovincialis.
Table 5. List of the 25 putative top down-regulated genes (ordered by p-value) in response to early concentrations of DSP toxins in the gill of M. galloprovincialis.
Sequence IDDescriptionLength (bp)baseMeanLog2FCFCp-ValueAdjusted p-Value
ci|000007038|Bact|Sample_MBC2|2low-density lipoprotein receptor-related protein 8 isoform X16891321.99−9.45−700.321.06 × 10−533.11 × 10−49
Contig3681NA8951561.88−10.25−1216.979.98 × 10−371.46 × 10−32
Contig11592NA798652.29−9.55−752.014.28 × 10−315.38 × 10−27
Contig1183NA581802.66−4.76−27.109.28 × 10−301.02 × 10−25
Contig8105NA1717199.06−6.43−86.513.62 × 10−272.45 × 10−23
ci|000015242|Bact|Sample_MGT2|2NA663429.13−7.24−151.364.36 × 10−251.92 × 10−21
ci|000005973|Bact|Sample_MBC1|2NA1860257.90−4.59−24.075.34 × 10−252.24 × 10−21
Contig10936* Oxidored_q1279718,810.53−1.67−3.183.68 × 10−241.47 × 10−20
ci|000000312|Bact|Sample_MBC1|2NA972208.97−7.32−159.954.63 × 10−241.77 × 10−20
Contig6277NA7022009.78−10.35−1303.406.01 × 10−232.04 × 10−19
ci|000016192|Bact|Sample_MBC1|2* Ldl_recept_a and PRKCSH-like9464008.01−2.65−6.269.71 × 10−233.17 × 10−19
Contig3876Predicted protein536295.82−8.43−344.345.01 × 10−221.52 × 10−18
Contig4774neurocalcin homolog1267339.45−4.83−28.391.49 × 10−214.22 × 10−18
ci|000000823|Bact|Sample_MBC2|2NA803325.92−3.99−15.911.06 × 10−202.59 × 10−17
Contig6059NA4865331.78−1.93−3.828.13 × 10−201.76 × 10−16
Contig7283cytochrome c oxidase subunit I2879132.24−5.09−34.069.51 × 10−201.99 × 10−16
ci|000015433|Bact|Sample_MBC1|2cytochrome c oxidase subunit I113617,057.57−3.52−11.493.00 × 10−195.74 × 10−16
ci|000004320|Bact|Sample_MBC1|2* Lipoxygenase1950564.59−9.46−703.803.48 × 10−196.24 × 10−16
ci|000001144|Bact|Sample_MBC2|2NA434221.23−8.61−389.884.06 × 10−197.15 × 10−16
ci|000008127|Bact|Sample_MBC2|2NA2503364.67−10.03−1043.767.20 × 10−191.22 × 10−15
ci|000005247|Bact|Sample_MBC1|2NA1246153.48−6.35−81.801.26 × 10−182.09 × 10−15
ci|000001610|Bact|Sample_MBC1|2NA1057341.22−9.94−982.301.93 × 10−183.11 × 10−15
ci|000000874|Bact|Sample_MBC2|2NA697153.63−5.02−32.425.21 × 10−187.77 × 10−15
ci|000002263|Bact|Sample_MBC1|2* Pfam-B_56821222151.49−4.77−27.376.23 × 10−188.99 × 10−15
ci|000003990|Bact|Sample_MGT2|2NA521310.73−9.76−866.012.46 × 10−173.23 × 10−14
FC: Fold Change. NA: No gene annotation for the transcript. * Pfam result: protein containing the specified domain.
Table 6. Pfam families significantly enriched (False Discovery Rate (FDR) adjusted p-value < 0.1) with seven or more differentially expressed genes in digestive gland.
Table 6. Pfam families significantly enriched (False Discovery Rate (FDR) adjusted p-value < 0.1) with seven or more differentially expressed genes in digestive gland.
CategoryNumber of Genesp-Value
PF04548.11//AIG1260.00248912
PF01926.18//MMR_HSR1250.0029543
PF00059.16//Lectin_C210.01366918
PF00100.18//Zona_pellucida160.00403746
PF13499.1//EF-hand_7140.00221868
PF13405.1//EF-hand_6140.00355134
PF00036.27//EF-hand_1130.00065889
PF13202.1//EF-hand_5130.02872925
PF13833.1//EF-hand_8120.02032022
PF00361.15//Oxidored_q1100.00489835
PF00119.15//ATP-synt_A80.00023995
PF10690.4//Myticin-prepro80.02237525
PF07679.11//I-set70.04744078
Table 7. Pfam families significantly enriched (FDR adjusted p-value < 0.1) with seven or more differentially expressed genes in gill.
Table 7. Pfam families significantly enriched (FDR adjusted p-value < 0.1) with seven or more differentially expressed genes in gill.
CategoryNumber of Genesp-Value
PF00386.16//C1q365.2 × 10−8
PF00036.27//EF-hand_1310.00035296
PF13499.1//EF-hand_7290.00014495
PF13405.1//EF-hand_6278.61 × 10−5
PF00147.13//Fibrinogen_C250.01665835
PF13202.1//EF-hand_5230.00079724
PF13833.1//EF-hand_8200.00015502
PF10690.4//Myticin-prepro130.01435613
PF00361.15//Oxidored_q1130.03222834
PF07679.11//I-set90.00010238
PF09458.5//H_lectin90.00621282
PF01607.19//CBM_1490.02592731
PB002965//Pfam-B_296590.03289021
PF13895.1//Ig_280.00039907
PF00119.15//ATP-synt_A80.01065053
PF13927.1//Ig_370.00090571
PF00092.23//VWA70.00272518
PF07686.12//V-set70.01729404
PF03281.9//Mab-2170.03056683
Table 8. Primers used in the real-time qPCR validation.
Table 8. Primers used in the real-time qPCR validation.
Gene NameAbbreviationReferenceEAmplicon Size (bp)Tm (°C)Primers 5′→3′
TropomyosinTPMab000907.11.9067F-55.3
R-57.1
F-GATGCTGAAAATCGTGCAAC
R-CGGTCTACTTCTTTTTGCAACTT
Ribosomal proteins S4rpS4Lozano et al. (2015)1.83138F-58.8
R-60.3
F-TGGGTTATCGAGGGCGTAG
R-TCCCTTAGTTTGTTGAGGACCTG
18S ribosomal RNA18SL33452.11.8660F-58.3
R-55.9
F-CCTGGAAAGGTCGGGTAAC
R-AATTACAAGCCCCAATCCCTA
18S ribosomal RNA18S-L33448Cubero-Leon et al. (2012)1.79114F-56.3
R-56.0
F-CATTAGTCAAGAACGAAAGTCAGAG
R-GCCTGCCGAGTCATTGAAG
Glyceraldehyde 3-phosphate-dehydrogenaseGAPDHLozano et al. (2015)1.92114F-59.4
R-58.4
F-AGGAATGGCCTTCAGGG
R-TCAGATGCTGCTTTAATGGCTG
Elongation Factor 1EF1Suarez-Ulloa et al. (2013)1.89106F-55.8
R-57.0
F-CCTCCCACCATCAAGACCTA
R-GGCTGGAGCAAAGGTAACAA
Big defensin 2BD2Contig378961.83110F-60.3
R-59.3
F-TCTGAGCAGGGAGTATCAACAG
R-TGGACAAAACAGCTACTAACAAGG
NADH dehidrogenase subunit 5NADH5Contig242661.8690F-53.7
R-56.5
F-GCAGTCATGCGCAAAAAG
R-ACCCGGTACAAATATGGCTAAA
Dynactin-subunit-6-likeDYNAContig145511.8960F-58.9
R-58.9
F-AGTATTCTCAGGCATGGTTTCTG
R-GGTTGTATAATTGGAGGCATGTG
GIY-YIG domain containing proteinGIY-YIG DCci|000000744|Bact|Sample_MBC1|21.8370F-57.6
R-55.3
F-AATCTACCAATTGCTTGTCTGTCA
R-CGAAACGTAGTGTGCGAAAA
KAZAL domain containing proteinKAZAL DCContig338321.9160F-53.2
R-60.3
F-ATAATCGGCAGTGCAAAACA
R-TTCCTTACTGAGTCAGTCG
Cytosolic phospholipase A-2 likeCPLA2ci|000001655|Bact|Sample_MBT1|21.8073F-61.6
R-57.1
F-CCTGTACTGTGAGATTAGGTTATTGC
R-CAGAAGGTTATTGACCGAAAGAA
Arachidonate 15-lipoxygenase B-likeALOX15Bci|000023941|Bact|Sample_MBT2|21.8194F-58.5
R-55.9
F-TGTTGTGAGTGAAGCAATAACTCTAA
R-CGGAATAAATCG AGAGAACCA
Alpha-L-fucosidase-likeFUCAci|000010451|Bact|Sample_MBT1|21.8774F-61.0
R-55.3
F-GGAATTCCAGTAGGAATCAGTAGC
R-TGGTAAATGCATACAAACCTGAA
H_Lectin domain containing proteinH_Lectin DCContig193411.8573F-56.5
R-55.3
F-CCCTTCTTTGCTTTAGATGCTT
R-TTGATGGCCAGATTACGACA
Fibrinogen_C domain containing proteinFibrinogen_C DCci|000024772|Bact|Sample_MBC1|21.8667F-57.3
R-59.4
F-AAGGTTGTCTCCAGCGTTTC
R-CGGTGATGCCTCTACCAACT
E: primer efficiency; F: forward; R: reverse.
Table 9. Rank of six candidate reference genes for real-time qPCR calculated by Normfinder and BestKeeper analyses.
Table 9. Rank of six candidate reference genes for real-time qPCR calculated by Normfinder and BestKeeper analyses.
RankNormfinderStabilityBestKeeperSDr
1rpS40.07rpS40.460.732
2TPM0.17TPM0.500.448
3GAPDH0.20GAPDH0.640.669
418S0.3718S0.710.827
518S-L334480.7618S-L334481.08
6EF11.78EF12.91
SD: standard deviation; r: coefficient of correlation between each gene and the BestKeeper index.

Share and Cite

MDPI and ACS Style

Prego-Faraldo, M.V.; Martínez, L.; Méndez, J. RNA-Seq Analysis for Assessing the Early Response to DSP Toxins in Mytilus galloprovincialis Digestive Gland and Gill. Toxins 2018, 10, 417. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins10100417

AMA Style

Prego-Faraldo MV, Martínez L, Méndez J. RNA-Seq Analysis for Assessing the Early Response to DSP Toxins in Mytilus galloprovincialis Digestive Gland and Gill. Toxins. 2018; 10(10):417. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins10100417

Chicago/Turabian Style

Prego-Faraldo, María Verónica, Luisa Martínez, and Josefina Méndez. 2018. "RNA-Seq Analysis for Assessing the Early Response to DSP Toxins in Mytilus galloprovincialis Digestive Gland and Gill" Toxins 10, no. 10: 417. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins10100417

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop