Next Article in Journal
Integrating Engineering, Manufacturing, and Regulatory Considerations in the Development of Novel Antivenoms
Next Article in Special Issue
Accumulation and Biotransformation of Dinophysis Toxins by the Surf Clam Mesodesma donacium
Previous Article in Journal
An Interview with Cesare Montecucco
Previous Article in Special Issue
Effect of Suspended Particulate Matter on the Accumulation of Dissolved Diarrhetic Shellfish Toxins by Mussels (Mytilus galloprovincialis) under Laboratory Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Detoxification- and Immune-Related Transcriptomic Analysis of Gills from Bay Scallops (Argopecten irradians) in Response to Algal Toxin Okadaic Acid

1
Laboratory of Aquatic Nutrition and Ecology, College of Animal Science and Technology, Nanjing Agricultural University, Weigang Road 1, Nanjing 210095, China
2
Laboratory of Aquatic Biomedicine, College of Veterinary Medicine and Research Institute for Veterinary Science, Seoul National University, Seoul 151742, Korea
*
Author to whom correspondence should be addressed.
Toxins 2018, 10(8), 308; https://doi.org/10.3390/toxins10080308
Submission received: 29 May 2018 / Revised: 24 July 2018 / Accepted: 26 July 2018 / Published: 28 July 2018
(This article belongs to the Special Issue Dinophysis Toxins: Distribution, Fate in Shellfish and Impacts)

Abstract

:
To reveal the molecular mechanisms triggered by okadaic acid (OA)-exposure in the detoxification and immune system of bay scallops, we studied differentially-expressed genes (DEGs) and the transcriptomic profile in bay scallop gill tissue after 48 h exposure to 500 nM of OA using the Illumina HiSeq 4000 deep-sequencing platform. De novo assembly of paired-end reads yielded 55,876 unigenes, of which 3204 and 2620 genes were found to be significantly up- or down-regulated, respectively. Gene ontology classification and enrichment analysis of the DEGs detected in bay scallops exposed to OA revealed four ontologies with particularly high functional enrichment, which were ‘cellular process’ (cellular component), ‘metabolic process’ (biological process), ‘immune system process’ (biological process), and ‘catalytic process’ (molecular function). The DEGs revealed that cyclic AMP-responsive element-binding proteins, acid phosphatase, toll-like receptors, nuclear erythroid 2-related factor, and the NADPH2 quinone reductase-related gene were upregulated. In contrast, the expression of some genes related to glutathione S-transferase 1, C-type lectin, complement C1q tumor necrosis factor-related protein, Superoxide dismutase 2 and fibrinogen C domain-containing protein, decreased. The outcomes of this study will be a valuable resource for the study of gene expression induced by marine toxins, and will help understanding of the molecular mechanisms underlying the scallops’ response to OA exposure.
Key Contribution: The Illumina platform was used for the first time to analyse gene expression in the gills of bay scallop exposed to OA. Detoxification- and immune-related genes and pathway enrichment following OA exposure were detected.

1. Introduction

Bivalves are among the most important commercially exploited marine species in China, sharing 75–80% of the total output of aquatic products in recent years [1]. Owing to their filter-feeding and sessile habits, worldwide distribution, and diversity of aquatic environments, bivalves are widely used as marine pollution bioindicators [2]. Scallop fisheries are mainly distributed along coastal areas of Japan, Korea, and North China [3]. In addition to their economic value, bivalves have always been studied as model species in toxicological investigation and as sentinel species in environmental monitoring programmes [4].
The frequent appearance of toxin-producing harmful algal blooms (HABs) in marine environments is a well-known worldwide problem [5]. HABs are well known for their potential to induce ecological damage, risk human health, and cause adverse effects to living marine resources [6,7]. Moreover, these HABs threaten aquaculture industries and may have deleterious effects on public health [8], because their phycotoxins may cause mass mortality of cultivated animals [9]. Shellfish toxins are the main marine phycotoxin, which includes amnaesic shellfish poisoning (ASP)-, paralytic shellfish poisoning (PSP)-, neurotoxic shellfish poisoning (NSP)-, diarrhetic shellfish poisoning (DSP)-, and azaspiracid shellfish poisoning (AZP) toxins [10]. These toxins may be taken up by humans eating shellfish contaminated with them, and lead to a series of neurological and gastrointestinal syndromes [6,7]. Okadaic acid (OA), representative of the DSP toxins, can be produced by species of the genera Dinophysis and Prorocentrum [11,12], and be accumulated in the shellfish adipose tissue [13]. This is the primary cause of acute DSP intoxication of human consumers, and harvesting bans causing huge economic losses to the shellfish aquaculture industry [14]. For example, Mouratidou et al. [15] reported maximum concentrations of 36 μg OA eq/g hepatopancreas in mussels from Thermaikos Gulf, Greece. OA is capable of binding to the active sites of protein phosphatases [16], inhibiting their activity and inducing tumorigenic and apoptotic processes [14,17]. Finally, it can lead to the hyperphosphorylation of many cellular proteins, metabolic deregulation, and genotoxic and cytotoxic damage [18]. When organisms are exposed to xenobiotics, short-term responses, such as changes in their immune response, and long-term effects on other biological parameters, including growth, ingestion and reproduction rates, and other metabolic processes may be observed [19]. Earlier investigations revealed that OA or P. lima exposure could induce haemocyte function damage and reduced survival in Ruditapes decussatus [20]. Huang et al. [11] reported that OA-producing P. lima caused oxidative stress, disorganization of cytoskeletons, and metabolic disturbance in mussels. In a previous work, we studied the toxic effects of OA exposure, up to 48 h, in bay scallops (Argopecten irradians). These included changes in glutathione (GSH), reactive oxygen species (ROS), malondialdehyde (MDA), and nitric oxide (NO) contents; lysozyme, acid phosphatase (ACP), lactate dehydrogenase (LDH), alkaline phosphatase (ALP), and superoxide dismutase (SOD) activity; total haemocyte counts (THC) and haemolymph total protein levels [8,12]. Overall, our previous work demonstrated that OA exposure increased oxidative stress, disrupted metabolism, modulated the immune response, and was toxic to physiological function in A. irradians. There are two resistance mechanisms that may counteract the effects of DSP in shellfish: detoxification pathways for the biotransformation or elimination of phycotoxins, and antioxidant metabolism to neutralize ROS induced by DSP exposure [21,22,23]. However, how scallops respond to OA toxicity, and the details of their detoxification process during acute OA exposure remain unclear, particularly the integral response at the transcriptional level. An understanding of the effects of OA exposure on the bay scallop is essential to establish effective measures to estimate its toxic potential. However, owing to the constraint of related genomic resources, a better understanding of the genetic and molecular mechanisms underlying the bay scallop response to sublethal concentrations of OA is yet to be elucidated.
De novo sequencing is an effective tool to obtain whole scallop transcriptome information. In this regard, the relatively low-cost/high-output Illumina HiSeq™ 4000 sequencing platform has found increasingly widespread use [24], having been applied to a growing number of aquatic organisms, including Oryzias melastigma [25], Crassostrea gigas [26], and Chlamys farreri [27], to study their responses to environmental stressors. Therefore, the aim of the present study was to obtain a better understanding of the molecular response of the bay scallop after exposure to OA. We specifically focused on the gill tissue of A. irradians, following exposure to 500 nM of OA for up to 48 h, since our previous studies found that this toxin induced oxidative stress, modulated the immune response, and was toxic to physiological function in A. irradians [8,12]. The gill was chosen as the target organ because it is the first organ in contact with OA during filtration [21]. Gills act as a defence barrier, because they play a crucial role in the filtration of suspended matter. Further, the gill was previously found to be directly affected by contact with toxic algae [21], and to have a high expression of putative immune-related genes [28]. Digital gene expression (DGE) analysis was performed with the Illumina HiSeq™ 4000 sequencing system, and then quantitative real-time PCR was conducted to verify differentially expressed genes (DEGs), which were selected according to the DGE analysis. The aim of the present work was to reveal the transcript abundance to facilitate a network of bay scallop genes enriched to regulate toxicological responses to OA exposure.

2. Results

2.1. Analysis of DGE Libraries

Two DGE libraries comprising DNA from the gills of control and OA-exposed scallops were analysed using the Illumina Hiseq 4000 sequencing system. We removed adaptors from the reads, poly N, and low-quality reads from the raw data, and then generated 9.14 Gb of totally clean bases, comprising 45.92 and 45.92 Mb clean reads for control and OA-exposed cDNA libraries, respectively. The Q20 and GC percentages of the clean reads in the two cDNA libraries were 98.21% and 98.17% and 39.13% and 39.24%, for control and OA-exposed cDNA libraries, respectively (Table S1). Clean sequences from each library were assembled by the Trinity tool, thereby producing a total of 78,510 and 77,330 transcripts in the control and OA-exposed groups, respectively, which had mean sizes of 675 with N50s of 1234 for the control group and 733 bp with N50s of 1451 bp for the OA-exposed groups, respectively (Table S1). Finally, 55 876 unigenes were further merged by transcript sets from the two libraries (Table 1). The size distribution of the unigenes was as follows: 67.58% (37,759) were between 300 and 1000 bp; 20.54% (11,477) were between 1000 and 3000 bp; and 6.24% (3488) had lengths greater than 3000 bp in length, as shown in Figure 1.

2.2. Functional Annotation and Species Distribution

After assembly, functional annotation was carried out through seven functional databases for unigenes. A total of 49.31% of the total unigenes (27,555 unigenes) were annotated, of which 24,521 unigenes (43.88%) were aligned to the Nr database; 10,466 unigenes (18.73%) to Nt; 19,220 unigenes (34.40%) to Swiss-Prot, 18,523 unigenes (33.15%) to Kyoto Encyclopedia of Genes and Genomes (KEGG); 8800 (15.75%) unigenes to Clusters of Orthologous Group (COG); 18,533 (33.17%) unigenes to Interpro; and 4027 unigenes (7.21%) to Gene Ontology (GO), respectively.
The distribution of annotated species was statistically analysed with NR annotation, as shown in Figure 2. For functional classification, 15 186 unigenes were totally annotated to the COG database (Figure 3). The most frequently functional classifications were the following: 20.70% (3143) accounted for general function; 8.52% (1294) related to recombinant and repair; translation, 8.49% (1289); transcription, 6.63% (1007); post-translational-modification-related, 6.26% (950); cell-cycle-control-related, 5.64% (856), and signal-transduction-related, 5.39% (819).

2.3. Differential Gene Expression Analysis

The unigene expression levels were calculated using the Fragments Per Kilobase Million (FPKM) method (Figure 4 and Figure 5) to identify the genes’ differential expression between the control and OA-treated groups. A total of 5825 unigenes with different expression levels (with over two-fold changes, and false discovery rate (FDR) ≤ 0.001) between the control and OA-exposed groups were identified. Of these, 3204 were upregulated genes, while 2621 were downregulated genes (Table S2).

2.4. Enrichment and Pathway Analysis

In order to identify their function, all the DEGs were mapped to the GO database. A total of 44 functional groups in the DEGs were substantially enriched compared with the genomic background (Figure 6). Genes in the OA-exposed scallop related to the terms ‘metabolic process’, ‘cellular process’, and ‘catalytic activity’ were dominant. Biological process and cellular components were found to be the most-represented known genes, followed by molecular function.
Markedly-enriched signal transduction and metabolic pathways were identified using KEGG enrichment analysis of the DEGs. A total of 3389 DEGs were aligned at 299 pathways in the KEGG database, and 74 metabolic pathways were significantly (corrected p value < 0.05) over-represented. The pathway classification results are shown in Figure 7, and the pathway functional enrichment results in Figure 8. Among these, the expression patterns of DEGs throughout OA exposure, which involved detoxification, and immunology in mechanisms against biotoxins were further analyzed on the bases of GO and KEEG analyses. The expression of genes related to the immunology and detoxification responses such as cyclic AMP-responsive element-binding proteins, acid phosphatase, toll-like receptors, nuclear factor erythroid 2-related factor, NADPH2: quinone reductase, cytochrome P450 3A64 and 3A80 increased under exposure to OA (Table 2). In contrast, the expression of some genes related to glutathione S-transferase 1, C-type lectin, complement C1q tumor necrosis factor-related protein, Superoxide dismutase 2 and fibrinogen C domain-containing protein decreased.

2.5. Identification of Genes Related to OA-Induced Stress Response

The real-time quantitative PCR (qPCR) technique was used to detect the relative expression levels of nine genes, which are immunology-, detoxification- and antioxidant-ability-related genes with high expression, from the DGE libraries. Four of these genes were suppressed and the others were induced. The melting-curve analysis of each gene performed by qPCR suggested a single product. The qPCR results were compared with those from the DGE analysis. As shown in Figure 9, nine genes followed a concurrent trend between qPCR analysis and DGE library, and the correlation coefficient was calculated as 0.95 (p value < 0.001).

3. Discussion

Okadaic acid (OA), as a representative of DSP toxins, can accumulate in bivalves and induce diarrheic shellfish poisoning in mammals [29]. OA has been reported to be cytotoxic in several cell lines (human monocytic U-937 cells; two epithelial tumour lines, HeLa and KB; neuroblastoma cell line Neuro-2a; neuroblastoma × glioma hybrid cell line NG108-15; breast cancer cell line MCF-7) as an efficient inhibitor of serine/threonine phosphatases [30,31,32]. Earlier, we reported that OA exposure could affect a variety of innate immune responses (e.g., THC, total protein level, ALP, ACP, and lysozyme activities,) and physiological responses (e.g., SOD and LDH activity, ROS, NO and MDA and GSH content) in the haemolymph of scallops, and can even induce oxidative stress and disrupt metabolism in bay scallops [8,12], rendering them sensitive to OA exposure. Previous studies have demonstrated the adverse impacts of the toxin OA on other marine bivalves [11,20]; however, the molecular response of these bivalves to OA is not well characterized. In the light of our earlier studies, the results of this transcriptome information could improve the description of the acute toxicity of high concentrations of OA for some physiological and biochemical processes and provide directions and insights for future studies involving biotoxicity models in scallops. Moreover, in the present study, the calculation and normalization methods of both analyses are different, although they report transcript abundances as fold-changes relative to the control [1]. The RNA-seq expression values employ Reads Per Kilobase Million (RPKM) for calculation [33], while qPCR fold-change values employ the mean normalized expressions method and incorporated reference gene to calculate [34]. In the present investigation, both methods were used for transcript quantification. The same directions of change and a similar magnitude of the fold-change in abundance confirmed the accuracy and reliability of the DGE data. To our knowledge, the present investigation is the first to reveal the transcriptomic responses of scallops after OA exposure using deep-sequencing technology.
Highly conserved heat shock proteins (HSPs), including HSP60, HSP70, and HSP90, could be synthesized or secreted rapidly by cells as soon as they experience stressed [3]. Therefore, HSPs have been widely considered as effective biomarkers of exogenous stimuli or as biomonitoring tools to identify the effects of environmental pollution in aquatic animals, including bivalves [3]. Our present investigation showed that the relative expression of HSP70, which was validated by qPCR, was strongly increased in the gills of bay scallop up to 48 h exposure to OA. Similarly, a previous study revealed that upregulated HSP70 expressed transcripts were identified in the mussel Mytilus galloprovincialis after exposure to OA stress [14]. In other investigations, the detection of HSP70 by immunoblotting and expression analysis of HSP70 mRNA was used to indicate marine contamination observed following exposure to heavy metals in Dreissena polymorpha [35], to hydrocarbon in Crassostrea gigas [36], to sub-lethal concentrations of quaternium-15 in M. galloprovincialis [37], and to cadmium in the gills of Ostrea edulis [38]. Therefore, in the present investigation, the upregulation of HSP70 mRNA in the gills of bay scallops also appears to be a helpful marker for toxic effects.
The genes encoding detoxification enzymes play crucial roles in bivalves after being stimulated by a variety of exogenous stimuli, including drugs, toxicants, and chemical carcinogens [3]. Among the DEGs detected in the present study, certain detoxification-related genes were identified. The cytochrome P450 (CYP450) family is an essential family of enzymes related to the primary or phase I metabolism of xenobiotics, including pesticides and toxins [3,39]. Many exogenous stimuli may impact the metabolism, and then activate or suppress the activity of CYP450 to clean exogenous stimuli [11]. A subset of cytochrome P450 enzymes, which are linked to detoxification and resistance, were involved in transforming liposoluble toxic chemicals into hydrosoluble substances that are easily eliminated [11,40]. Our results clearly showed that OA provoked the differential expression of CYP1A5 and CYP3A24, which were downregulated, whereas CYP3A4 and CYP3A80 were upregulated. This is consistent with a previous study that reported OA-exposure-induced expression of CYP450 mussel gills, which suggests that CYP450 participates in the process of OA elimination [3]. Guo et al. also [41] reported that human recombinant cytochrome CYP3A4 could eliminate OA by generating oxidized products. Accordingly, CYP3A4 and CYP3A80 may participate in the process of accelerating the biotransformation of OA and facilitating its excretion in bay scallops when exposed to OA. ATP-binding cassette (ABC) transporters are a family of transmembrane proteins that can transport a variety of strGSTucturally diverse substrates across biological membranes in an ATP-dependent manner [11]. In mammalian tumor cells, they are responsible for a multidrug resistance phenotype. Moreover, in aquatic organisms, they are responsible for a multixenobiotic resistance phenotype by exporting xenobiotics out of the cells or by facilitating the sequestration of toxins within specialized cells or organelles, effectively segregating them away from vulnerable protein and DNA targets [11]. In our present study, we found that ABCB10, ABCC5, and ABCC1 were upregulated in bay scallops after 500 nM OA exposure. These results are consistent with a previous study showing that ABC transporters in mussels were upregulated after exposure to P. lima. Huang et al [42] also found that the expression level of a P-glycoprotein gene (P-gp), belonging to the family of ATP-binding cassette (ABC) transporters in the gills of Perna viridis, increased significantly after exposure to P. lima. These phenomena suggest the possible role of ABC transporters in OA detoxification.
Nicotinamide adenine dinucleotide phosphate-oxidases (NADPH-oxidases) are enzymes completely devoted to ROS production [43]. The family of NADPH-oxidases comprises trans-membrane proteins that transfer electrons across biological membranes. Owing to their involvement in ROS production, NADPH-oxidases play crucial roles in various physiological mechanisms which include host defence, gene expression, cellular signalling, apoptosis, and oxidative stress [44]. The NADPH oxidase is composed of six homologues of the cytochrome subunit (NOX1, NOX3, NOX4, NOX5, DUOX1, and DUOX2), and increased NOX activity also induces a series of pathologies [44]. Cai et al. [1] found that benzo[a]pyrene (BaP) exposure caused the upregulation of NADPH transcript in Chlamys farreri after three days. The findings of the present study indicated a greater abundance of NOX-3 transcripts in the gills of scallops exposed to OA, suggesting that it induces the activation of the NADPH oxidases, thereby generating more ROS and even cell damage.
The detoxification and biotransformation of exogenous compounds also rests on Phase II and Phase III reactions [1]. Glutathione-S-transferase (GST), which is a kind of Phase II enzyme, could catalyse the endogenous and exogenous compounds combining with glutathione (GSH) [1,45]. Our previous field studies have shown that GSH levels in the haemolymph of A. irradians exposed to 500 nM OA decreased sharply at 48 hpe [12]. Consistently, in the present study, the expression of GST mRNAs, including GST1, GST2, GST-A, GST-Theta-1, GST-Omeaga, and GST-Kappa, in the DGE library decreased in the gills of A. irradians exposed to OA compared to the control group. This is consistent with previous studies showing that the expression of GST-pi was significantly down-regulated in the digestive gland of M. galloprovincialis in response to toxic dinoflagellate Prorocentrum lima (1000 cells/L) for 48 h [46]. These results suggest that the expression level of GST was attenuated by 500 nM OA exposure, which weakened the detoxification or antioxidant capacity of the OA-exposed scallops. SOD is a crucial gene belonging to the antioxidant defence system. It can eliminate the ROS, which can induce lipid peroxidation processes and ultimately lead to DNA damage [47]. We previously reported that Mn SOD expression levels in the haemolymph of OA-exposed bay scallops decreased significantly after 48 hours post-exposure [8]. These observations are in agreement with the results of the present study, verified by qPCR, showing that the SOD2 expression levels in gills were downregulated after 48 h exposure to OA. However, we found that the expression of the Cu/Zn SOD mRNA was clearly induced, indicating that OA exposure could induce the expression level of Cu/Zn SOD in the gills when the scallops are exposed to up to 48 h 500 nM levels of OA. Additionally, it might be plausible that the downregulation of GST is partially compensated by the upregulation of Cu/Zn SOD, since both enzymes use the same substrate [46,48,49].
Cyclic adenosine monophosphate responsive element binding-protein (CREB) plays a pivotal role in the immune response. OA stimulation was found to enhance the levels of phosphorylated-CREB [50]. The expression of these genes is essentially regulated by the phosphorylation state of CREB, since phosphorylation is necessary for CREB to bind to the cAMP response element in the promoter of several early response genes [50]. This result is in accordance with a previous study showing that OA was able to induce CREB expression in mussels [50]. Acid phosphatase (ACP) is a kind of essential hydrolytic enzyme in phagocytic lysosomes [51]. In the present research, we found that the ACP mRNA expression increased in the gills of OA-exposed bay scallops. Nevertheless, in an earlier investigation, we demonstrated that OA exposure suppressed the ACP levels in the haemolymph of bay scallops, indicating that although OA could induce ACP expression, it might also affect the assembly, folding, or modification of the ACP, leading to a deficiency in the elimination of pathogens or phagocytized microorganisms in the OA-exposed gills.
In conclusion, we present here broader research into the OA-responsive genes, such as the Toll-like receptor, ATP-binding cassette, cyclic AMP-responsive element-binding protein, cytochrome P450 and Cu-Zn superoxide dismutase related genes, that show differential expression in the bay scallop, suggesting participation in the resistance to OA toxicity. These genes are related to a series of detoxification and immune processes in the response to OA. The present investigation not only reveals the transcriptional complexity of the response to OA stimulation in scallops, but also suggests the possibility of identifying the genes implicated in regulating the bivalves’ tolerance or the elimination of algal toxin stress. However, it remains unclear whether these immune responses are directly stimulated by abiotic factors or whether OA exposure just facilitates the opportunistic attack of pathogens present in the scallops’ microbiota [14]. Illumina next-generation sequencing technology provides a good resource to explain the immune- and detoxification-associated molecular mechanisms triggered in the bay scallop to endure the toxic effects of OA. Furthermore, it supplements and reinforces the results from our previous investigations, from which a strong cause and effect relationship between OA and the differential expression of immune- and detoxification-associated factors in the bay scallop were established. These results will be useful to develop potential countermeasures to manage the toxic effects of OA on exploited bivalve resources.

4. Materials and Methods

4.1. Maintenance of Scallops

Bay scallops A. irradians (weight: 46.02 ± 2.67 g; shell length: 60–70 mm) were procured in a wholesale market in Seoul, South Korea. To acclimate them to laboratory conditions, these scallops were kept for 2 weeks in 800-L tanks containing filtered and aerated seawater, with a temperature of 10 ± 1 °C and a salinity of 30 ± 0.1 psu [8]. They were fed with a commercial shellfish diet (Instant Algae® Shellfish Diet, Campbell, CA, USA) at a rate of approximately 1.2 × 1010 algae cells/scallop/day [8]. Half the seawater volume was daily renewed.

4.2. Okadaic Acid Exposure and RNA Extraction

In total, 120 scallops were divided in two groups, i.e., control and OA-exposure groups. Each group consisted of 60 scallops distributed in 3 replicate tanks with 20 scallops each. Okadaic acid (OA) (92–100% HPLC purified) was purchased from Sigma-Aldrich, USA and stored at 4 °C until use. To prepare the stock solution, OA was dissolved in 1 mL of dimethyl sulfoxide (DMSO; Sigma-Aldrich, St. Louis, MO, USA). The final concentration of OA in the OA-treated group was kept at 500 nM [8]. The scallops in the control group were treated with an equal volume of DMSO, with a final concentration of 0.0125‱DMSO in each tank [8]. After 48 h of OA exposure, six scallops were collected from each tank (i.e., 6 scallops × 3 replicates = 18 scallops per group) and maintained on ice. Scallop gills were dissected, stored in 1 mL TRIzol reagent (Invitrogen, Waltham, MA, USA), and frozen at −80 °C until use. Samples from 6 scallops were pooled for each replicate for RNA extraction [14]. Total RNA was extracted using TRIzol (Invitrogen, Waltham, MA, USA) following the manufacturer’s instruction. Agarose gels (1%) electrophoresis was preformed to monitor the RNA contamination and degradation [52]. The RNA purity and contamination was checked with a NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA) and a Qubit® RNA Assay Kit and a Qubit® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) respectively [52]. RNA integrity was measured using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) [52].

4.3. Library Preparation and Illumina Sequencing

After treating the total RNA extract sample with DNase I, 200 ng were purified with oligo-dT beads. In brief, total RNA and RNA Purification Beads (Illumina, San Diego, CA, USA) were incubated and resuspended in Elution Buffer (Illumina, San Diego, CA, USA). The mRNA was eluted from the beads, and then incubated to rebind the beads after adding Bead Binding Buffer (Illumina, San Diego, CA, USA). Finally, Fragment Buffer was used to fragment poly (A)-containing mRNA into small pieces. The mRNA fragments were used as templates during the cDNA synthesis. First-strand cDNA was synthesized by reverse transcription using First Strand Master Mix (Illumina, San Diego, CA, USA) and Super Script II (Invitrogen, Waltham, MA, USA). The conditions for the reverse transcription reaction were: 25 °C for 10 min; 42 °C for 50 min and 70 °C for 15 min. Next, the second-strand cDNA was synthesized at 16 °C for 1 h using Second Strand Master Mix (Illumina, San Diego, CA, USA). Then, the ds cDNA was separated from the second strand using AMPure XP beads (Agencourt, Beverly, MA, USA). The remaining overhangs were converted into blunt ends using an End Repair Mix. Next, after adding the A-Tailing Mix, the mixture was incubated at 37 °C for 30 min. The Adenylate 3′Ends DNA, RNA Index Adapter and Ligation Mix were combined and the ligate reaction incubated at 30 °C for 10 min to perform the A ligation reaction. AMPure XP Beads were used to purify the end-repaired DNA. In order to enrich the cDNA fragments, several rounds of PCR amplification were performed by adding PCR Primer Cocktail and PCR Master Mix. The AMPure XP Beads were used to purify the library fragments to select cDNA fragments of 260 bp in length. The final library quantified (qPCR) by loading 1 µL of resuspended construct on an Agilent Technologies 2100 Bioanalyzer using a DNA-specific chip (Agilent DNA 1000). For cluster generation, the qualified and quantified libraries were first amplified within the flow cell on the cBot instrument (HiSeq® 4000 PE Cluster Kit, Illumina, San Diego, CA, USA).
For paired-end sequencing, the clustered flow cell was then loaded onto the HiSeq 4000 Sequencer (HiSeq® 4000 SBS Kit, Illumina, San Diego, CA, USA) with 100 bp which was the recommended read length. The library preparation and Illumina sequencing were performed by the Beijing Genomics Institute (BGI) (Hong Kong, China).

4.4. De Novo Transcriptome Assembly

In order to remove adaptors from the reads, low-quality reads, and reads in which unknown bases (N) comprised more than 5% of the read, raw Illumina paired-end reads were filtered using the SOAPnuke software (version: v.1.5.6, Beijing Genomics Institute, Shenzhen, China, https://github.com/BGI-flexlab/SOAPnuke,). Post-filtered reads were stored in the FASTQ format [53]. To obtain unigenes, clean reads were assembled using the Trinity software (version: v2.0.6, Trinity Software, Arlington, Tx, USA) [54]. The resulting sequences assembled using Trinity were referred to as transcripts. Gene family clustering was then carried out using TGICL (TIGR Gene Indices clustering tools) to obtain the final unigenes, which were classified to two categories: clusters and singletons. The former were labeled by the prefix ‘CL’, followed by the cluster ID. The latter were indicated by the prefix ‘unigene’.

4.5. Gene Annotation and Analysis

Identification and functional annotation of all unigene sequences were carried out in seven functional databases (e-value < 10−5): Nr, Nt, GO, COG, KEGG, Swiss-Prot, and Interpro databases. Blast (version: v2.2.23, NCBI, Bethesda, MD, USA, http://blast.ncbi.nlm.nih.gov/Blast.cgi) [55] was used to align the unigenes to NT, NR, COG, KEGG, and SwissProt to obtain annotations. Blast2GO (version: v2.5.0, BioBam, Valencia, Spain, https://www.blast2go.com) [56] used NR annotations to obtain GO annotations, and InterProScan5 (version: v5.11-51.0, EMBL-EBI, Hinxton, UK, https://code.google.com/p/interproscan/wiki/Introduction) to obtain InterPro annotations.

4.6. Differential Expression Analysis

Bowtie v2.2.5 was devoted to map the high-quality reads to the reference unigene sequences [57], and then calculate the gene expression levels, which were determined using RSEM (version: v1.2.12, http://deweylab.biostat.wisc.edu/RSEM) [58]. DEGs were detected based on a Poisson distribution using PossionDis, as described by Audic and Claverie [59]. The unigene expression level was calculated following the fragments per kilobase million (FPKM) formula. A false discovery rate (FDR) of 0.001 and a two-fold change were selected as the thresholds for significantly differential expression.

4.7. GO and KEGG Enrichment Analysis of Differentially Expressed Genes

DEGs were classified according to the official classification on the basis of the GO annotation results. Pathway functional enrichment was also carried out by the R-function phyper. The p value calculating formula in the hypergeometric test was as Equation (1):
P = 1 i = 0 m 1 ( M i ) ( N M n i ) ( N n )
FDR was calculated for each p value, and in general, the terms for which FDR did not exceed 0.001 were defined as significantly enriched.

4.8. Quantitative Real-Time PCR Validation

The expression of nine genes, which were singled out for the validation of the DGE data, was performed by qPCR. β-actin was used as a house-keeping gene [8]. cDNA synthesis was performed with 500 ng of DNase-treated RNA by using a PrimeScriptTM RT Reagent Kit (TaKaRa Bio, Kyoto, Japan). All qPCR reactions were carried out using SYBR Premix Ex TaqTM Perfect Real-Time Kits (TaKaRa Bio, Japan) with a QiagenRotor-Gene Q RT-PCR Detection System (Qiagen, Hilden, Germany). PCR primers, listed in Table S3, were designed using the Primer 5 software (version: v.5, PREMIER Biosoft, Palo Alto, CA, USA) based on transcriptome sequences. The reaction mixture consisted of 1 µL cDNA (50 ng), 1 µL of the forward and reverse primers (10 µM), and 6.25 µL of SYBR Premix Ex TaqTM. To ensure that the final volume of the reaction mixture was 12 µL, ultra-pure water was added. The following reaction conditions were maintained for extension: 94 °C for 2 min, followed by 40 cycles of 94 °C for 20 s, 58 °C for 30 s, and 72 °C for 40 s [8]. In order to eliminate the possibility of primer dimer formation or non-specific amplifications, a melting curve analysis was carried out after the amplification phase [8]. A standard curve was constructed from serial dilutions of the cDNA sample and drawn by plotting the natural log of the threshold cycle (Ct) against the number of molecules [8]. Standard curves for each gene were prepared in duplicate and triplicate to obtain a reliable measure of the amplification efficiency [8]. The amplification efficiencies were between 90% and 110%, and the correlation coefficients (R2) of all standard curves were >0.99. The relative expression ratios of the target genes were calculated using the method described by M.W. Pfaffl [60]. In all cases, PCR was carried out in triplicate. Statistical analysis was carried out using the statistical software SPSS 19.0 (version: 19.0, IBM Corp., Armonk, NY, USA, 2017). The differences were determined using the LSD test, with p-values < 0.05 indicating statistical significance. Values were expressed as the arithmetic mean ± standard deviation (SD).

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-6651/10/8/308/s1, Table S1: Summary of sequencing reads after filtering, and quality metrics of transcripts, Table S2: The differential expression unigenes (with higher than two-fold changes, and FDR≤ 0.001) between the control and the OA-treated groups. Table S3: All primers used in the validation analysis, The accession number for our raw dataset in the GEO database is: GSE116508.

Author Contributions

C.C. and S.C.P. conceived and designed the experiments; C.C. and S.S.G. performed the experiments; C.C. and J.W.J. analyzed the data; H.J.K., S.W.K., J.W.K. contributed with reagents/materials/analysis tools; C.C. wrote the manuscript.

Funding

This research was funded by the Basic Science Research Program through the National Research Foundation of Korea (NRF) [2017R1C1B2004616] and the supportive managing project of the Center for Companion Animals Research of the Korean government [PJ013877].

Conflicts of Interest

The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Cai, Y.; Pan, L.; Hu, F.; Jin, Q.; Liu, T. Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo [a] pyrene. Gene 2014, 551, 261–270. [Google Scholar] [CrossRef] [PubMed]
  2. Goldberg, E.D.; Bowen, V.T.; Farrington, J.W.; Harvey, G.; Martin, J.H.; Parker, P.L.; Risebrough, R.W.; Robertson, W.; Schneider, E.; Gamble, E. The mussel watch. Environ. Conserv. 1978, 5, 101–125. [Google Scholar] [CrossRef]
  3. Hu, F.; Pan, L.; Cai, Y.; Liu, T.; Jin, Q. Deep sequencing of the scallop Chlamys farreri transcriptome response to tetrabromobisphenol A (TBBPA) stress. Mar. Genom. 2015, 19, 31–38. [Google Scholar] [CrossRef] [PubMed]
  4. Liu, N.; Pan, L.; Wang, J.; Yang, H.; Liu, D. Application of the biomarker responses in scallop (Chlamys farreri) to assess metals and PAHs pollution in Jiaozhou Bay, China. Mar. Environ. Res. 2012, 80, 38–45. [Google Scholar] [CrossRef] [PubMed]
  5. De Jesús Romero-Geraldo, R.; 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]
  6. Zingone, A.; Enevoldsen, H.O. The diversity of harmful algal blooms: A challenge for science and management. Ocean. Coast. Manag. 2000, 43, 725–748. [Google Scholar] [CrossRef]
  7. Anderson, D.M. Approaches to monitoring, control and management of harmful algal blooms (HABs). Ocean. Coast. Manag. 2009, 52, 342–347. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. 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]
  9. Glibert, P.M.; Anderson, D.M.; Gentien, P.; Granéli, E.; Sellner, K.G. The global, complex phenomena of harmful algal blooms. Oceanography 2005, 18, 136–147. [Google Scholar] [CrossRef]
  10. Basti, L.; Hégaret, H.; Shumway, S.E. Harmful algal blooms and shellfish. In Harmful Algal Blooms: A Compendium Desk Reference; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2018; pp. 135–190. [Google Scholar]
  11. 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]
  12. Chi, C.; Giri, S.S.; Jun, J.W.; Kim, H.J.; Yun, S.; Kim, S.G.; Park, S.C. Marine Toxin Okadaic Acid Affects the Immune Function of Bay Scallop (Argopecten irradians). Molecules 2016, 21, 1108. [Google Scholar] [CrossRef] [PubMed]
  13. Espiña, B.; Louzao, M.; Cagide, E.; Alfonso, A.; Vieytes, M.R.; Yasumoto, T.; Botana, L.M. The methyl ester of okadaic acid is more potent than okadaic acid in disrupting the actin cytoskeleton and metabolism of primary cultured hepatocytes. Br. J. Pharmacol. 2010, 159, 337–344. [Google Scholar] [CrossRef] [PubMed]
  14. 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]
  15. Mouratidou, T.; Kaniou-Grigoriadou, I.; Samara, C.; Kouimtzis, T. Detection of the marine toxin okadaic acid in mussels during a diarrhetic shellfish poisoning (DSP) episode in Thermaikos Gulf, Greece, using biological, chemical and immunological methods. Sci. Total Environ. 2006, 366, 894–904. [Google Scholar] [CrossRef] [PubMed]
  16. Stonik, V.A.; Stonik, I.V. Toxins Produced by Marine Microorganisms: A Short Review. Mar. Freshw. Toxins 2016, 3–21. [Google Scholar] [CrossRef]
  17. 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 Part A 2015, 78, 814–824. [Google Scholar] [CrossRef] [PubMed]
  18. Mello, D.F.; Proença, L.A. d. O.; Barracco, M.A. Comparative study of various immuneparameters in three bivalve species during a natural bloom of Dinophysis acuminata in Santa Catarina Island, Brazil. Toxins 2010, 2, 1166–1178. [Google Scholar] [CrossRef] [PubMed]
  19. Burgos-Aceves, M.A.; Faggio, C. An approach to the study of the immunity functions of bivalve haemocytes: Physiology and molecular aspects. Fish Shellfish Immunol. 2017, 67, 513–517. [Google Scholar] [CrossRef] [PubMed]
  20. Prado-Alvarez, M.; Flórez-Barrós, F.; Méndez, J.; Fernandez-Tajes, J. Effect of okadaic acid on carpet shell clam (Ruditapes decussatus) haemocytes by in vitro exposure and harmful algal bloom simulation assays. Cell. Biol. Toxicol. 2013, 29, 189–197. [Google Scholar] [CrossRef] [PubMed]
  21. Fabioux, C.; Sulistiyani, Y.; Haberkorn, H.; Hégaret, H.; Amzil, Z.; Soudant, P. Exposure to toxic Alexandrium minutum activates the detoxifying and antioxidant systems in gills of the oyster Crassostrea gigas. Harmful Algae 2015, 48, 55–62. [Google Scholar] [CrossRef] [PubMed]
  22. Kim, C.S.; Lee, S.G.; Lee, C.K.; Kim, H.G.; Jung, J. Reactive oxygen species as causative agents in the ichthyotoxicity of the red tide dinoflagellate Cochlodinium polykrikoides. J. Plankton Res. 1999, 21, 2105–2115. [Google Scholar] [CrossRef]
  23. Flores, H.S.; Wikfors, G.H.; Dam, H.G. Reactive oxygen species are linked to the toxicity of the dinoflagellate Alexandrium spp. to protists. Aquat. Microb. Ecol. 2012, 66, 199–209. [Google Scholar] [CrossRef]
  24. Reuter, J.A.; Spacek, D.V.; Snyder, M.P. High-throughput sequencing technologies. Mol. Cell 2015, 58, 586–597. [Google Scholar] [CrossRef] [PubMed]
  25. Huang, Q.; Dong, S.; Fang, C.; Wu, X.; Ye, T.; Lin, Y. Deep sequencing-based transcriptome profiling analysis of Oryzias melastigma exposed to PFOS. Aquat. Toxicol. 2012, 120, 54–58. [Google Scholar] [CrossRef] [PubMed]
  26. Zhao, X.; Yu, H.; Kong, L.; Li, Q. Transcriptomic responses to salinity stress in the Pacific. oyster Crassostrea gigas. PLoS ONE 2012, 7, e46244. [Google Scholar] [CrossRef]
  27. Fu, X.; Sun, Y.; Wang, J.; Xing, Q.; Zou, J.; Li, R.; Wang, Z.; Wang, S.; Hu, X.; Zhang, L. Sequencing-based gene network analysis provides a core set of gene resource for understanding thermal adaptation in Zhikong scallop Chlamys farreri. Mol. Ecol. Resour. 2014, 14, 184–198. [Google Scholar] [CrossRef] [PubMed]
  28. Philipp, E.E.; Kraemer, L.; Melzner, F.; Poustka, A.J.; Thieme, S.; Findeisen, U.; Schreiber, S.; Rosenstiel, P. Massively parallel RNA sequencing identifies a complex immune gene repertoire in the lophotrochozoan Mytilus edulis. PLoS ONE 2012, 7, e33091. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Svensson, S.; Särngren, A.; Förlin, L. Mussel blood cells, resistant to the cytotoxic effects of okadaic acid, do not express cell membrane p-glycoprotein activity (multixenobiotic resistance). Aquat. Toxicol. 2003, 65, 27–37. [Google Scholar] [CrossRef]
  30. Ravindran, J.; Gupta, N.; Agrawal, M.; Bhaskar, A.B.; Rao, P.L. Modulation of ROS/MAPK signaling pathways by okadaic acid leads to cell death via, mitochondrial mediated caspase-dependent mechanism. Apoptosis 2011, 16, 145–161. [Google Scholar] [CrossRef] [PubMed]
  31. Von Zezschwitz, C.; Vorwerk, H.; Tergau, F.; Steinfelder, H.J. Apoptosis induction by inhibitors of Ser/Thr phosphatases 1 and 2A is associated with transglutaminase activation in two different human epithelial tumour lines. FEBS Lett. 1997, 413, 147–151. [Google Scholar] [CrossRef] [Green Version]
  32. Soliño, L.; Sureda, F.X.; Diogène, J. Evaluation of okadaic acid, dinophysistoxin-1 and dinophysistoxin-2 toxicity on Neuro-2a, NG108-15 and MCF-7 cell lines. Toxicol. Vitro 2015, 29, 59–62. [Google Scholar] [CrossRef] [PubMed]
  33. Marioni, J.C.; Mason, C.E.; Mane, S.M.; Stephens, M.; Gilad, Y. RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008, 18, 1509–1517. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Simon, P. Q-Gene: Processing quantitative real-time RT–PCR data. Bioinformatics 2003, 19, 1439–1440. [Google Scholar] [CrossRef] [PubMed]
  35. Clayton, M.E.; Steinmann, R.; Fent, K. Different expression patterns of heat shock proteins hsp 60 and hsp 70 in zebra mussels (Dreissena polymorpha) exposed to copper and tributyltin. Aquat. Toxicol. 2000, 47, 213–226. [Google Scholar] [CrossRef]
  36. Boutet, I.; Tanguy, A.; Moraga, D. Response of the Pacific oyster Crassostrea gigas to hydrocarbon contamination under experimental conditions. Gene 2004, 329, 147–157. [Google Scholar] [CrossRef] [PubMed]
  37. Faggio, C.; Pagano, M.; Alampi, R.; Vazzana, I.; Felice, M.R. Cytotoxicity, haemolymphatic parameters, and oxidative stress following exposure to sub-lethal concentrations of quaternium-15 in Mytilus galloprovincialis. Aquat. Toxicol. 2016, 180, 258–265. [Google Scholar] [CrossRef] [PubMed]
  38. Piano, A.; Valbonesi, P.; Fabbri, E. Expression of cytoprotective proteins, heat shock protein 70 and metallothioneins, in tissues of Ostrea edulis exposed to heat andheavy metals. Cell. Stress Chaperones 2004, 9, 134–142. [Google Scholar] [CrossRef] [PubMed]
  39. Anzenbacher, P.; Anzenbacherová, E. Cytochromes P450 and metabolism of xenobiotics. Cell. Mol. Life Sci. 2001, 58, 737–747. [Google Scholar] [CrossRef] [PubMed]
  40. Liu, D.; Pan, L.; Cai, Y.; Li, Z.; Miao, J.J. Response of detoxification gene mRNA expression and selection of molecular biomarkers in the clam Ruditapes philippinarum exposed to benzo[a]pyrene. Environ. Pollut. 2014, 189, 1–8. [Google Scholar] [CrossRef] [PubMed]
  41. 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]
  42. Huang, L.; Wang, J.; Chen, W.C.; Li, H.Y.; Liu, J.S.; Jiang, T.; Yang, W.D. P-glycoprotein expression in Perna viridis after exposure to Prorocentrum lima, a dinoflagellate producing DSP toxins. Fish Shellfish Immunol. 2014, 39, 254–262. [Google Scholar] [CrossRef] [PubMed]
  43. Donaghy, L.; Hong, H.K.; Jauzein, C.; Choi, K.S. The known and unknown sources of reactive oxygen and nitrogen species in haemocytes of marine bivalve molluscs. Fish Shellfish Immunol. 2015, 42, 91–97. [Google Scholar] [CrossRef] [PubMed]
  44. Bedard, K.; Krause, K.H. The NOX family of ROS-generating NADPH oxidases: Physiology and pathophysiology. Physiol. Rev. 2007, 87, 245–313. [Google Scholar] [CrossRef] [PubMed]
  45. Liu, D.; Pan, L.; Li, Z.; Cai, Y.; Miao, J. Metabolites analysis, metabolic enzyme activities and bioaccumulation in the clam Ruditapes philippinarum exposed to benzo [a] pyrene. Ecotoxicol. Environ. Saf. 2014, 107, 251–259. [Google Scholar] [CrossRef] [PubMed]
  46. 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]
  47. Jin, Q.; Pan, L.; Liu, T.; Hu, F. RNA-seq based on transcriptome reveals differ genetic expressing in Chlamys farreri exposed to carcinogen PAHs. Environ. Toxicol. Pharmacol. 2015, 39, 313–320. [Google Scholar] [CrossRef] [PubMed]
  48. Regoli, F.; Benedetti, M.; Giuliani, M.E. Antioxidant defenses and acquisition of tolerance to chemical stress. Toler. Environ. Contam. 2011, 153–173. [Google Scholar]
  49. Regoli, F.; Giuliani, M.E.; Benedetti, M.; Arukwe, A. Molecular and biochemical biomarkers in environmental monitoring: A comparison of biotransformation and antioxidant defense systems in multiple tissues. Aquat. Toxicol. 2011, 105, 56–66. [Google Scholar] [CrossRef] [PubMed]
  50. 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]
  51. Chen, M.Y.; Yang, H.S.; Delaporte, M.; Zhao, S.J.; Xing, K. Immune responses of the scallop Chlamys farreri after air exposure to different temperatures. J. Exp. Mar. Biol. Ecol. 2007, 345, 52–60. [Google Scholar] [CrossRef]
  52. Zhou, J.; Xiong, Q.; Chen, H.; Yang, C.; Fan, Y. Identification of the spinal expression profile of non-coding RNAs involved in neuropathic pain following spared nerve injury by sequence analysis. Front. Mol. Neurosci. 2017, 10, 91. [Google Scholar] [CrossRef] [PubMed]
  53. Cock, P.J.; Fields, C.J.; Goto, N.; Heuer, M.L.; Rice, P.M. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2010, 38, 1767–1771. [Google Scholar] [CrossRef] [PubMed]
  54. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Altschul, S.F.; Gish, W.; Miller, W.; Myers, E.W.; Lipman, D.J. Basic local alignment search tool. J. Mol. Biol. 1990, 215, 403–410. [Google Scholar] [CrossRef]
  56. 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]
  57. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef] [PubMed]
  59. Audic, S.; Claverie, J.-M. The significance of digital gene expression profiles. Genome Res. 1997, 7, 986–995. [Google Scholar] [CrossRef] [PubMed]
  60. 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]
Figure 1. Distribution of all-unigenes in the bay scallop transcriptome.
Figure 1. Distribution of all-unigenes in the bay scallop transcriptome.
Toxins 10 00308 g001
Figure 2. Annotated species and their distribution.
Figure 2. Annotated species and their distribution.
Toxins 10 00308 g002
Figure 3. COG functional classification of All-unigenes.
Figure 3. COG functional classification of All-unigenes.
Toxins 10 00308 g003
Figure 4. Gene transcription profile of the control (CN) and the OA-exposed group (OA) libraries. Blue points represent downregulated genes. Red points represent upregulated genes. Black points represent non-differential expression genes.
Figure 4. Gene transcription profile of the control (CN) and the OA-exposed group (OA) libraries. Blue points represent downregulated genes. Red points represent upregulated genes. Black points represent non-differential expression genes.
Toxins 10 00308 g004
Figure 5. M (log ratio) and A (mean average) (MA) plot of DEGs of the control (CN) and the OA-exposed group (OA) libraries. X-axis represent value A (log2 mean expression level). Y-axis represents value M (log2 transformed fold change). Red points represent upregulated DEG. Blue points represent downregulated DEG. Black points represent non-DEGs.
Figure 5. M (log ratio) and A (mean average) (MA) plot of DEGs of the control (CN) and the OA-exposed group (OA) libraries. X-axis represent value A (log2 mean expression level). Y-axis represents value M (log2 transformed fold change). Red points represent upregulated DEG. Blue points represent downregulated DEG. Black points represent non-DEGs.
Toxins 10 00308 g005
Figure 6. GO classification of differentially expressed gene (DEGs). X-axis represent the GO term.
Figure 6. GO classification of differentially expressed gene (DEGs). X-axis represent the GO term.
Toxins 10 00308 g006
Figure 7. Pathway classification of DEGs. The X-axis shows the number of DEGs. The Y-axis shows the pathway name.
Figure 7. Pathway classification of DEGs. The X-axis shows the number of DEGs. The Y-axis shows the pathway name.
Toxins 10 00308 g007
Figure 8. Enrichment of DEGs and pathways. The X-axis indicates enrichment factor and the Y-axis indicates the pathway name. Coloring indicates the q value (high: white, low: blue), the lower q value indicates the more significant enrichment. The point size indicates the DEG number (more: big, less: small).
Figure 8. Enrichment of DEGs and pathways. The X-axis indicates enrichment factor and the Y-axis indicates the pathway name. Coloring indicates the q value (high: white, low: blue), the lower q value indicates the more significant enrichment. The point size indicates the DEG number (more: big, less: small).
Toxins 10 00308 g008
Figure 9. Results of the qPCR analysis. The y-axis represents the gene expressed log2 (fold change) and the x-axis is the gene name. SOD2 = superoxide dismutase 2, Cu/Zn SOD = copper and zinc superoxide dismutase, GST = glutathione-S-transferase, ACP = acid phosphatase, NOX-3 = NADPH oxidase 3, HSP70 = heat shock protein 70, CYP3A24 = cytochrome P450 3A24, CREB = cyclic adenosine monophosphate responsive element binding protein.
Figure 9. Results of the qPCR analysis. The y-axis represents the gene expressed log2 (fold change) and the x-axis is the gene name. SOD2 = superoxide dismutase 2, Cu/Zn SOD = copper and zinc superoxide dismutase, GST = glutathione-S-transferase, ACP = acid phosphatase, NOX-3 = NADPH oxidase 3, HSP70 = heat shock protein 70, CYP3A24 = cytochrome P450 3A24, CREB = cyclic adenosine monophosphate responsive element binding protein.
Toxins 10 00308 g009
Table 1. Quality metrics of unigenes.
Table 1. Quality metrics of unigenes.
SampleTotal NumberTotal LengthMean LengthN50N70N90GC(%)
Control51,46541,105,722798141170430239.48
OA-treated49,45343,129,157872164680331839.63
All-unigene55,87653,465,429956184096034539.42
N50: a weighted median statistic within which 50% of the Total Length is contained in unigenes greater than or equal to this value. GC (%): the percentage of G and C bases in all unigenes.
Table 2. Detoxification and immune-related differentially expressed genes (DEGs) in bay scallop gills regulated after up to 48 h exposure to 500 nM OA.
Table 2. Detoxification and immune-related differentially expressed genes (DEGs) in bay scallop gills regulated after up to 48 h exposure to 500 nM OA.
FunctionTranscriptLog2 (Fold Change) (RNAseq)Regulation
Immune systemC-type lectin superfamily 17 member A−4.255Down
C-type lectin domain family 4 member E−3.507Down
Complement C1q tumor necrosis factor-related protein 2−4.791Down
Fibrinogen C domain-containing protein 1−2.100Down
Toll-like receptor 42.880Up
Toll-like receptor 131.347Up
Acid phosphatase2.238Up
NADPH oxidase 32.493Up
DetoxificationATP-binding cassette, subfamily C, member 11.773Up
ATP-binding cassette sub-family B member 101.165Up
ATP-binding cassette, sub-family C member 51.280Up
Cyclic AMP-responsive element-binding protein1.953Up
Nuclear factor erythroid 2-related factor 21.231Up
NADPH2:quinone reductase1.677Up
Cytochrome P450 3A801.207Up
Cytochrome P450 3A641.783Up
Cytochrome P450 1A5−1.686Down
Cytochrome P450 3A24−2.315Down
Superoxide dismutase Cu-Zn family1.139Up
Superoxide dismutase 2−1.126Down
Glutathione S-transferase 1−1.552Down
Glutathione S-transferase 2−2.511Down
Glutathione S-transferase omega−1.775Down
Glutathione S-transferase theta-1−1.254Down
Glutathione S-transferase A−1.218Down
Glutathione S-transferase kappa−2.356Down

Share and Cite

MDPI and ACS Style

Chi, C.; Giri, S.S.; Jun, J.W.; Kim, S.W.; Kim, H.J.; Kang, J.W.; Park, S.C. Detoxification- and Immune-Related Transcriptomic Analysis of Gills from Bay Scallops (Argopecten irradians) in Response to Algal Toxin Okadaic Acid. Toxins 2018, 10, 308. https://doi.org/10.3390/toxins10080308

AMA Style

Chi C, Giri SS, Jun JW, Kim SW, Kim HJ, Kang JW, Park SC. Detoxification- and Immune-Related Transcriptomic Analysis of Gills from Bay Scallops (Argopecten irradians) in Response to Algal Toxin Okadaic Acid. Toxins. 2018; 10(8):308. https://doi.org/10.3390/toxins10080308

Chicago/Turabian Style

Chi, Cheng, Sib Sankar Giri, Jin Woo Jun, Sang Wha Kim, Hyoun Joong Kim, Jeong Woo Kang, and Se Chang Park. 2018. "Detoxification- and Immune-Related Transcriptomic Analysis of Gills from Bay Scallops (Argopecten irradians) in Response to Algal Toxin Okadaic Acid" Toxins 10, no. 8: 308. https://doi.org/10.3390/toxins10080308

APA Style

Chi, C., Giri, S. S., Jun, J. W., Kim, S. W., Kim, H. J., Kang, J. W., & Park, S. C. (2018). Detoxification- and Immune-Related Transcriptomic Analysis of Gills from Bay Scallops (Argopecten irradians) in Response to Algal Toxin Okadaic Acid. Toxins, 10(8), 308. https://doi.org/10.3390/toxins10080308

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