2 Shandong Marine Resource and Environment Research Institute, Shandong Provincial Key Laboratory of Restoration for Marine Ecology, Yantai 264006, China;
3 Weihai Hongrun Marine Technology Co. Ltd., Weihai 264200, China
Polybrominated diphenyl ethers (PBDEs) comprise a series of aromatic compounds containing bromine atoms, where there are 209 congeners according to the degree of bromination. They are used widely as brominated flame retardants (BFRs) in electronics, textiles, plastics, resins, and building materials, and they are readily leaked into the environment (Law et al., 2006). Thus, PBDEs have been detected in many types of environmental media and organisms. They can enter the ocean via atmospheric deposition and surface runoff, where they accumulate in sediments and marine organisms. The PBDEs levels are as high as 1 800 ng/g in sediments and 720 ng/g lipid in mussels around Laizhou Bay in China, which is surrounded by BFR manufacturers. 2, 2′, 4, 4′-Tetrabromodiphenyl ether (BDE-47), a type of tetra-PBDE, is most common in these samples (Jin et al., 2008). Similar results have been obtained in human tissues and wild animal samples (Schecter et al., 2007; Barja-Fernández et al., 2013). PBDEs are considered to be an increasing pollution problem in the coastal waters of China.
The toxicological effects of PBDEs on aquatic organisms are a major concern due to their persistent and bioaccumulative properties. PBDEs can cause immune, developmental, and reproduction dysfunctions, as well as affect the endocrine systems offish, mollusks, copepods, and other marine organisms (Breitholtz and Wollenberger, 2003; Yu et al., 2015). However, the toxicological mechanism responsible for the effects of PBDEs remains unclear.
Bivalves are used throughout the world as sentinel species to monitor coastal environments due to their ease of sampling, sessile habit, high filtration rate, and great efficiency at bioaccumulating pollutants (Santovito et al., 2005, 2015; Irato et al., 2007). Many studies have used bivalves to determine the impact of PBDEs pollution on marine organisms and the responses of marine organisms to PBDE pollution. In the mussel Mytilus edilus, exposure to BDE-47 for 3 weeks increased the frequency of micronuclei, as well as inducing bi-nucleated, fragmented apoptotic cells, and nuclear buds, damaging ovarian follicles and ovocytes in females, inducing spawning in males, and decreasing the phospho-protein levels in both sexes (Aarab et al., 2006; Baršienė et al., 2006). In another study, metabolomic analysis of the mussel Mytilus galloprovincialis indicated that BDE-47 mainly disrupted energy metabolism in males but osmotic regulation and energy metabolism in females, where proteomic analysis demonstrated that BDE-47 disturbed protein homeostasis in males and proteolysis in females (Ji et al., 2013). mRNA expression occurs earlier than protein expression. Thus, in order to identify changes in gene transcription caused by BDE-47, suppression subtractive hybridization (SSH) was used successfully in the manila clam Ruditapes philippinarum (Miao et al., 2014). However, the genes that are not represented by specific probes and with low expression levels cannot be detected by SSH. Thus, the ability of SSH to identify novel transcripts and elucidate the mechanism of toxicity is limited at the mRNA level. Due to its obvious advantages in terms of the sequencing depth, sequencing accuracy, detection sensitivity, genome coverage, and molecular pathway analysis, highthroughput transcriptome sequencing has been used widely to detect transcripts in aquatic organisms, thereby obtaining a clearer understanding of the mechanism of xenobiotics action and to provide molecular biomarkers for pollution monitoring. Recently, the transcriptomic profiles of zebrafish larvae exposed to BDE-209 were determined using high-throughput transcriptome sequencing (Chen et al., 2016). The estrogenic effects of BDE-209 were demonstrated and numerous genes that responded to BDE-209 were explored.
The clam Mactra veneriformis, a type of infaunal suspension-feeding bivalve, is one of the major economic mollusks in China. In our previous study, we demonstrated that M. veneriformis is a good model for toxicological research and a good indicator of marine pollution (Fang et al., 2010, 2013). Furthermore, several studies have reported that high concentrations of PBDEs accumulate in this clam. Therefore, we employed high-throughput transcriptomic profiling of M. veneriformisafter BDE-47 exposure to improve current understanding of the toxic effects and mechanism of action for BDE-47 in aquatic organisms.2 MATERIAL AND METHOD 2.1 Animals and BDE-47 exposure
Adult M. veneriformis clams (average shell length: 3.46±0.24 cm; height: 3.06±0.21 cm; width: 2.38±0.32 cm; mass: 12.45±0.82 g) were collected from the coast at Dongying and acclimated in aerated tanks at 16℃ for two weeks before the exposure experiment.
After acclimatization, the clams were transferred randomly to glass aquaria filled with 30 L seawater and separated into two groups (control and treatment). In the treatment group, the clams were exposed to BDE-47 for 3 days at a nominal concentration of 5 μg/L, which is the sublethal BDE-47 concentration used in bivalves. It has been reported that nuclear injuries can be detected in mussels exposed to 5 μg/L BDE-47 for 3 weeks (Baršienė et al., 2006). Moreover, differentially expressed genes (DEGs) were identified in the clam R. philippinarum after exposure to 5 μg/L BDE-47 for 15 days (Miao et al., 2014). BDE-47 (CAS No. 5436-43-1, purity 98.5%, AccuStandard Inc.) was first dissolved in dimethyl sulfoxide (DMSO) and added to the aquarium water to obtain a final DMSO concentration of 0.05% v/v. In the control group, DMSO was added directly to the aquarium water at a final concentration of 0.05% v/v. Three replicates were performed for the control and treatment, where each replicate comprised 30 clams. During the experiment, the water was aerated constantly and changed daily. The clams were fed 2 h before the water was changed, and the appropriate concentrations of BDE-47 and DMSO were added after changing the water.
The digestive gland of mollusks is the main tissue involved with detoxification and the elimination of xenobiotics (Moore and Allen, 2002). The digestive gland was used successfully in our previous studies to investigate the molecular response of M. veneriformis to environmental pollution (Fang et al., 2010, 2013). Thus, the digestive gland was sampled after 3 days exposure in this study. The digestive glands were dissected from 30 clams in both the control and treatment groups. For each clam, half of the digestive gland was used for high-throughput transcriptomic sequencing and the other half was used for real-time PCR validation. The samples were frozen in liquid nitrogen immediately after dissection and stored at -80℃ until subsequent analysis. The samples used for high-throughput sequencing were pooled together for RNA extraction. The samples used for real-time PCR analysis were assigned to five replicates (digestive glands from six clams per replicate) for RNA extraction.2.2 Total RNA extraction and Illumina sequencing
TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract the total RNA from the control and treatment samples according to the manufacturer's protocol. After RNase-free DNase (Promega, Madison, Wisconsin, USA) treatment, the purity of the total RNA was determined using an ultraviolet spectrophotometer (Thermo Scientific NanoDrop 8000). The RNA concentration was measured using a Qubit®RNA Assay Kit in Qubit®2.0 Flurometer (Life Technologies, Carlsbad, California, USA). RNA integrity was checked using an RNA Nano 6000 Assay Kit for the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, California, USA). cDNA library construction was performed by the Beijing Genomics Institute (Shenzhen, China) using aNEBNext® Ultra™ RNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, Massachusetts, USA). Briefly, mRNA was separated from the total RNA using poly-T oligo-attached magnetic beads and then disrupted into small pieces. The short mRNA was transcribed into first-strand cDNA using reverse transcriptase and random hexamer primers. After adding dNTPs, RNaseH, and DNA polymerase Ⅰ, the second-strand cDNA was synthesized using the first-strand cDNA as the template. Next, the double-stranded cDNA was purified with a QiaQuick PCR extraction kit (Qiagen, Hilden, Germany), before repairing the cohesive ends, and ligating with sequencing adaptors. After screening the target fragments with agarose gel electrophoresis, PCR amplification was performed using the recovered fragments as the templates. The quality and output of the cDNA library were tested using an Agilent 2100 Bioanalyzer and an ABI StepOne Plus Real-time PCR System. Library sequencing was then performed using the IlluminaHiseq 2000 platform.2.3 Bioinformatics analyses
The raw reads were filtered to obtain clean reads by discarding reads containing adaptors or of low quality. De novo assembly of the transcriptomes was conducted using Trinity (release-20120608), a short read assembly program (Grabherr et al., 2011). BLASTx was employed to annotate the unigenes obtained by using the NR (non-redundant protein database, NCBI), NT (non-redundant nucleotide database, NCBI), Swiss-Prot, KEGG (Kyoto Encyclopedia of Genes and Genomes), and COG (clusters of orthologous groups) databases. The sequence orientation of unigenes was determined based on the best results. If the results from different databases conflicted, the priority order was NR, Swiss-Prot, KEGG, and then COG. The software tool ESTScan (Iseli et al., 1999) was used to determine the sequence directionality if a unigene did not align with any of the databases mentioned above. Functional annotations were processed using the Blast2Go program (Conesa et al., 2005) and then classified with WEGO software (Ye et al., 2006). The metabolic and signaling pathways of unigenes were also annotated by the Blastall program based on the KEGG annotations (Kanehisa et al., 2008).
Unigene expression was calculated and normalized by the fragments per kb per million fragments (FPKM) method (Mortazavi et al., 2008). The differentially expressed unigenes in the control group (CG) and BDE-47 group (BDEG) were analyzed usingmodified Audic's method (Audic and Claverie, 1997), as described in a previous study (Liu et al., 2011). DEGs were identified based on a false discovery rate ≤0.001 and absolute value of log2Ratio ≥1.
In order to determine the main functional classifications of DEGs, we performed GO functional classification for the DEGs as well as GO and pathway enrichment analysis. After mapping the DEGs to GO terms or pathways, we calculated the gene numbers for every term or pathway and performed a hypergeometric test, where the significantly enriched GO terms and pathways were identified based on a Bonferroni corrected P value ≤0.05 and Q value ≤0.05, respectively.2.4 Quantitative real-time reverse-transcription (RT)-PCR validation
We employed 26 genes related to detoxification, apoptosis, signal transduction, immune defense, disease, and antioxidant defense to validate the Illumina sequencing results by real-time PCR. Primer 5.0 was used to design the primers for these genes based on the sequences obtained by Illumina sequencing. The specific primers are shown in Table 1. β-Actin was employed as an internal control. The samples employed for RT-PCR analysis (as described in Section 2.1) were used for RNA extraction according to the same method employed for Illumina sequencing preparation. RNA was reverse transcribed into cDNA using M-MLV Reverse Transcriptase (Promega) and oligod (T) primers after digestion with DNase (Promega). DEPC-water was used as a negative control to replace the cDNA template. SYBR Green® real-time PCR assay (SYBR Prime Script TM RT-PCR Kit Ⅱ; TaKaRa, Dalian, China) was used to determine the mRNA expression levels with an Eppendorf Mastercycler® eprealplex S (Eppendorf, Hamburg, Germany). The amplification volume was 20 μL, which comprised 10 μL of SYBR Green Master Mix (TaKaRa), 0.4 μL of each forward and reverse primer (10 mmol/L), 1 μL of 1:10 diluted cDNA, and 8.2 μL of DEPC-treated water. The reaction conditions were as follows: 95℃ for 5 min followed by 40 cycles at 95℃ for 15 s and at 60℃ for 30 s. Melting-curve analysis of the amplification products was performed at the end of each PCR reaction to confirm that only one PCR product was amplified and detected. The PCR products were also confirmed by electrophoresis and some were sequenced. The 2-ΔΔCt method was used to analyze the relative expression levels. For both the CG and BDEG groups, five biological pools (n=5) were analyzed, where each pool was a mixture of six different individual clams (n=6). All of the data were expressed as the mean±standard deviation (n=5). We employed one-way analysis of variance followed by Tukey's test using the SPSS statistical package (version 18.0; SPSS Inc., Chicago, IL, USA) to test for significant differences between the CG and BDEG. Significant differences were accepted when P < 0.05.3 RESULT 3.1 Sequencing and read assembly
After filtering the raw reads, 50 142 012 nt and 50 152 900 nt clean reads were generated from CG and BDEG by Illumina RNA-seq deep sequencing, respectively, and they have been submitted to NCBI (accession No. GSE88918). The de novo assembly yielded 247 603 and 301 96 contigs with mean lengths of 282 bp and 280 bp in CG and BDEG, respectively (Table 2). The contigs of the two groups were resolved into 127 648 and 154 225 unigenes with average lengths of 530 bp and 529 bp, respectively. In CG, 33 092 were clusters and 94 556 were singletons, and in BDEG, 44 730 were clusters and 109 495 were singletons. In addition, 124 985 unigenes with a mean length of 733 bp were obtained from the combined clean reads from CG and BDEG (CG & BDEG) using the same assembly strategy. There were 47 425 clusters and 7 560 singletons in CG & BDEG.
Similar length distributions were found for the CG unigenes and BDEG unigenes, where sequences comprising 200–800 bp accounted for the highest proportion, followed by sequences > 3 000 bp (Fig. 1). Furthermore, sequences of 200–300 bp were not found in the CG & BDEG unigenes because of the improved length distribution (Fig. 1c). The minimum length of the CG & BDEG unigenes was 300 bp (Fig. 1c), whereas it was 200 bp for the CG unigenes or BDEG unigenes (Fig. 1a, b). In general, compared with the CG unigenes or BDEG unigenes, the quality and the identity of the known proteins determined in NR was better for the CG & BDEG unigenes.3.2 Structural and functional annotations
We annotated 34 255 unigenes among the CG & BDEG unigenes using BLASTX and ESTscan with the NR, NT, Swiss-Prot, KEGG, and COG databases. Most of the unigenes (32 176 of 34 255) were annotated using the NR database. Compared with the NR database, the distributions of the E-value and similarity showed that 31.13% (E-value of 0 to 1e-60) and 24.33% (100%–60% similarity), respectively, of the sequences possessed high homology (Fig. 2a, b). On a species basis, the highest proportion of matching sequences in the NR database was derived from Crassostrea gigas (54.6%), followed by Capitella sp. Ⅰ (7.95%, Fig. 2c).
COG analysis was used to annotate the unigenes in order to examine the integrity of the transcriptome library and the effectiveness of the annotation process. We categorized 24 401 unigenes into 25 functional COG clusters (Fig. 3). The five largest categories were: (1) general function predictions only (19.3%); (2) recombination and repair (8.8%); (3) translation, ribosomal structure, and biogenesis (7.6%); (4) transcription (7.0%); and (5) cell cycle control, cell division, and chromosome partitioning (6.5%). The M. veneriformis unigenes were categorized into three main GO categories: biological process (55 230; 53.2%), cellular component (33 434; 32.2%), and molecular function (15 151; 14.6%). These GO terms were further subdivided into 56 sub-categories and they corresponded to the categories determined by COG analysis (Fig. 4).
To analyze the potential involvement of the consensus sequences in cellular metabolic pathways, the annotated unigenes were mapped to the KEGG database and 21 749 unigenes were assigned to 259 KEGG pathways. Metabolic pathways accounted for the largest number of unigenes (2 555; 11.8%, ko01100), followed by focal adhesion (1 135; 5.2%, ko04510), regulation of actin cytoskeleton (1 061; 4.9%, ko04810), pathways in cancer (1.035; 4.8%, ko05200), vascular smooth muscle contraction (894; 4.1%, ko04270), and tight junction (711 members; 3.3%, ko04530) (Table 3).3.3 DEGs in CG and BDEG
We explored 17 625 DEGs using arigorous algorithm and the FPKM method as described in Section 2.3. We found that 10 028 genes were significantly upregulated and 7 597 genes were significantly downregulated after BDE-47 exposure (Fig. 5). The key DEGs after exposure are shown in Table 4.3.4 GO ontology analysis of DEGs
To understand their functions, all of the DEGs were mapped onto terms in the GO database. GO ontology analysis assigned GO categories to 9 543 DEGs (Fig. 6), where "cellular process" was the dominant term in the ontology for "biological process". Stress-related processes such as "metabolic process", "biological regulation", "regulation of biological process", and "response to stimulus" accounted for the highest proportion. The dominant terms in the ontology for "cellular component" were "cell" and "cell part". In the ontology for "molecular function", "binding" was dominant, and "catalytic activity" and "transporter activity" were also abundant. Moreover, the significantly enriched GO terms were screened based on a corrected P < 0.05. One and seven significantly enriched GO terms were classified as "biological process" and "molecular function", respectively (Table 5). No significantly enriched GO term was found for "cellular component".3.5 Pathway enrichment analysis of DEGs
Enrichment analysis of DEGs was performed to further explore the biological pathways involved in the response to BDE-47. We mapped 2 795 DEGs into 255 metabolic pathways and signal transduction pathways, where 63 pathways were significantly enriched. The top 20 significantly enriched pathways are listed in Table 6. As expected, some pathways that respond to environmental stress, such as cell adhesion molecules, the NF-kappa B signaling pathway, ATP binding cassette (ABC) transporters, and cytokinecytokine receptor interaction, were significantly enriched. Pathways related to the immune system were also detected, such as the T cell receptor signaling pathway (ko04660), complement and coagulation cascades (ko04610), and Fc gamma R-mediated phagocytosis (ko04666). Interestingly, we found some DEGs that are implicated in many disease pathways, such as primary immunodeficiency (ko05340), prion diseases (ko05020), toxoplasmosis (ko05145), and systemic lupus erythematosus (ko05322).3.6 Verification of DEGs by RT-PCR
We used RT-PCR to confirm the relative mRNA expression levels of 26 DEGs identified by highthroughput sequencing, i.e., nineteen upregulated and seven downregulated genes. These DEGs were related to detoxification of xenobiotics, apoptosis, signal transduction, immune defense, disease, and antioxidant defense. The fold change differences between the RTPCR and high-throughput sequencing results are shown in Fig. 7. Except for ABC transporter G family member 22, all of the verified genes exhibited similar trends in their relative mRNA expression levels according to the two methods. The mRNA expression level of ABC transporter G family member 22 in BDE-47 treated group was significantly upregulated according to high-throughput sequencing but it did not change significantly when measured by RT-PCR. The correlation coefficient between the two series of data measured by the two methods was 0.987 (P < 0.000 1).4 DISCUSSION
Due to growing evidence for their toxicity and persistence in the environment and organisms, the production of penta-and octa-BDEs was listed in the "Stockholm Convention" and it ceased in Europe in 1998, and in USA and Canada in 2004 (Stockholm Convention, 2015). Deca-BDE is still used in many countries, especially in Asia (Jin et al., 2008). However, pollution by PBDEs is ubiquitous, where they have been detected in air, soil, sediment, aquatic organisms, plants, wildlife, human tissues, and even the polar bear. BDE-47 is dominant in these samples and it accounts for as much as 70% of the total PBDEs present in biota (Wollenberger et al., 2005). Thus, BDE-47 is attracting increasing attention as a pollutant.
M. veneriformis is a newly emerging model for marine eco-toxicological studies and transcriptomic information can be obtained for molecular analyses using this bivalve model. The powerful de novo assembly approach based on Illumina RNA-seqdeep sequencing is increasingly applied to obtain sequence information. Thus, in this study, we conducted deep sequencing-based transcriptome profiling analysis using M. veneriformis to elucidate the mechanisms responsible for the toxic effects of BDE-47 in bivalves. We obtained 127 648, 154 225, and 124 985 unigenes by assembling the CG reads, BDEG reads, and CG & BDEG reads, separately. This massive amount of sequencing data yielded a DEG catalog, which can be used for toxicology research at the molecular level.
RNA-seq transcriptional analysis identified DEGs in the CG and BDEG libraries, where 10 028 genes were significantly upregulated and 7 597 genes were downregulated. Using SSH, Miao et al. (2014) only identified 75 DEGs in R. philippinarum after exposure to BDE-47. Compared with the traditional method, RNA-seq has many advantages for detecting novel transcripts. In addition, the reliability and accuracy of the DEG data were verified by our RT-PCR results. GO and KEGG analysis were performed to clarify the functions of the DEGs, which were mainly related to metabolism and detoxification of xenobiotics, responses to stimuli, apoptosis, signal transduction, immune defense, and antioxidant defense.4.1 Genes associated with metabolism and detoxification of xenobiotics
Three main classes (phases) of proteins are responsible for the metabolism and detoxification of xenobioticsin organisms. The cytochrome p450 superfamily of heme-containing monooxygenases is typical phase Ⅰ enzymes that are responsible for organic pollutant biotransformation. Several MVCYP 450 genes (CYP1A1, CYP3A56, CYP2C8, CYP10, CYP2J6, CYP4F2, CYP3A13, and CYP3A31) were upregulated whereas CYP2F2 and CYP4 were downregulated in response to BDE-47. Given the important role of the metabolism of xenobiotics, CYP1A has been used as a biomarker in assessments of aquatic contamination (Brammell et al., 2010). Transcription of the CYP1A1 gene was upregulated 10.75-fold by BDE-47 in the present study. Giannetti et al. (2012) and Fossi et al. (2008) reported that expression of the CYP1A protein exhibited dosedependent induction in skinbiopsy slices from Caretta caretta and the Mediterranean cetacean Tursiops truncates when treated with polyaromatic hydrocarbons (PAHs), organochlorines, and PBDEs. Higher induction was found after treatment with PBDEs, thereby indicating the high potential of these chemicals for interfere with bioactivating enzyme systems. A field study showed that the mRNA expression level of CYP1A was elevated significantly in redeye mullet Liza haematocheila collected from a polluted area compared with a reference site, where the level correlated with the pollutant contents of their tissues, particularly with those of polychlorinated biphenyls (PCBs) and PBDEs (An et al., 2011). Nevertheless, it was reported that PBDEs are not capable of inducing CYP1A1 activity in cynomolgus monkey Macaca fascicularis hepatocytes (Peters et al., 2006). The different responses of species to the various congeners of PBDEs used in these studies may have yielded inconsistent results. In mammals, the CYP1 gene expression level is regulated by the aryl hydrocarbon receptor (Ah receptor), which is a cytosolic transcriptional factor that has a high affinity for planar compounds, such as PAHs and PCBs. The Ah receptor forms a complex with Ah receptor nuclear translocator (Arnt), which can bind to xenobioticresponsive elements in the promoters of CYP1 genes (Peters et al., 2006). PBDEs share structural similarity with PCBs. In the present study, we identified one Ah receptor gene and two Arnt genes in M. veneriformis. No obvious changes in the MVArnt gene expression levels were found after BDE-47 exposure. However, the expression level of the Ah receptor gene was about 73% lower in the BDEG group than the CG group. It would be interesting to clarify whether the induction of CYP1A1 by BDE-47 in M. veneriformis was a result of Ah receptor binding, which will be addressed in our future research.
We found that the mRNA expression levels of CYP3A56, CYP3A13, and CYP3A31 in BDEG group were 10.27, 2.38, and 2.01 times higher than those in the CG group, respectively. Similarly, a previous study reported that the CYP3A gene expression level was significantly increased in HepG2 cells exposed to BDE-47 (Hu et al., 2014). In mammals, expression of the CYP3A gene can be regulated by a wide range of chemicals, such as PBDEs, pregnanes, and rifampicin (Hu et al., 2014). The regulation mechanism may allow these chemicals to activate the pregnane X receptor (PXR), which is a nuclear transcription factor characterized by a ligand-binding domain and a DNAbinding domain. PXR can regulate CYP3A gene expression by binding to its ligand, translocating to the nucleus, and interacting with aresponse element in the promoter (Watkins et al., 2001). In invertebrates, the CYP3A gene is upregulated in scallops Chlamys farreri when exposed to benzo[a]pyrene, a type of polycyclic aromatic hydrocarbon (Cai et al., 2014). It is not known how these organic pollutants regulate the expression of CYP3A in invertebrates, and whether receptors exist in invertebrates that can regulate CYP3A expression in the same manner as PXR in mammals. These interesting problems will be addressed in our next study.
Glutathione-S-transferase (GST), an important phase Ⅱ enzyme, can catalyze the conjugation of reduced glutathione (L-γ-glutamyl-L-cysteinylglycine) to yield reactive electrophiles of endogenous and exogenous compounds. Thus, the activity level of GST has been used widely as a biomarker of exposure to chemical substances in aquatic organisms. The GST gene transcript level is significantly induced by exposure to PBDEs in razor clams Solen grandis, thereby indicating its involvement in cellular xenobiotic defenses against PBDEs (Yang et al., 2012). PBDEs can also elevate the GST enzyme activity in rock cod Trematomus bernacchii, which suggests that GST is involved with the detoxification of PBDEs (Ghosh et al., 2013). In this study, we detected responses to BDE-47 by five GST isoenzymes, i.e., A-GST, theta 1-GST, pi-GST, sigmaGST, and rho-GST. The mRNA expression levels of these isoenzymes increased after BDE-47 exposure, except for rho-GST, which indicates that several GST isozymes cooperate in BDE-47 detoxification in M.veneriformis.
ABC transporters are phase Ⅲ enzymes, which actively pump the substrates obtained from phases Ⅰ and Ⅱ out of cells. Multi-drugresistance (MDR) proteins are classical ABC transporters, which comprise cell membrane proteins known as P-glycoproteins. In mammals, numerous studies have shown that the failure of cancer chemotherapy is due mainly to drug efflux via pumping by MDR proteins, which is regulated by the activation of AhR, the constitutive androstane receptor, and PXR (Geick et al., 2001). In addition, MDR genes have been identified in several aquatic organisms and their involvement in general biological defenses against environmental toxicants has been demonstrated (Veldhoen et al., 2009; Milan et al., 2016). In our study, the mRNA expression of MvMDR1 and MvMDR3 in BDE-47 group decreased to 54.9% and 74.1%, respectively, of that in the control group. To the best of our knowledge, the present study is the first to describe the change of MDR gene expression, including MvMDR1 and MvMDR3, in response to BDE-47 in invertebrates, thereby indicating their roles in BDE-47 metabolism and detoxification processes. The regulation of the MDR gene by PXR was shown to be conserved in the zebrafish Danio rerio (Bresolin et al., 2005). The receptor that serves as an activator of MDR gene expression in M.veneriformis will be investigated in our next study.4.2 Genes associated with antioxidant defense
Oxyradical production has been identified as a common pollution-mediated mechanism of toxicity in aquatic organisms. It has been demonstrated that BDE-47 markedly enhances the production of oxygen radicals in aquatic animals which suggests that the cell antioxidant defense pathways are overwhelmed (Ghosh et al., 2013). In the present study, we found that the transcripts of some well-known antioxidant enzyme genes, such as copper/zinc superoxide dismutase, selenium-dependent glutathione peroxidase, and chorion peroxidase, were upregulated in response to BDE-47, which indicates their involvement in the defense against oxidative stress caused by BDE-47. Oxidation is regarded as a PBDE toxicity mechanism and reactive oxygen species (ROS) are often produced through the NF-kappa B signaling pathway (Lv et al., 2015). We found that the NF-kappa B signaling pathway (ko: 04064) related to 95 DEGs was significantly enriched. This may be the signaling pathway that allowed BDE-47 to induce oxidative stress in M. veneriformis.
Heat shock proteins (HSPs) are molecular chaperones, which stabilize protein folding and prevent indiscriminate protein interactions by sequestering unfolded proteins. HSPs are stressrelated proteins, which are induced by various chemical, physical, and pathogenic stresses, and they are used as general indicators of non-specific stress. HSP90 (83–90 kDa), HSP70 (66–78 kDa), HSP60, HSP47, and small heat shock protein (sHSP; 15– 43 kDa) are the main HSP families. In particular, the HSP70 family members participate in protein repair, transport, and protection from oxidative stress (Contardo-Jara et al., 2010), and they are used extensively as early biomarkers of environmental stress in diverse organisms. We found that the mRNA expression levels of MvHSP70 12A and MvHSP70 12B were significantly induced by BDE-47. This finding agrees with the induction of HSP 70 gene expression by PBDEs in skin slices of striped dolphins Stenella coeruleoalba (Fossi et al., 2014). Similarly, a field investigation showed that HSP70 transcripts increased in horse mussels Modiolus modiolus (L.) located near a marine municipal waste water outfall with high PBDE concentrations compared with a reference site (Veldhoen et al., 2009). HSP STI1 plays a role in mediating the heat shock response of some HSP70 genes. We detected the upregulation of MVHSP STI1 by BDE-47, which indicates the possible mediatory role of MVHSP STI1 in the HSP70 response to BDE-47, which should be verified in further research. In addition, sHSPs such as HSP40 and HSP22 are considered to belong to a heterogeneous HSPs family, which can counter heat, oxidative, and contamination stress in the clam R. philippinarum, and Corbicula fluminea (Chen et al., 2013; Liu et al., 2015). We detected the downregulation of MVsHSP 40A mRNA expression in response to BDE-47. The different responses of MVHSP family members to BDE-47 suggest their participation in the resistance to BDE-47 stress.4.3 Genes associated with the immune response and apoptosis
Immune system disruption caused by exposure to PBDEs has been demonstrated in several aquatic vertebrates (Ye et al., 2012). However, the immunotoxicity of PBDEs has been reported rarely in aquatic invertebrate. In bivalves, innate immunity is the primary defense against invasion by pathogens because they lack well-developed acquired immunity. We found that the mRNA expression levels of many innate immunity factors, such as lectins, integrins, toll-like receptors, and complement, were regulated by BDE-47, thereby suggesting the participation of these factors in the immune defense against BDE-47 stress in M. veneriformis. Changes in the transcript levels of these immunity factors have also been observed in bivalves in the presence of several environmental contaminants (Cai et al., 2014). In aquatic invertebrates, Ye et al. (2012) detected genderspecific responses of complement genes to BDE-47 in the marine medaka Oryzias melastigma, and they suggested the use of these responses in risk assessments for environmental pollution. Thus, gender must also be considered when determining the immunotoxicity mechanism and in risk assessments for BDE-47.
Apoptosis is a gene-controlled cellular autonomous death mode that maintains homeostasis in organisms. Caspase family members play vital roles in apoptosis. We detected very large increases in caspase transcripts as well as significant changes in the expression levels of other apoptosis-related genes when M. veneriformis was exposed to BDE-47, thereby indicating the response of apoptosis to BDE-47. In murine peritoneal macrophages, BDE-47 was shown to induce concentration-and time-dependent increases in caspase-3 mRNA and protein expression levels by real-time PCR and western blotting, respectively (Lv et al., 2015). In addition, 10 μmol/L BDE-47 caused significant increases in cell apoptosis according to flow cytometry. However, apoptosis was not detected in harbor seal Phoca vitulina cells when exposed to 12 μmol/L BDE-47, BDE-99, and BDE-153 (Frouin et al., 2010). These different results may be attributable to the distinct cell lines employed and the PBDEs concentrations applied. The mechanisms of PBDEinduced apoptosis should clearly be investigated. The involvements of ROS and oxidant stress may provide insights into the mechanism responsible for PBDEinduced apoptosis (Lv et al., 2015). We did not test for ROS, but the mRNA expression levels of antioxidant enzymes exhibited obvious responses to BDE-47 in M. veneriformis. Furthermore, changes in the immune response are possible adverse effects caused by apoptosis. Phagosomes are a key innate immunity factor, which can remodel tissue, eliminate apoptotic cells, and limit the spread of intracellular pathogens. Phagosome pathway (ko: 04145) links immune defense and apoptosis. Based on the pathway analysis in this study, many immune system DEGs such as Toll-like receptors and integrins, and apoptosis DEGs such as cathepsin were detected in the phagosome pathway, which was significantly enriched in M. veneriformis after exposure to BDE-47. Further research would be carried on the apoptosis of immune cells under BDE-47 challenge.5 CONCLUSION
Using Illumina sequencing, we analyzed the transcriptome of M. veneriformis and its response to BDE-47 exposure for the first time in this study. Numerous DEGs related to detoxification, antioxidant defense, immune response, and apoptosis were detected in the clams exposed to BDE-47, thereby indicating the transcriptional complexity in M. veneriformis in response to BDE-47. We provided insights into the toxic effects of BDE-47 as well as its toxicity mechanism. Moreover, many DEGs could be developed as biomarkers to monitor the molecular response to PBDEs by using as M. veneriformis indicator. The regulation and interaction of the DEGs and the toxicity mechanism of PBDEs when combined with multiple chemicals will be investigated in our future research.
|Aarab N, Lemaire-Gony S, Unruh E, Hansen P D, Larsen B K, Andersen O K, Narbonne J F, 2006. Preliminary study of responses in mussel (Mytilus edilus) exposed to bisphenol A, diallyl phthalate and tetrabromodiphenyl ether. Aquatic Toxicology, 78(S): S86–S92.|
|An L H, Hu J Y, Yang M, Zheng B H, Wei A, Shang J J, Zhao X R, 2011. CYP1A mRNA expression in redeye mullets(Liza haematocheila) from Bohai Bay, China. Marine Pollution Bulletin, 62(4): 718–725. Doi: 10.1016/j.marpolbul.2011.01.019|
|Audic S, Claverie J M, 1997. The significance of digital gene expression profiles. Genome Research, 7(10): 986–995. Doi: 10.1101/gr.7.10.986|
|Barja-Fernández S, Míguez J M, Álvarez-Otero R, 2013. Histopathological effects of 2, 2', 4, 4'-tetrabromodiphenyl ether (BDE-47) in the gills, intestine and liver of turbot(Psetta maxima). Ecotoxicology and Environmental Safety, 95: 60–68. Doi: 10.1016/j.ecoenv.2013.05.028|
|Baršienė J, Šyvokienė J, Bjornstad A, 2006. Induction of micronuclei and other nuclear abnormalities in mussels exposed to bisphenol A, diallyl phthalate and tetrabromodiphenyl ether-47. Aquatic Toxicology, 78(S): S105–S108.|
|Brammell B F, McClain J S, Oris J T, Price D J, Birge W J, Elskus A A, 2010. CYP1A expression in caged rainbow trout discriminates among sites with various degrees of polychlorinated biphenyl contamination. Archives of Environmental Contamination and Toxicology, 58(3): 772–782. Doi: 10.1007/s00244-009-9368-x|
|Breitholtz M, Wollenberger L, 2003. Effects of three PBDEs on development, reproduction and population growth rate of the harpacticoid copepod Nitocra spinipes. Aquatic Toxicology, 64(1): 85–96. Doi: 10.1016/S0166-445X(03)00025-0|
|Bresolin T, de Freitas Rebelo M, Bainy A C D, 2005. Expression of PXR, CYP3A and MDR1 genes in liver of zebrafish. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 140(3-4): 403–407.|
|Cai Y F, Pan L Q, Hu F X, Jin Q, Liu T, 2014. Deep sequencingbased transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene. Gene, 551(2): 261–270. Doi: 10.1016/j.gene.2014.09.003|
|Chen H H, Zha J M, Liang X F, Bu J H, Wang M, Wang Z J, 2013. Sequencing and de novo assembly of the Asian Clam (Corbicula fluminea) transcriptome using the illumine GAⅡx method. Plos One, 8(11): e79516. Doi: 10.1371/journal.pone.0079516|
|Chen L G, Zhu B R, Guo Y Y, Xu T, Lee J S, Qian P Y, Zhou B S, 2016. High-throughput transcriptome sequencing reveals the combined effects of key e-waste contaminants, decabromodiphenyl ether (BDE-209) and lead, in zebrafish larvae. Environmental Pollution, 214: 324–333. Doi: 10.1016/j.envpol.2016.04.040|
|Conesa A, Götz S, García-Gómez J M, Terol J, Talón M, Robles M, 2005. Blast2GO:a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics, 21(18): 3 674–3 676. Doi: 10.1093/bioinformatics/bti610|
|Contardo-Jara V, Pflugmacher S, Nützmann G, Kloas W, Wiegand C, 2010. The β-receptor blocker metoprolol alters detoxification processes in the non-target organism Dreissena polymorpha. Environmental Pollution, 158(6): 2 059–2 066. Doi: 10.1016/j.envpol.2010.03.012|
|Fang Y, Yang H S, Liu B Z, Zhang L B, 2013. Transcriptional response of lysozyme, metallothionein, and superoxide dismutase to combined exposure to heavy metals and bacteria in Mactra veneriformis. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 157(1): 54–62.|
|Fang Y, Yang H S, Wang T M, Liu B Z, Zhao H L, Chen M Y, 2010. Metallothionein and superoxide dismutase responses to sublethal cadmium exposure in the clam Mactra veneriformis. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 151(3): 325–333.|
|Fossi M C, Casini S, Bucalossi D, Marsili L, 2008. First detection of CYP1A1 and CYP2B induction in Mediterranean cetacean skin biopsies and cultured fibroblasts by Western blot analysis. MarineEnvironmental Research, 66(1): 3–6.|
|Fossi M C, Casini S, Maltese S, Panti C, Spinsanti G, Marsili L, 2014. An "ex vivo" model to evaluate toxicological responses to mixtures of contaminants in cetaceans:integumentum biopsy slices. Environmental Toxicology, 29(10): 1 107–1 121. Doi: 10.1002/tox.v29.10|
|Frouin H, Lebeuf M, Hammill M, Masson S, Fournier M, 2010. Effects of individual polybrominated diphenyl ether(PBDE) congeners on harbour seal immune cells in vitro. Marine Pollution Bulletin, 60(2): 291–298. Doi: 10.1016/j.marpolbul.2009.09.006|
|Geick A, Eichelbaum M, Burk O, 2001. Nuclear receptor response elements mediate induction of intestinal MDR1 by rifampin. The Journal of Biological Chemistry, 276(18): 14 581–14 587. Doi: 10.1074/jbc.M010173200|
|Ghosh R, Lokman P M, Lamare M D, Metcalf V J, Burritt D J, Davison W, Hageman K J, 2013. Changes in physiological responses of an Antarctic fish, the emerald rock cod(Trematomus bernacchii), following exposure to polybrominated diphenyl ethers (PBDEs). Aquatic Toxicology, 128-129: 91–100. Doi: 10.1016/j.aquatox.2012.11.019|
|Giannetti M, Casini S, Marsili L, Maltese S, Campani T, Carletti L, Fossi M C, 2012. First evidence of induction of CYP1A in Caretta caretta skin biopsy slices treated with polycyclic aromatic hydrocarbons, organochlorines and polybrominated diphenyl ethers. Comparative Biochemistry and Physiology Part A:Molecular & Integrative Physiology, 163(S): S18–S19.|
|Grabherr M G, Haas B J, Yassour M, Levin J Z, Thompson D A, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q D, Chen Z H, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren B W, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A, 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology, 29(7): 644–652. Doi: 10.1038/nbt.1883|
|Hu X N, Zhang J Q, Jiang Y S, Lei Y X, Lu L G, Zhou J, Huang H Y, Fang D K, Tao G H, 2014. Effect on metabolic enzymes and thyroid receptors induced by BDE-47 by activation the pregnane X receptor in HepG2, a human hepatoma cell line. Toxicology in Vitro, 28(8): 1 377–1 385. Doi: 10.1016/j.tiv.2014.07.004|
|Irato P, Piccinni E, Cassini A, Santovito G, 2007. Antioxidant responses to variations in dissolved oxygen of Scapharca inaequivalvis and Tapes philippinarum, two bivalve species from the lagoon of Venice. Marine Pollution Bulletin, 54(7): 1 020–1 030. Doi: 10.1016/j.marpolbul.2007.01.020|
|Iseli C, Jongeneel C V, Bucher P. 1999. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. In: Park M ed. Proceedings of the Seventh International Conference on Intelligent Systems for Molecular Biology. AAAI Press, Menlo Park. p. 138-148. http://ci.nii.ac.jp/ncid/BA45988031|
|Ji C L, Wu H F, Wei L, Zhao J M, Yu J B, 2013. Proteomic and metabolomic analysis reveal gender-specific responses of mussel Mytilus galloprovincialis to 2, 2', 4, 4'-tetrabromodiphenyl ether (BDE 47). Aquatic Toxicology, 140-141: 449–457. Doi: 10.1016/j.aquatox.2013.07.009|
|Jin J, Liu W Z, Wang Y, Tang X Y, 2008. Levels and distribution of polybrominated diphenyl ethers in plant, shellfish and sediment samples from Laizhou Bay in China. Chemosphere, 71(6): 1 043–1 050. Doi: 10.1016/j.chemosphere.2007.11.041|
|Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y, 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Research, 36(D1): D480–D484.|
|Law R J, Allchin C R, de Boer J, Covaci A, Herzke D, Lepom P, Morris S, Tronczynski J, de Wit C A, 2006. Levels and trends of brominated flame retardants in the European environment. Chemosphere, 64(2): 187–208. Doi: 10.1016/j.chemosphere.2005.12.007|
|Liu B, Jiang G F, Zhang Y F, Li J L, Li X J, Yue J S, Chen F, Liu H Q, Li H J, Zhu S P, 2011. Analysis of transcriptome differences between resistant and susceptible strains of the citrus red mite Panonychus citri (Acari:Tetranychidae). PLoS One, 6(12): e28516. Doi: 10.1371/journal.pone.0028516|
|Liu T, Pan L Q, Jin Q, Cai Y F, 2015. Differential gene expression analysis of benzo(a)pyrene toxicity in the clam, Ruditapes philippinarum. Ecotoxicology and Environmental Safety, 115: 126–136. Doi: 10.1016/j.ecoenv.2015.02.007|
|Lv Q Y, Wan B, Guo L H, Zhao L X, Yang Y, 2015. In vitro immune toxicity of polybrominated diphenyl ethers on murine peritoneal macrophages:apoptosis and immune cell dysfunction. Chemosphere, 120: 621–630. Doi: 10.1016/j.chemosphere.2014.08.029|
|Miao J J, Pan L Q, Zhang W H, Liu D, Cai Y F, Li Z, 2014. Identification of differentially expressed genes in the digestive gland of manila clam Ruditapes philippinarum exposed to BDE-47. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 161: 15–20.|
|Milan M, Matozzo V, Pauletto M, Di Camillo B, Giacomazzo M, Boffo L, Binato G, Marin M G, Patarnello T, Bargelloni L, 2016. Can ecological history influence response to pollutants? Transcriptomic analysis of Manila clam collected in different Venice lagoon areas and exposed to heavy metal. Aquatic Toxicology, 174: 123–133. Doi: 10.1016/j.aquatox.2016.02.024|
|Moore M N, Allen J I, 2002. A computational model of the digestive gland epithelial cell of marine mussels and its simulated responses to oil-derived aromatic hydrocarbons. Marine Environmental Research, 54(3-5): 579–584. Doi: 10.1016/S0141-1136(02)00166-6|
|Mortazavi A, Williams B A, McCue K, Schaeffer L, Wold B, 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods, 5(7): 621–628. Doi: 10.1038/nmeth.1226|
|Peters A K, Sanderson J T, Bergman Å, van den Berg M, 2006. Antagonism of TCDD-induced ethoxyresorufin-Odeethylation activity by polybrominated diphenyl ethers(PBDEs) in primary cynomolgus monkey (Macaca fascicularis) hepatocytes. Toxicology Letters, 164(2): 123–132. Doi: 10.1016/j.toxlet.2005.12.002|
|Santovito G, Boldrin F, Irato P, 2015. Metal and metallothionein distribution in different tissues of the Mediterranean clam Venerupis philippinarum during copper treatment and detoxification. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 174-175: 46–53.|
|Santovito G, Piccinni E, Cassini A, Irato P, Albergoni V, 2005. Antioxidant responses of the Mediterranean mussel, Mytilus galloprovincialis, to environmental variability of dissolved oxygen. Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology, 140(3-4): 321–329.|
|Schecter A, Johnson-Welch S, Tung K C, Harris T R, Päpke O, Rosen R, 2007. Polybrominated diphenyl ether (PBDE) levels in livers of U. S. human fetuses and newborns. Journal of Toxicology and Environmental Health, Part A, 70(1): 1–6.|
|Stockholm Convention. 2015. http://chm.pops.int.com.|
|Veldhoen N, Lowe C J, Davis C, Mazumder A, Helbing C C, 2009. Gene expression profiling in the deep water horse mussel Modiolus modiolus (L.) located near a marine municipal wastewater outfall. Aquatic Toxicology, 93(2-3): 116–124.|
|Watkins R E, Wisely G B, Moore L B, Collins J L, Lambert M H, Williams S P, Willson T M, Kliewer S A, Redinbo M R, 2001. The human nuclear xenobiotic receptor PXR:structural determinants of directed promiscuity. Science, 292(5525): 2 329–2 333. Doi: 10.1126/science.1060762|
|Wollenberger L, Dinan L, Breitholtz M, 2005. Brominated flame retardants:activities in a crustacean development test and in an ecdysteroid screening assay. Environmental Toxicology and Chemistry, 24(2): 400–407. Doi: 10.1897/03-629.1|
|Yang J L, Wei X M, Xu J, Yang D L, Liu X Q, Yang J M, Fang J H, Hu X K, 2012. A sigma-class glutathione S-transferase from Solen grandis that responded to microorganism glycan and organic contaminants. Fish & Shellfish Immunology, 32(6): 1 198–1 204.|
|Ye J, Fang L, Zheng H K, Zhang Y, Chen J, Zhang Z J, Wang J, Li S T, Li R Q, Bolund L, Wang J, 2006. WEGO:a web tool for plotting GO annotations. Nucleic Acids Research, 34(S2): W293–W297.|
|Ye R R, Lei E N Y, Lam M H W, Chan A K Y, Bo J, van de Merwe J P, Fong A C C, Yang M M S, Lee J S, Segner H E, Wong C K C, Wu R S S, Au D W T, 2012. Genderspecific modulation of immune system complement gene expression in marine medaka Oryzias melastigma following dietary exposure of BDE-47. Environmental Science and Pollution Research, 19(7): 2 477–2 487. Doi: 10.1007/s11356-012-0887-z|
|Yu L Q, Han Z H, Liu C S, 2015. A review on the effects of PBDEs on thyroid and reproduction systems in fish. General and Comparative Endocrinology, 219: 64–73. Doi: 10.1016/j.ygcen.2014.12.010|