Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene

GENE-39929; No. of pages: 10; 4C: Gene xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

3Q4

Yuefeng Cai, Luqing Pan ⁎, Fengxiao Hu, Qian Jin, Tong Liu

4

Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao 266003, PR China

5

a r t i c l e

6 7 8 9 10

Article history: Received 5 June 2014 Received in revised form 11 August 2014 Accepted 3 September 2014 Available online xxxx

11 12 13 14 15

Keywords: Chlamys farreri Digital gene expression Benzo[a]pyrene Detoxification mechanism

O R O

i n f o

F

2

Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene

a b s t r a c t

E

D

P

Whole-genome transcriptome measurements are pivotal for characterizing molecular mechanisms of chemicals and predicting toxic classes, such as genotoxicity and carcinogenicity, from in vitro and in vivo assays. We analyzed the dynamic defense transcriptome responsive to Chlamys farreri upon exposure to benzo[a]pyrene (BaP) using a digital gene expression (DGE) approach. Following exposure, 251 and 177 genes were upregulated, and 142 and 300 genes were down-regulated at 3 days post-exposure and 10 days post-exposure, respectively. The differentially expressed genes were related to toxicological response, oxidative stress and the metabolism of proteins and fats. Of these genes, most genes up-regulated at the early stage of exposure tended to be constantly down-regulated at the later stage whereas the landscape of the up- or down-regulated genes differed significantly at the two time points investigated. Functional enrichment analyses show that RNA-seq yields more insight into the biological mechanisms related to the toxic effects caused by BaP, i.e., two to fivefold more affected pathways and biological processes. Besides, we observed a change in the expression of ten genes which are important and differentially-expressed detoxification-related genes, and this was subsequently confirmed via quantitative real-time PCR. Our results provide evidence that RNA-seq is a powerful tool for toxicology and is capable of generating novel and valuable information at the transcriptome level for characterizing deleterious effects caused by BaP. © 2014 Published by Elsevier B.V.

C

T

1

31

E

35 33 32 34

1. Introduction

37 38

Polycyclic aromatic hydrocarbons (PAHs) are aromatic hydrocarbons with two or more fused benzene rings with natural as well as anthropogenic sources. They are common coastal environmental organic pollutants derived from petroleum, ocean shipping and emissions from fossil fuel utilization and conversion processes (Menzie et al., 1992). Due to their ubiquitous occurrence, recalcitrance, bioaccumulation potential and carcinogenic activity, the PAHs have gathered more attention in the field of aquatic environment. For instance, the toxicity and carcinogenicity of PAHs can be a direct consequence of metabolic activation by cytochrome P450, or indirectly modulated by the increased production of reactive oxygen species (ROS) and transcriptional regulation of several gene systems (Meyer et al., 2002). Bivalve is one of the most important marine commercial species in China, contributing 75–80% of the total aquaculture production in the last three years (Yue and Wang, 2012). Due to their sessile nature,

43 44 45 46 47 48 49 50 51 Q7

R

N C O

41 42

U

39 40

R

36

Abbreviations: BaP, benzo[a]pyrene; DGE, digital gene expression; PAHs, polycyclic aromatic hydrocarbons; ROS, reactive oxygen species; DMSO, dimethyl sulfoxide; Cq, cycle threshold; RPKM, reads per kilobase of per million mapped reads; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. ⁎ Corresponding author at: Fisheries College, Ocean University of China, Yushan Road 5, Qingdao 266003, China. E-mail address: [email protected] (L. Pan).

filter-feeding habits, and pollutant bioconcentration, bivalves are widely used as marine pollution sentinels (Goldberg et al., 1978). There are persistent concerns regarding the release of PAHs into aquatic environments and the potential harmful effects on bivalves, even when released at low levels (Geffard et al., 2003; Frouin et al., 2007). Nevertheless, many of these toxic responses (e.g. oxidative stress, enzymatic activation/inhibition, genotoxicity) are common to several contaminants, including diol epoxides or ROS (Baussant et al., 2009). However, the detailed detoxification processes of these pollutants on bivalve is unclear, especially overall response at the level of transcription. Understanding the effects of PAHs on bivalves is required to effectively establish methods of assessing and predicting the toxicity of PAHs. Nevertheless, due to the lack of available genomic resources such as genome and transcriptome sequences, an overall understanding of the mechanism of action of PAHs on bivalves is a high priority. Benzo[a]pyrene (BaP), one of carcinogenic PAHs, is especially interesting for transcriptomic studies, as its genotoxic and nongenotoxic properties. The analysis of changes in gene expression represents a powerful tool to characterize immediate cell responses to stressors, constituting an early warning of the effect of contaminants. These changes also represent a useful complement to existing monitoring methods to study the effects of toxicants at the biochemical level (Faria et al., 2009, 2010). Gene expression profiling is currently widely applied for unraveling mechanisms of toxic properties of chemicals as well as for

http://dx.doi.org/10.1016/j.gene.2014.09.003 0378-1119/© 2014 Published by Elsevier B.V.

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

16 17 18 19 20 21 22 23 24 25 26 27 28 Q5 29 30

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 Q8 75

112 113

Adult C. farreri with the shell length of 7.25 ± 0.5 cm, were obtained from Pacific Corner (Yellow Sea, Qingdao, China). All the scallops were acclimated to laboratory condition in aquarium (1 L water per scallop) for one week before formal experiments. They were randomly distributed and cultured under filtered seawater at a salinity of 30‰ with a temperature of 17 ± 1 °C and a pH of 8.1 in 30 × 40 × 50 cm3 tanks. Each replicate contained 30–40 scallops (three replicates were performed). Water was renewed completely and the scallops were fed with dried powder of Spirulina platensis (30 mg for each individually) daily. BaP (CAS No. 50-32-8; 97.5% purity, SIGMA, USA) was first dissolved in dimethyl sulfoxide (DMSO), and then added to seawater to achieve a final DMSO concentration of 0.001% (v/v). The water solubility of BaP is 16 nM (Schirmer et al., 1998) and the concentration of BaP in natural seawater from Qingdao was 0.155 ng/L (Pan et al., 2008). But BaP in Daya Bay (Shenzhen, China) reached 655 ng/L (Qiu et al., 2004). The pollutant nominal exposure concentrations (control, 1 μg/L) were used for exposure referring to their water solubility and the concentration in the polluted offshore. Experimental conditions (salinity, pH, temperature, scallop density, feeding) were the same as those used for acclimation, and exposure media were renewed completely daily. The exposure experiment lasted for ten days. Six scallops for each replicates were sampled at days three and ten. Digestive gland tissue was collected from each group and immediately frozen in liquid nitrogen and then stored at −80 °C until use. The samples were named as CfX1_1\CfX1_1_1(control and treated after 3 days) and CfX1_2\CfX1_2_1(control and treated after 10 days).

101 102 103 104 105 106 107 108

114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137

C

99 Q9 100

E

97 98

R

95 96

R

93 94

O

91 92

C

89 90

N

87 88

U

85 86

2.3. Sequence Reads, Mapping and Assembly

175

Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, the clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N and low quality reads from raw data. At the same time, quality parameters of clean data including Q30, GC-content and sequence duplication level were used for data filtering. All the succeeding analyses were carried out using high quality clean data. Reference transcriptome and unigene model annotation files were built and provided by our laboratory directly. Index of the reference unigene was built using Bowtie v0.12.9 (Langmead et al., 2009) and single-end clean reads were aligned to the reference transcriptome using RSEM (Li and Dewey, 2011). Clean reads were aligned to the reference transcriptome through RSEM, then duplicated reads and multimapped reads were filtered from the alignment results in order to eliminate the PCR interference and ambiguous mapping. The clean reads have been submitted to NCBI Short Read Archive under the accession number of SRP034886.

176 177

2.4. Quantification and Differential Expression Analysis of Transcripts

193

RSEM was used to count the read numbers mapped to each unigene. And then the RPKM of each transcript was calculated based on the length of the transcript and read count mapped to this transcript. RPKM, reads per kilobase of per million mapped reads, considers the

194

F

2.1. Scallops and Exposure

83 84

139 140

O

111

82

Trizol Reagent was used to isolate total RNA from tissues according to the manufacturer's instructions (Invitrogen, USA). RNA degradation and contamination were assessed on 1% agarose gels. RNA concentration was measured using Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, CA, USA). RNA purity and integrity were checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively. A total amount of 3 μg RNA was used as input material for the RNA sample preparations. Finally, four samples with RNA integrity number (RIN) values above 8 were used for library construction. Sequencing libraries were generated using the IlluminaTruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, USA) following the manufacturer's recommendations and four index codes were added to attribute sequence to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in Illumina proprietary fragmentation buffer. First-strand cDNA was synthesized using random oligonucleotides and SuperScript II. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities and enzymes were removed. After adenylation of 3′ ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. In order to select cDNA fragments of 200 bp in length the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 10 cycle PCR reaction. Products were purified (AMPure XP system) and quantified using the Agilent high sensitivity DNA assay on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq2000 platform and 100 bp single-end reads were generated.

R O

2. Materials and Methods

80 81

138

P

110

78 79

2.2. RNA Extraction, Library Preparation and Sequencing

T

109

predicting toxicity of chemical compounds (Kim et al., 2005; Mathijs et al., 2009; McHale et al., 2010; Paules et al., 2011; Tian et al., 2013). Brown et al. firstly used SSH and a nylon microarray to study Mytilus edulis gene expression profiling in response to BaP exposure, results showed that 112/250 up-regulated and 25/55 down-regulated clones were positive for differential expression and the transcripts isolated were from a diverse range of genes involved in general stress, oxidative stress, cell adhesion, transcriptional and translational regulation etc. (Brown et al., 2006). However, the gene chip and SSH techniques have become challenging when dealing with more inventories of gene species, e.g. the genes represented by unspecific probe sets and with low expression levels cannot be readily detected. Next-generation sequencing (NGS) technologies provide powerful alternative strategies for transcriptome analysis. RNA sequencing is increasingly being used for global gene expression profiling because it allows unbiased quantification of expression levels of transcripts with a higher sensitivity and broader genome coverage than microarrays (Mortazavi et al., 2008). RNA-seq and DGE technology have been applied to stress-related gene and signaling pathway analysis of several bivalves such as yesso scallop (Patinopecten yessoensis), pacific oyster (Crassostrea gigas) and mussels (Mytilus galloprovincialis) (Figueras et al., 2012; Gomes et al., 2013; Wang et al., 2013; Zhao et al., 2012). In the present study, we examined the dynamic transcriptome responses to BaP exposure of Chlamys farreri for the first time using the Illumina's sequencing technology. Considering individual monitoring of the scallop responses to BaP exposure, four libraries, 3 day samples (control and BaP exposure) and 10 day samples (control and BaP exposure) were established from the digestive gland of scallops. The goal of this study is to: (1) establish a platform to characterize gene expression on a larger scale that can facilitate detailed characterization on genes regulating crop toxicological responses to BaP; (2) quantify the transcript abundance in the species under BaP exposure; and (3) uncover the networks of genes enriched for regulating C. farreri resistance to the PAHs.

D

76 77

Y. Cai et al. / Gene xxx (2014) xxx–xxx

E

2

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174

178 179 180 181 182 183 184 185 186 187 188 189 190 191 192

195 196 197

Y. Cai et al. / Gene xxx (2014) xxx–xxx

Name of gene

Primer sequences used for qRT-PCR (5′–3′)

t1:4

Heat shock protein 70 (hsp70)

F: GCCGACAGCCCAGTAAAT R: AAACCTGTGATTTCTTGGTAGG F: GTATTTGCCCATCGGAGAG R: CCCGTTCGTCCAAGAATACA F: ATACCTTCCTGGTGACATGTTC R: TGTGTCCCCCTCCTGTGA F: CGGGATGGTGGTGACTGT R: TAAAGCAAGCCGCATAGC F: GAACCCAGTCAGACCCATTAC R: ACTAGATAGTCATACCCAACCTCTT F: AGTTTGGTTTGGCGGGAG R: AGCAAATGTTGTTAGAACTGAATC F: GCTACTGGACCTCCGTGTCT R: GCAACAGTCAAAGCCACGT F: AGTTCGGACATCAGGAGAACA R: CATGAAGGAGGTTGGATCGT F: GACTTCAGCGTCGTTATCTTG R: GATCCCATCAACCTGAAACAG F: ACGGGCAGGCAGAACTAC R: CTTGGTCTGTAGAGTAACCAGTGT F: ATTACGTCCCTGCCCTTTGTAC R: ATGATCCTTCCGCAGGTTCA

2.5106 2.2638 2.0111 2.1286 1.9934

2.0695 2.0486

3. Results 2.0647

2.5. GO and KEGG Enrichment Analysis of Differentially Expressed Genes

C

E

205 206

T

208

203 204

209 210

222

U

N C O

R

R

Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the GOseq R package (Young et al., 2010), in 211 which gene length bias was corrected. GO terms with corrected p value 212 Q10 less than 0.05 were considered significantly enriched by differentially 213 expressed genes. 214 Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database re215 source for understanding high-level functions and utilities of the biolog216 ical system, such as the cell, the organism and the ecosystem, from 217 molecular-level information, especially large-scale molecular datasets 218 Q11 generated by genome sequencing and other high-throughput experi219 mental technologies (http://www.genome.jp/kegg/). We used KOBAS 220 2.0 (Mao et al., 2005) software to test the statistical enrichment of dif221 ferential expression genes in KEGG pathways.

223

2.6. Validation of Differentially Expressed Genes by Quantitative Real-Time PCR

224 225

Differentially expressed genes identified by the above described method were validated using quantitative real-time PCR (qPCR). The

t2:1 t2:2

Table 2 Summary of RNA-seq and mapping results.

226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244

3.1. Digital Gene Expression (DGE) Library Sequencing and Mapping

245

Solexa/Illumina DGE analysis was performed to obtain a global view of C. farreri transcriptome after being exposed to BaP. Four C. farreri DGE libraries were sequenced corresponding to 3 days and 10 days after BaP treatment and control groups. The total number of reads per library ranged from 12.6 to 17.3 million (Table 2). After removing the low quality reads, the number of clean reads ranged from 12.49 to 17.15 million (98.90–99.07%). The clean data sets are available at the NCBI SRA with the accession number: SRP034886. To elucidate the molecular events underlying DGE profiles, all clean reads were aligned to the reference transcriptome database of C. farreri established by our laboratory. Gene annotation was performed by tag mapping analysis using the 92,999 unigenes from RNA-seq based on transcriptome analysis as the reference transcript database (raw data SRA accession number: SRP018044). Approximately 74.53% to 77.89% of the distinct reads were mapped uniquely to reference sequences (Table 2). Finally, read mapping to the database generated 12.8, 9.7, 9.7, and 9.8 million mapped clean reads for the four libraries, respectively (Table 2). Saturation of gene detection occurred when the sequencing counts reached 8 M or higher. In this case, the gene expression level was determined by calculating the number of mapped reads for each gene locus and then normalizing it to the number of RPKM (Fig. 1).

246

3.2. Effects of Two Points of BaP Exposure

267

Neither three days of exposure nor ten days of exposure caused lethality to C. farreri, and there were no obvious signs of distress such as an inability to attach its position in the water column or changes in feeding behavior. Therefore, changes in abundances of transcripts in digestive glands of scallops exposed to BaP might not be indicative of acute toxicity but might be indicators of the effects of BaP or PAHs in coastal environment. To calculate the expression abundance of each gene between CfX1_1 and CfX1_1_1 or CfX1_2 and CfX1_2_1, the DEGseq package was used (Anders and Huber, 2010). After the data were filtered using strict

268 269

P

2.0314

207

201 202

O

1.9853

effect of sequencing depth and transcript length for the read count at the same time (Mortazavi et al., 2008). Prior to differential gene expression analysis, for each sequenced library, the read counts were standardized by TMM, which removes the potential of differentially expressed genes. Differential expression analysis of two conditions was performed using the DESeq R package (1.12.0) (Sun et al., 2013). The p values were adjusted using the Benjamini Hochberg method. |log2(FoldChange)| N1 and qvalue b 0.005 were set as the threshold for significantly differential expression.

199 200

F

2.0025

D

198

Monoamine oxidase (moA) Cytochrome P450 2D10 (CYP 2D10) Cytochrome P450 3A2 (CYP 3A2) Sulfide:quinone oxidoreductase (NQO1) Glutathione S-transferase (GST) Superoxide dismutase (SOD) Glutathione peroxidase 2 (GPx 2) Glutathione peroxidase 3 (GPx 3) Multidrug resistance protein 1-like (MRP1) 18 s rRNA

Efficiency

R O

t1:3

t1:5 t1:6 t1:7 t1:8 t1:9 t1:10 t1:11 t1:12 t1:13 t1:14 t1:15 t1:16 t1:17 t1:18 t1:19 t1:20 t1:21 t1:22 t1:23

RNA sample was reverse-transcribed using PrimeScript RT-PCR Kit (TaKaRa, Dalian, China). 18s rRNA was used as a reference control. The qRT-PCR was performed using 2× SYBR Green master mix (Takara, Dalian, China) on the Real-Time Thermal Cycler 5100 (Thermo Scientific, Finland). The reaction was performed using the following conditions: denaturation at 95 °C for 3 min, followed by 40 cycles of amplification (95 °C for 10 s, 57 °C for 20 s, and 72 °C for 30 s). The primer sequences of the tested genes used in the qRT-PCR study are presented in Table 1. Cq values were determined based on two biological replicates each with two technical replicates. The relative quantification of a target molecule relative to the reference gene was quantified according to the Mean Normalized Expression (MNE) method of Simon (2003). The values of expression were imported to Microsoft Excel and data analyses were performed using SigmaPlot 12.5. All data were presented as means ± S.D. (n = 3) of the relative mRNA expression. Differences were considered significant at p b 0.05. Significant differences were tested by a one-way ANOVA with a Tukey test (a post hoc test). The level of p b 0.05 was considered to be significant.

Table 1 Primers used in SYBR quantitative RT-PCR analysis.

E

t1:1 t1:2

3

t2:3

Sample name

Raw reads

Clean reads

Error rate (%)

Q30 (%)

GC content (%)

Total mapped

Mapped ratio (%)

t2:4 t2:5 t2:6 t2:7

CfX1_1 CfX1_1_1 CfX1_2 CfX1_2_1

17,308,487 12,611,437 12,711,432 12,685,692

17,147,258 12,489,691 12,571,017 12,558,021

0.03 0.03 0.03 0.03

92.13 92.07 91.94 92.04

42.39 43.11 42.92 43.07

12,779,920 9,718,659 9,688,809 9,781,005

74.53 77.81 77.07 77.89

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266

270 271 Q12 272 273 274 Q13 275 276 277 Q14

Y. Cai et al. / Gene xxx (2014) xxx–xxx

283 284 285 286 287 288 289 290

3.3. Functional Analysis of Differentially Expressed Genes

292 293

To investigate changes in the patterns of gene expression after BaP exposure, the percentage of genes in GO categories was analyzed for CfX1_1 vs. CfX1_1_1 and CfX1_2 vs. CfX1_2_1. The significantly enriched GO terms are shown in Fig. 2. The oxidoreductase activity is the significantly enriched GO term of molecular function ontology after exposure to BaP. For the ontology of the cellular component, 13 GO terms are significantly enriched after three days of exposure, such as the mitochondrial matrix, succinate–CoA ligase complex, caspase complex, endoplasmic reticulum chaperone complex, gap junction and the focal adhesion. But there are only 6 GO terms significantly enriched after ten days of exposure, as heterogeneous nuclear ribonucleoprotein complex, Arp2/3 protein complex and protein–lipid complex. To better understand the functions of differentially expressed genes, GO category enrichment analysis was performed using a p-value of ≤0.05 as the cutoff. Table 4 lists the GO enrichment analysis results. According to biological function, some stress-related progresses were enriched during BaP exposure. The most representative GO terms included “xenobiotic metabolic process”, “response to xenobiotic stimulus”, “cellular response to xenobiotic stimulus”, “xenobiotic catabolic

303 304 305 306 307 308 309 310 311

t3:1 t3:2 t3:3

R

R

O

301 302

C

299 300

N

297 298

U

295 296

E

291

294

Table 3 Significantly differentially expressed unigenes with BLAST identities based on hits against the reference (|log2(FoldChange)| N 1 and q-value b 0.005).

t3:4 Q1 t3:5 t3:6 t3:7

F

O

R O

P

281 282

T

279 280

statistical criteria (|log2(FoldChange)| N 1 and q-value b 0.005), a total of 870 significantly changed gene entities were observed following exposure to BaP for 3 days and 10 days (Table 3). After being exposed to BaP for three days, there were 251 up-regulated genes and 142 downregulated genes (Table 3). As Table 3 shows, ten days of exposure to BaP resulted in greater abundances of 177 transcripts and lesser abundances of 300 transcripts. A down-regulated trend in the number of differentially expressed genes was observed from 3 days to 10 days. Most of the genes were up-regulated in 3 days libraries and up-regulated genes showed decreasing expression patterns compared with the number of 10 day libraries (Table 3). In general, there was very little overlap between the two points in transcripts that were of greater or lesser abundances.

C

278

D

Fig. 1. Abundant distribution of unigenes. The Y-axis represents the gene number; the X-axis represents RPKM (reads per kilobase of exon per million mapped sequence reads) range of genes. 0–10: Low expression genes; 10–500: moderate expression genes; N500: high-expression genes.

process”, “glutathione metabolic process”, “immune response” and “inflammatory response”. Furthermore, some GO terms of biological functions changed gradually during BaP exposure. “Response to hormone stimulus”, “G1 DNA damage checkpoint”, and “regulation of lymphocyte proliferation” were the predominant enrichment GO terms in the 3 day library, which were absent in the 10 day library. The GO terms “protein acetylation”, and “protein heterooligomerization” were enriched only in the 10 day library. Some other biological processes, such as “juvenile hormone mediated signaling pathway”, “vitamin metabolic process”, “deoxyribonucleoside diphosphate metabolic process”, “methylglyoxal metabolic process”, and “oxidation–reduction process” were enriched during BaP exposure. According to molecular function, the most represented GO terms were “binding”, “catalytic activity”, and “transporter activity”, including bHLH transcription factor binding, activating transcription factor binding, p53 binding, SMAD binding, dioxygenase activity, transaminase activity, alkaline phosphatase activity, sterol esterase activity, NADPH: quinone reductase activity, oxidoreductase activity, acting on NAD(P) H, quinone or similar compound as acceptor, thioredoxin-disulfide reductase activity, glutamate dehydrogenase [NAD(P)+] activity, lactoylglutathione lyase activity, O-methyltransferase activity, acetylCoA C-acetyltransferase activity, norepinephrine transmembrane transporter activity, norepinephrine:sodium symporter activity, neurotransmitter transporter activity, organic hydroxy compound transmembrane transporter activity, and glycine hydroxymethyltransferase activity. Some GO terms of molecular function also changed during BaP exposure. Further, KEGG enrichment was performed according to C. farreri transcriptome annotated information. Among all the genes with KEGG annotation, 350 differentially expressed genes were observed between CfX1_1 and CfX1_1_1. Most of them were involved in PI3K-Akt signaling pathway, focal adhesion, ECM-receptor interaction, p53 signaling pathway. Between CfX1_2 and CfX1_2_1, 480 differentially expressed genes were annotated with KEGG pathway annotation. The main enriched metabolism pathways concentrated on focal adhesion, PI3K-Akt signaling pathway, ECM-receptor interaction, biosynthesis of secondary metabolites, protein digestion and absorption. Cellular process including lysosome and peroxisome were also enriched (Table 4). Between CfX1_1_1 and CfX1_2_1, 480 differentially expressed genes were annotated with KEGG pathway annotation. The main enriched metabolism pathways concentrated on focal adhesion, PI3K-Akt signaling pathway, ECM-receptor interaction, biosynthesis of secondary metabolites, protein digestion and absorption. Cellular processes including lysosome and peroxisome were also enriched (Table 5).

E

4

Up-regulated Down-regulated Total

3 days

10 days

3 days vs. 10 days

251 142 393

177 300 477

75 136 211

312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356

3.4. Expression Patterns of Genes Related Organic Compounds 357 Detoxification 358 Based on GO and KEGG analysis, the expression patterns of genes differentially expressed throughout BaP stress which involved defense in mechanisms against xenobiotic contamination were further analyzed. Table 3 showed that most of the expression of genes increased, and then decreased under more serious oxidative stress induced by BaP exposure, such as quinone oxidoreductase, dual oxidase 1, O-methyltransferase, Thioredoxin reductase 3, glutathione peroxidase, glutamate dehydrogenase 1, major vault protein, and amino acid transporter. In contrast, some metallothionein, UDP-glucuronic acid decarboxylase 1 and thyroglobulin genes were up-regulated all the time (Table 3).

359

3.5. Confirmation of Differentially Expressed Genes by qRT-PCR

370

360 361 362 363 364 365 366 367 368 369

Ten genes were selected for qRT-PCR analysis to validate the DGE 371 data. Expressions of all ten genes detected by qRT-PCR were the same 372 with the DGE data. The results showed that the mRNA level of 373

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

5

N C O

R

R

E

C

T

E

D

P

R O

O

F

Y. Cai et al. / Gene xxx (2014) xxx–xxx

Fig. 2. Functional categorization of differentially expressed genes based on known genes in the Uniprot database. The Y-axis shows the 2nd level term of Gene Ontology; X-axis shows number of genes in differentially expressed genes. A: CfX1_1\CfX1_1_1 (control and treated after 3 days), B: CfX1_2\CfX1_2_1 (control and treated after 10 days).

374

380

4. Discussion

381

Marine bivalve, such as C. farreri, is an important marine commercial species that accounts for 15–20% of the total commercial production in China. They have a set of unique capabilities to adapt to the complicated conditions owing to their habitats, living habits and feeding ways. Simultaneously, bivalves are widely used as model organisms in disease control as well as in monitoring programs because of their ecological habits. The transcriptome information can provide the molecular basis for this scallop model.

376 377

382 383 384 385 386 387 388

U

378 379

CYP2D10 and GPx2 was significantly up-regulated after three days of exposure (Fig. 3A), but NQO1 and CYP3A2 were increased after ten days of exposure (Fig. 3B). In contrast, the expression level of GPx3 and MRP1 decreased after three days of exposure, and CYP2D10, GST3, SOD, moA and hsp70 were significantly down-regulated after ten days of exposure (Fig. 3C).

375

4.1. Biotransformation and Detoxification of BaP

389

Benzo[a]pyrene (BaP), a member of the polycyclic aromatic hydrocarbon (PAH) family, is a powerful mutagen and carcinogen. PAHs pollution becomes more serious in the marine environment due to the development of offshore petroleum industry and steamboat transportation. Pan et al. measured the average content of BaP in bottom water of the Jiaozhou Bay which reached 0.155 ng/L (Pan et al., 2008) and in Daya Bay (Shenzhen, China) it reached 655 ng/L (Qiu et al., 2004). In this study, the effects of BaP on C. farreri were revealed at the transcriptome levels using DGE analysis. The different expression genes were involved in the various biological functions. The genes involved in the amino acids, glycerolipid, glycerophospholipid metabolism and glycolysis/gluconeogenesis were significantly enriched. This indicated that BaP can affect the basic metabolism in the scallop. Changes also appeared in the metabolite levels in the scallop and clam after exposure to BaP (Liu et al., 2012; Wang et al., 2011).

390 391

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

392 Q15 393 394 395 396 Q16 397 398 399 400 401 402 403 404

6

Table 4 Fold-changes in abundances of transcripts of genes after exposure to BaP. Gene names and log2 (fold change) in abundance of transcripts determined by RNAseq are given. Transcript

3 days-log2 (fold change) (RNAseq)

10 days-log2 (fold change) (RNAseq)

Cytochrome P450 2B5 Cytochrome P450 2U1 Cytochrome P450 2D10 Cytochrome P450 3A2 Catechol-O-methyltransferase family 3 Glutathione S-transferase (GST3) Lactoylglutathione lyase Multidrug resistance protein 1-like Major vault protein Aldehyde reductase Monoamine oxidase Nuclear factor erythroid 2-related factor 1 cAMP-responsive element-binding protein-like 2 CREB-binding protein Vitellogenin Acetylcholinesterase Glyoxalase I Cu/Zn superoxide dismutase Glutathione peroxidase Glutathione S-transferase 6-Phosphogluconate dehydrogenase Thioredoxin reductase (NADPH) Metallothionein Peroxidasin Cathepsin D Cathepsin L Cathepsin L1-like, transcript variant 2 Cathepsin L1 Complement C1q-like protein 2 Complement C1q-like protein 4 Complement C1q tumor necrosis factor-related protein 3 Complement C1q and tumor necrosis factor-related protein 9B-like Complement C1q subcomponent subunit B Thioester-containing protein-F CR1-3 Lysozyme 3 Lysozyme 2 Immunodominant MHC-associated peptides (SIMP) Lysosomal aspartic protease Lysosomal acid lipase/cholesteryl ester hydrolase Mucin-2 precursor C-type lectin Thioredoxin reductase 3 Fibrillin-2-like

−1.2546 −4.8378 7.4606 / 4.8469 / 1.4924 −1.074 6.1455 1.9666 1.6441 1.4475 1.6893 / −3.0269 / 1.4924 / −1.4455 / 1.4260 1.3005 1.0413 3.2219 / −1.1539 / / 1.3062 1.5639 1.2471 / / 1.1485 6.1577 −2.1983 −1.0071 −2.0231

/ / −4.6892 1.2754 −4.2835 −1.0169 −1.4837 / −5.3772 −1.5564 / −1.0595 / −1.3318 −8.085 1.4390 −1.4837 −3.1520 −4.4590 −1.0169 / / 1.6784 −1.6647 1.2352 1.5317 1.2007 1.3878 1.7276 / / 1.1371 1.2647 −1.3332 −4.3178 3.7512

Metabolism of xenobiotic compounds

t4:5

Oxidative stress

t4:6

Apoptosis

t4:7

Immune system

R

R

E

C

T

E

D

P

t4:4

R O

Q2

405 406

U

N

C

O

The mRNA expression levels of oxidoreductase, dehydrogenase, cytochrome 450 and monoamine oxidase were all significantly changed 407 after exposure to BaP. These genes consisted of a redox complex that 408 belonged to Phase I metabolism. They were involved in the oxidative 409 phosphorylation and drug metabolism which are two of the significant410 ly enriched pathways after scallop exposure to BaP. The change was also 411 observed in the enzymatic activity levels which were reported in the ju412 venile white shrimp Litopenaeus vannamei exposed to 3 μg/L BaP (Ren 413 et al., 2014). Our further qRT-PCR data showed that after three days of 414 BaP exposure increasing of CYP2D10, NQO1, SOD and GPx2 transcripts 415 Q17 were also observed in the rat model (Xia et al., 2011). CYP450 is 416 a major family of enzymes involved in the primary or Phase I metabo417 lism of xenobiotics, including drugs, carcinogens and pesticides 418 (Anzenbacher and Anzenbacherova, 2001). We clearly observed that 419 BaP provoked differential expression of cytochrome P450 genes, includ420 ing CYP2B5, CYP 2D10, CYP 3A2, and CYP 2U1 (Table 3). Among these 421 genes, CYP 2D10 transcriptionally up-regulated 7.5 fold, while CYP 422 Q18 2B5 and CYP2U1 displayed an opposite response at 3 days. Kolaja and 423 Kramer (2002) show that the BaP single exposure not only included in424 creases in CYP1A1 and CYP1B1, but also included increases in non-Ah 425 Q19 regulated genes such as CYP2D10 expression. It was reported that BaP 426 may not probably induce CYP1A via the aromatic hydrocarton receptor 427 (AhR) pathway like PAHs (Cherng et al., 2001; Peters et al., 2004), which

F

t4:3

O

t4:1 t4:2

Y. Cai et al. / Gene xxx (2014) xxx–xxx

−2.0191 / 1.3005 −3.3561

1.2352 1.3171 2.4646 1.9035 / −5.8657

was supported by our data. To our knowledge, previous studies only reported that genes CYP2C9, CYP3A1, and CYP3A4 were involved in the BaP metabolism in human (Ding and Kaminsky, 2003), while the gene CYP3A2 was only linked to the BaP metabolism (Lupp et al., 2001). This study is the first to demonstrate greater abundance of transcripts of a CYP3A2 in scallop exposed to BaP in two points. In addition to greater abundances of transcripts of CYPs, abundances of transcripts of monoamine oxidase (moA) were greater after being exposed for three days, but was less at ten days point (Table 3). It is not known if enzyme products of these transcripts biotransform organic compounds in BaP and what effects this might have on the toxicity of BaP. NAD(P)H:quinone oxidoreductase-1 (NQO1) is a phase II enzyme that catalyzes the two-electron reduction of quinones (Lind, 1985), thereby diminishing the toxicity of quinone-structured compounds (Joseph and Jaiswal, 1994). The expression of NQO1 is highly inducible by BaP (Kann et al., 2005). In our study, NQO1 was also induced during BaP exposure. Biotransformation and detoxification of chemicals are also dependent on Phase II and Phase III reactions. Phase II enzymes include glutathione-S-transferase (GST) that catalyzes conjugation of reduced glutathione (L-γ-glutamyl-L-cysteinyl-glycine, GSH) to reactive electrophiles of endogenous and exogenous compounds. Therefore, GST activity has been widely used as a biomarker of exposure to chemical substances in molluscs (Fitzpatrick et al., 1997; Liu et al., 2014).

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

428 Q20 429 430 431 Q21 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450

Y. Cai et al. / Gene xxx (2014) xxx–xxx Table 5 KEGG classification of C. farreri genes involved in metabolism of xenobiotic pathways.

t5:6 Q3

3 days vs. 10 days

5 (2.86%) 3 (1.90%) 2 (4.17%) 1 (0.97%) 1 (0.86%) 2 (1.33%) 5 (1.43%) 5 (1.25%) 3 (1.71%) 3 (1.90%) 1 (1.25%) 1 (0.47%) 1 (0.43%) 1 (0.57%) 2 (0.54%) 3 (0.86%) 1 (1.57%) 3 (2.89%) 11 (1.43%)

Hoarau et al. (2001) and Hoarau et al. (2002) reported that, in the Med-

T

451

1 (0.66%)

F

10d-DGE pathways

1 (0.97%) 1 (0.86%) 2 (1.75%) 4 (4.21%) 4 (1.14%) 2 (0.50%) 1 (1.25%)

O

t5:5

Drug metabolism Fatty acid metabolism Glycerolipid metabolism Glycolysis/gluconeogenesis Insulin signaling pathway MAPK signaling pathway Metabolism of xenobiotics by cytochrome P450 NOD-like receptor signaling pathway p53 signaling pathway PPAR signaling pathway Citrate cycle (TCA cycle) Drug metabolism Fatty acid metabolism Glycerophospholipid metabolism Insulin signaling pathway MAPK signaling pathway p53 signaling pathway PPAR signaling pathway Metabolism of xenobiotics by cytochrome P450 Apoptosis Peroxisome p53 signaling pathway Lysosome Insulin signaling pathway Steroid hormone biosynthesis Vitamin digestion and absorption Focal adhesion

R O

3d-DGE pathways

Unigenes in each pathway

P

t5:4

Pathway

D

t5:3

E

t5:1 t5:2

7

C

452 iterranean clam Ruditapes decussatus, the expression of the GST seems 453 Q22 to be regulated by inducing BaP. In addition, the GST related mRNA ex-

pression of 3 μg/L BaP in gill tissue of white shrimp L. vannamei was found to start increasing at three days, reaching maximum levels at 456 six days, while returning to the control level during the last day (Ren 457 Q23 et al., 2014). The same trend was obtained in our study, the expression 458 of GST3 was induced at three days but less (−1.02 folds) in the diges459 tive gland of scallops exposed to BaP for ten days (Table 3). Substrates 460 of Phases I and II metabolism are excreted from cells by multidrug resis461 tance associated protein (MRP) (Borst et al., 1999). Abundance of tran462 scripts of MRP1 was less in digestive gland from scallops exposed to BaP 463 (Table 3). Global cellular stress response has been defined as all proteins 464 over-produced due to environmental stress. Following BaP exposure, 465 several protein families including cAMP responsive element binding 466 protein (CREB), CREB binding protein (CBP), nuclear factor erythroid 467 2-related factor 1 (Nrf1), HSP70, CYP families, GST, MRP and GPx were 468 observed and over transcribed in 3 days of exposure or 10 days of 469 Q24 exposure of C. farreri. The present study found that CREB at 3 days of 470 exposure was over transcribed 1.68 fold change (q-value 3.00E-03) 471 and CBP was found significantly down-regulated (−1.33 fold change, 472 q-value 2.13E-05) at 10 days of exposure.

U

N C O

R

R

E

454 455

473

4.2. Oxidative Stress

Fig. 3. The expression level of differently expressed genes validated by qPCR. A: The results of genes related to metabolism of BaP after 3 days of exposure; B: the results of genes related to metabolism of BaP after 3 days of exposure; C: the relative quantitative expression of all the related genes after 3 days and 10 days of exposure.

474 475

In order to survive in its environment, cells have developed a sophisticated system of enzymes and small molecules that limit the formation of highly reactive oxygen species and effectively detoxify them. This system includes superoxide dismutases, glutathione peroxidases, catalase, thioredoxins, and peroxiredoxins in various cell organelles. In addition, water-soluble antioxidants such as glutathione and lipid-soluble antioxidants such as tocopherols and carotenoids are critical cofactors or direct scavengers of various oxidants. Another strategy to limit oxidant-induced tissue damage is to keep the levels of unbound

transition metals in the cell extremely low through effective chelation by storage and transport proteins such as ferritin, transferrin, and ceruloplasmin (Jaeschke, 2010). In this research, after BaP exposure we have observed mostly antioxidant-related genes, which changed in a regular pattern and identified at different time points. Exposure to BaP might have caused oxidative stress in the digestive glands of scallops. A minimal concentration of reactive oxygen species (ROS) that is required for normal cellular function is maintained by antioxidant defense enzymes and a pool of antioxidants. Concentrations of

476 477 478 479 480 481 482

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

483 484 485 486 487 488 489 490 491

Y. Cai et al. / Gene xxx (2014) xxx–xxx

556 557

4.4. RNAseq vs qPCR

566

F

organisms exposed to BaP. There were tissue and dose- and timedependent changes in expression of mRNA of the pro-inflammatory genes in scallop exposed to BaP (Zhang et al., 2009). In general, BaP attenuated the effects on expression of genes of the complement system. This study is the first to report effects of BaP on the complement system. In conclusion, this study demonstrated that exposure of scallops to BaP leads to significant inhibition of immunity functions in the initial phase, and then the immunocompetence has a direct enhancement during the exposure, which supports the view that a relationship exists between BaP and immunomodulation in molluscs.

C

E

R

R

O

N

C

542

4.3. Immune System

U

541

The immune response mechanisms of polycyclic aromatic hydrocarbons (with benzopyrene as representative) on bivalve have been reported. The immune defense system of bivalve molluscs includes 545 cellular immunity and humoral immunity, and circulating hemocytes 546 play an important role (Zhang et al., 2005). Hemocytes can kill or clear 547 bacteria through phagocytosis. Hemocytes can also produce various hu548 moral factors, such as lysozyme, antimicrobial peptides and so on, 549 Q28 which are involved in humoral immunity (Liu and Mai, 2002). Abun550 dances of transcripts of a novel protein similar to vertebrate comple551 ment component complement C1q like protein 2, complement C1q 552 like protein 4, complement C1q tumor necrosis factor-related protein 553 3, and complement C1q and tumor necrosis factor-related protein 9B554 like were greater in digestive glands from scallops exposed to BaP 555 (Table 3). Studies have reported impaired immune responses in 543 544

558 559 560 561 562 563 564 565

567

5. Conclusion

596

In conclusion, we presented an extensive survey of BaP responsive genes, showing differential expression in scallop at the different time. These genes are involved in many metabolic processes in response to BaP. Our results not only demonstrate transcriptional complexity in scallop with BaP, but also open up a possibility for identification of the genes in regulating bivalve tolerance to xenobiotic stress. Furthermore, the high-throughput sequencing platform described here serves as a powerful approach for developing biomarkers monitoring molecular response to the organic hazards.

597

6. Uncited references

606 Q31

P

R O

O

A subset of transcripts that was determined to be of greater abundance by RNAseq was quantified by qPCR. Both methods report abundances of transcripts as fold-change relative to the control but the calculation and normalization methods are not the same. The RNA-seq expression values were calculated on the basis of RPKM (Mortazavi et al., 2008), and qPCR fold-change values were calculated by use of the MNE method and incorporated reference gene (Simon, 2003). For transcripts quantified by both methods, directions of change were always the same and the magnitude of the fold-change in abundance was very similar (Fig. 3A,B). The data support previous studies where RNA-seq was an effective method for determining changes in abundances of transcripts (Marioni et al., 2008). One major advantage of the high throughput sequencing is the direct estimation of gene expression. Moreover, millions of generated scallop sequences allowed us to analyze enrichment of various genes. The identified genes within the libraries displayed varied abundance from one another. According to the GO and KEGG distributions, the majority of genes were enriched to involve metabolic processes (e.g. genes for metabolism of carbohydrates, organic acids, amino acids, and secondary metabolites) and regulation of cellular process (Table 4). These results indicate that scallop exposure to BaP altered many basic biological pathways associated with metabolism of carbohydrates, proteins, lipids, and nucleic acids. We also retrieved many important genes grouped into the three categories including response to intracellular and environmental stimulus, glutathione transferase activity, and oxidoreductase activity. Using the database categories, we identified 194/161 (3 days/10 days) genes for stimulus-responsive protein and 64/48 (3 days/10 days) genes for oxidoreductases, all of which were differentially regulated by BaP.

T

GSH, the most abundant antioxidant in eukaryotic cells, are maintained in part by the actions of glutathione synthase (GS), glutathione reduc494 tase (GR), and glutathione peroxidase (GPx) (Circu and Aw, 2010). 495 The transcripts of glyoxalase1 and GPx were differentially expressed 496 in digestive glands from scallops exposed to BaP at three and ten days 497 (Table 3). When GSH is utilized to reduce ROS a molecule of NADPH 498 supplied from the pentose phosphate pathway is used to reduce the 499 resulting disulfide-oxidized GSH. Abundances of transcripts of enzymes 500 of the pentose phosphate pathway, including 6-phosphate dehydroge501 nase (6-pgdh) was greater in digestive glands of scallops exposed to 502 BaP three days (Table 3). In mammals, expression of 6-PGDH and GPx 503 Q25 is regulated by the Nrf2 signaling pathway (Thimmulappa et al., 2006) 504 that is activated by metabolites from CYP catalyzed metabolism of or505 ganic compounds. These metabolites might have originated from me506 tabolism of BaP by CYP enzymes. These results are consistent with a 507 study where concentrations of ROS were greater in the hemolymph of 508 the temperate scallop Pecten maximus after PAH phenanthrene expo509 sure (Hannam et al., 2010). 510 Oxidative stress can disrupt redox conditions required for proper 511 folding of proteins. Consequently, the glutaredoxin and thioredoxin sys512 tems, which are important for proper folding of proteins, are important 513 in the response to oxidative stress (Berndt et al., 2008). The thioredoxin 514 system is comprised of NADPH, thioredoxin (Trx), thioredoxin reduc515 tase (TrxR), and protein disulfide isomerases (PDI) (Wilkinson and 516 Gilbert, 2004). Abundances of transcripts of TrxR were greater in diges517 tive glands of scallops exposed to BaP (Table 3). In addition to the −1.4518 fold and −4.45 less abundance of transcripts of glutathione peroxidase, 519 the function of Grx5 is unknown, but loss of function in yeast leads to 520 oxidative damage and sensitization to ROS (Rodrıguez-Manzaneque 521 et al., 2002). 522 Metallothioneins (MTs) are low-molecular weight, cysteine rich 523 stress proteins. The presence of cysteine in MTs confers substantial 524 metal binding capabilities. MTs play a role in the metabolism of the rel525 Q26 atively non-toxic essential metals (zinc and copper) (Dunn et al., 1987) 526 as well as in the detoxification of toxic metals such as cadmium 527 (Viarengo and Nott, 1993). During oxidative stress, synthesis of MTs 528 may increase several-fold (Thornalley and Vašák, 1985) to protect the 529 cells from cytotoxicity (Aschner et al., 1998) and DNA damage (Cai 530 et al., 1995). The BaP treatment could act on the metallothionein level 531 by stimulating ROS (reactive oxygen species) production (increased 532 MFO (mixed function oxidase) activity). As a matter of fact, many xeno533 biotics (such as BaP) may produce oxyradicals by inducing cytochrome 534 P450 monooxygenase in relation with NADH cytochrome b5 reductase 535 a stress syndrome in the animal and subsequently increase MT levels. 536 Q27 Bremner (1986) hypothesized that MTs may act as free radical scaven537 gers, the complex system for regulation of MT synthesis in animals sub538 jected to a wide variety of stress strongly suggested that it forms part of 539 a host-defense mechanism and is of some functional significance in cel540 lular metabolism (Roméo et al., 1997).

D

492 493

E

8

English and Storey, 2003 Haritash and Kaushik, 2009 Hayashi et al., 2003 Jaiswal, 2004 Maher et al., 2008 Wang et al., 2007

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

568 569 570 571 572 573 574 575 576 577 578 Q29 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 Q30 595

598 599 600 601 602 603 604 605

607 608 609 610 611 612

Y. Cai et al. / Gene xxx (2014) xxx–xxx

Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biol. 11, R106. Anzenbacher, P., Anzenbacherova, E., 2001. Cytochromes P450 and metabolism of xenobiotics. Cell Mol. Life Sci. 58, 737–747. Aschner, M., Conklin, D.R., Yao, C.P., Allen, J.W., Tan, K.H., 1998. Induction of astrocyte metallothioneins (MTs) by zinc confers resistance against the acute cytotoxic effects of methylmercury on cell swelling, Na+ uptake, and K+ release. Brain Res. 813, 254–261. Baussant, T., Bechmann, R.K., Taban, I.C., Larsen, B.K., Tandberg, A.H., Bjørnstad, A., Torgrimsen, S., Nævdal, A., Øysæd, K.B., Jonsson, G., Sanni, S., 2009. Enzymatic and cellular responses in relation to body burden of PAHs in bivalve molluscs: a case study with chronic levels of North Sea and Barents Sea dispersed oil. Mar. Pollut. Bull. 58, 1796–1807. Berndt, C., Lillig, C.H., Holmgren, A., 2008. Thioredoxins and glutaredoxins as facilitators of protein folding. Biochimica et Biophysica Acta (BBA)-Molecular. Cell Res. 1783, 641–650. Borst, P., Evers, R., Kool, M., Wijnholds, J., 1999. The multidrug resistance protein family. Biochim. Biophys. Acta Biomembr. 1461, 347–357. Bremner, I., 1986. Nutritional and physiological significance of metallothionein. Experientia Suppl. 52, 81–107. Brown, M., Davies, I.M., Moffat, C.F., Craft, J.A., 2006. Application of SSH and a macroarray to investigate altered gene expression in Mytilus edulis in response to exposure to benzo[a]pyrene. Mar. Environ. Res. 62, S128–S135. Cai, L., Koropatnick, J., Cherian, M.G., 1995. Metallothionein protects DNA from copperinduced but not iron-induced cleavage in vitro. Chem. Biol. Interact. 96, 143–155. Cherng, S.H., Lin, P., Yang, J.L., Hsu, S.L., Lee, H., 2001. Benzo[g, h, i]perylene synergistically transactivates benzo[a]pyrene-induced CYP1A1 gene expression by aryl hydrocarbon receptor pathway. Toxicol. Appl. Pharmacol. 170, 63–68. Circu, M.L., Aw, T.Y., 2010. Reactive oxygen species, cellular redox systems, and apoptosis. Free Radic. Biol. Med. 48, 749–762. Ding, X., Kaminsky, L.S., 2003. Human Extrahepatic Cytochromes P450: function in xenobiotic metabolism and tissue-selective chemical toxicity in the respiratory and gastrointestinal tracts. Annu. Rev. Pharmacol. Toxicol. 43, 149–173. Dunn, M.A., Blalock, T.L., Cousins, R.J., 1987. Metallothionein. Exp. Biol. Med. 185, 107–119. English, T.E., Storey, K.B., 2003. Freezing and anoxia stresses induce expression of metallothionein in the foot muscle and hepatopancreas of the marine gastropod Littorina littorea. J. Exp. Biol. 206, 2517–2524. Faria, M., Carrasco, L., Diez, S., Riva, M.C., Bayona, J.M., Barata, C., 2009. Multibiomarker responses in the freshwater mussel Dreissena polymorpha exposed to polychlorobiphenyls and metals. Comp. Biochem. Physiol. C Toxicol. Pharmacol. 149, 281–288. Figueras, A., Costa, M.M., Novoa, B., 2012. Applications of functional genomics in molluscs aquaculture. Funct. Genom. Aquacult. 377. Fitzpatrick, P.J., O'Halloran, J., Sheehan, D., Walsh, A.R., 1997. Assessment of a glutathione S-transferase and related proteins in the gill and digestive gland of Mytilus edulis (L.), as potential organic pollution biomarkers. Biomarkers 2, 51–56. Frouin, H., Pellerin, J., Fournier, M., Pelletier, E., Richard, P., Pichaud, N., Rouleau, C., Garnerot, F., 2007. Physiological effects of polycyclic aromatic hydrocarbons on soft-shell clam Mya arenaria. Aquat. Toxicol. 82, 120–134. Geffard, O., Geffard, A., His, E., Budzinski, H., 2003. Assessment of the bioavailability and toxicity of sediment-associated polycyclic aromatic hydrocarbons and heavy metals applied to Crassostrea gigas embryos and larvae. Mar. Pollut. Bull. 46, 481–490. Goldberg, E.D., Bowen, V.T., Farrington, J.W., Harvey, G., Martin, J.H., Parker, P.L., Robert, W.R., William, R., Eric, S., Gamble, E., 1978. The mussel watch. Environ. Conserv. 5, 101–125. Gomes, T., Pereira, C.G., Cardoso, C., Bebianno, M.J., 2013. Differential protein expression in mussels Mytilus galloprovincialis exposed to nano and ionic Ag. Aquat. Toxicol. 136, 79–90. Hannam, M.L., Bamber, S.D., Galloway, T.S., John Moody, A., Jones, M.B., 2010. Effects of the model PAH phenanthrene on immune function and oxidative stress in the haemolymph of the temperate scallop Pecten maximus. Chemosphere 78, 779–784. Haritash, A.K., Kaushik, C.P., 2009. Biodegradation aspects of polycyclic aromatic hydrocarbons (PAHs): a review. J. Hazard. Mater. 169, 1–15. Hayashi, A., Suzuki, H., Itoh, K., Yamamoto, M., Sugiyama, Y., 2003. Transcription factor Nrf2 is required for the constitutive and inducible expression of multidrug resistance-associated protein1 in mouse embryo fibroblasts. Biochem. Biophys. Res. Commun. 310, 824–829. Hoarau, P., Gnassia-Barelli, M., Romeo, M., Girard, J.P., 2001. Differential induction of glutathione S-transferases in the clam Ruditapes decussatus exposed to organic compounds. Environ. Toxicol. Chem. 20, 523–529.

F

621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691

O

References

C

E

R

R

N C O

U

618

R O

620

616 617

P

619

This work was supported by National Natural Sciences Foundation of China (research grant no. 30972237) and the State Ocean Administration Special Public Project of China (research grant no. 201105013). We thank the staff at the Laboratory of Environmental Physiology of Aquatic Animals for their help with sampling and taking care of the scallops.

D

614 615

Hoarau, P., Garello, G., Gnassia-Barelli, M., Romeo, M., Girard, J.P., 2002. Purification and partial characterization of seven glutathione S-transferase isoforms from the clam Ruditapes decussatus. Eur. J. Biochem. 269, 4359–4366. Jaeschke, H., 2010. Antioxidant defense mechanisms, In: McQueen, Charlene A. (Ed.), Comprehensive Toxicology, 2nd ed. Elsevier, Oxford, pp. 319–337. Jaiswal, A.K., 2004. Nrf2 signaling in coordinated activation of antioxidant gene expression. Free Radic. Biol. Med. 36, 1199–1207. Joseph, P., Jaiswal, A.K., 1994. NAD (P) H: quinone oxidoreductase1 (DT diaphorase) specifically prevents the formation of benzo[a]pyrene quinone-DNA adducts generated by cytochrome P4501A1 and P450 reductase. Proc. Natl. Acad. Sci. 91, 8413–8417. Kann, S., Huang, M.Y., Estes, C., Reichard, J.F., Sartor, M.A., Xia, Y., Puga, A., 2005. Arseniteinduced aryl hydrocarbon receptor nuclear translocation results in additive induction of phase I genes and synergistic induction of phase II genes. Mol. Pharmacol. 68, 336–346. Kim, J.Y., Kwon, J., Kim, J.E., Koh, W.S., Chung, M.K., Yoon, S., Song, C.W., Lee, M., 2005. Identification of potential biomarkers of genotoxicity and carcinogenicity in L5178Y mouse lymphoma cells by cDNA microarray analysis. Environ. Mol. Mutagen. 45, 80–89. Kolaja, K.L., Kramer, J.A., 2002. Toxicogenomics: an opportunity to optimise drug development and safety evaluation. Expert Opin. Drug Saf. 1, 275–286. Langmead, B., Trapnell, C., Pop, M., Salzberg, S.L., 2009. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinforma. 12, 323. Lind, C., 1985. Formation of benzo[a]pyrene-3,6-quinol mono-and diglucuronides in rat liver microsomes. Arch. Biochem. Biophys. 240, 226–235. Liu, S., Mai, K., 2002. The progress of studies on molluscs immunological system and mechanism—a review. Acta Oceanol. Sin. 25, 95–105. Liu, N., Pan, L., Wang, J., Yang, H., Liu, D., 2012. Application of the biomarker responses in scallop Chlamys farreri to assess metals and PAHs pollution in Jiaozhou Bay, China. Mar. Environ. Res. 80, 38–45. Liu, D., Pan, L., Li, Z., Cai, Y.F., Miao, J., 2014. Metabolites analysis, metabolic enzyme activities and bioaccumulation in the clam Ruditapes philippinarum exposed to benzo[a] pyrene. Ecotoxicol. Environ. Saf. 107, 251–259. Lupp, A., Danz, M., Müller, D., 2001. Morphology and cytochrome P450 isoforms expression in precision-cut rat liver slices. Toxicology 16, 53–66. Maher, J.M., Aleksunes, L.M., Dieter, M.Z., Tanaka, Y., Peters, J.M., Manautou, J.E., Klaassen, C.D., 2008. Nrf2- and PPARα-mediated regulation of hepatic Mrp transporters after exposure to perfluorooctanoic acid and perfluorodecanoic acid. Toxicol. Sci. 106, 319–328. Mao, X., Cai, T., Olyarchuk, J.G., Wei, L., 2005. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. Mathijs, K., Brauers, K.J., Jennen, D.G., Boorsma, A., van Herwijnen, M.H., Gottschalk, R.W., Kleinjans, J.C., van Delft, J.H., 2009. Discrimination for genotoxic and nongenotoxic carcinogens by gene expression profiling in primary mouse hepatocytes improves with exposure time. Toxicology 112, 374–384. McHale, C.M., Zhang, L., Hubbard, A.E., Smith, M.T., 2010. Toxicogenomic profiling of chemically exposed humans in risk assessment. Mutat. Res. 705, 172–183. Menzie, C.A., Potocki, B.B., Santodonato, J., 1992. Exposure to carcinogenic PAHs in the environment. Environ. Sci. Technol. 26, 1278–1284. Meyer, R.P., Podvinec, M., Meyer, U.A., 2002. Cytochrome P450 CYP1A1 accumulates in the cytosol of kidney and brain and is activated by heme. Mol. Pharmacol. 62, 1061–1067. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. Pan, L.Q., Miao, J.J., Wang, J., Liu, J., 2008. AHH activity, tissue dose and DNA damage in different tissues of the scallop Chlamys farreri exposed to benzo[a]pyrene. Environ. Pollut. 153, 192–198. Paules, R.S., Aubrecht, J., Corvi, R., Garthoff, B., Kleinjans, J.C., 2011. Moving forward in human cancer risk assessment. Environ. Health Perspect. 119, 739–743. Peters, A.K., Van Londen, K., Bergman, Å., Bohonowych, J., Denison, M.S., Van den Berg, M., Sanderson, J.T., 2004. Effects of polybrominateddiphenyl ethers on basal and TCDDinduced ethoxyresorufin activity and cytochrome P450-1A1 expression in MCF-7, HepG2, and H4IIE cells. Toxicology 82, 488–496. Qiu, Y.W., Zhou, J.L., Maskaoui, K., Hong, H.S., Wang, Z.D., 2004. Distribution of polycyclic aromatic hydrocarbons in water and sediments from Daya Bay and their ecological hazard assessment. J. Trop. Oceanogr. 23, 72–80. Ren, X., Pan, L., Wang, L., 2014. Metabolic enzyme activities, metabolism-related genes expression and bioaccumulation in juvenile white shrimp Litopenaeus vannamei exposed to benzo[a]pyrene. Ecotoxicol. Environ. Saf. 104, 79–86. Rodrıguez-Manzaneque, M.T., Tamarit, J., Bellı́, G., Ros, J., Herrero, E., 2002. Grx5 is a mitochondrial glutaredoxin required for the activity of iron/sulfur enzymes. Mol. Biol. Cell 13, 1109–1121. Roméo, M., Cosson, R.P., Gnassia-Barelli, M., Risso, C., Stien, X., Lafaurie, M., 1997. Metallothionein determination in the liver of the sea bass Dicentrarchus labrax treated with copper and B(a)P. Mar. Environ. Res. 44, 275–284. Schirmer, K., Dixon, D.G., Greenberg, B.M., Bols, N.C., 1998. Ability of 16 priority PAHs to be directly cytotoxic to a cell line from the rainbow trout gill. Toxicology 127, 129–141. Simon, P., 2003. Q-Gene: processing quantitative real-time RT-PCR data. Bioinformatics 19, 1439–1440. Sun, J., Nishiyama, T., Shimizu, K., Kadota, K., 2013. TCC: an R package for comparing tag count data with robust normalization strategies. BMC Bioinforma. 14, 219. Thimmulappa, R.K., Lee, H., Rangasamy, T., Reddy, S.P., Yamamoto, M., Kensler, T.W., Biswal, S., 2006. Nrf2 is a critical regulator of the innate immune response and survival during experimental sepsis. J. Clin. Investig. 116, 984–995.

E

Acknowledgments

T

613

9

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777

10

transcriptomic analysis with Yesso scallop (Patinopecten yessoensis). PLoS One 8, e63927. Wilkinson, B., Gilbert, H.F., 2004. Protein disulfide isomerase. Biochim. Biophys. Acta Protein Proteomics 1699, 35–44. Xia, X., Yan, J., Shen, Y., Tang, K., Yin, J., Zhang, Y., Yang, D., Liang, H., Ye, J., Weng, J., 2011. Berberine improves glucose metabolism in diabetic rats by inhibition of hepatic gluconeogenesis. PLoS One 6, e16556. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Method gene ontology analysis for RNA-seq: accounting for selection bias. PMC free article. , (PubMed). Yue, D.D., Wang, L.M., 2012. Relationship mariculture shellfish production and carbon sequestration in China. J. Jiangsu Agric. Sci. 40, 246–248. Zhang, W.Z., Wu, X.Z., Li, D.F., Sun, J.F., 2005. Immunological functions of blood cells in the scallop Chlamys farreri. Acta Zool. Sin. 51, 669–677. Zhang, L., Pan, L., Liu, J., 2009. Immunotoxicity effect of benzo[α]pyrene on scallop Chlamys farreri. J. Ocean Univ. China 8, 89–94. Zhao, X., Yu, H., Kong, L., Li, Q., 2012. Transcriptomic responses to salinity stress in the pacific oyster Crassostrea gigas. PLoS One 7, e46244.

N

C

O

R

R

E

C

T

E

D

P

R O

O

F

Thornalley, P.J., Vašák, M., 1985. Possible role for metallothionein in protection against radiation-induced oxidative stress. Kinetics and mechanism of its reaction with superoxide and hydroxyl radicals. Biochim. Biophys. Acta Protein Struct. Mol. Enzymol. 827, 36–44. Tian, S., Pan, L., Sun, X., 2013. An investigation of endocrine disrupting effects and toxic mechanisms modulated by benzo[a]pyrene in female scallop Chlamys farreri. Aquat. Toxicol. 144, 162–171. Viarengo, A., Nott, J.A., 1993. Mechanisms of heavy metal cation homeostasis in marine invertebrates. Comp. Biochem. Physiol. C Comp. Pharmacol. 104, 355–372. Wang, X., Tomso, D.J., Chorley, B.N., Cho, H.Y., Cheung, V.G., Kleeberger, S.R., Bell, D.A., 2007. Identification of polymorphic antioxidant response elements in the human genome. Hum. Mol. Genet. 16, 1188–1200. Wang, L., Pan, L., Liu, N., Liu, D., Xu, C., Miao, J., 2011. Biomarkers and bioaccumulation of clam Ruditapes philippinarum in response to combined cadmium and benzo[α]pyrene exposure. Food Chem. Toxicol. 49, 3407–3417. Wang, S., Hou, R., Bao, Z., Du, H., He, Y., Su, H., Li, Y., Du, H., Hu, J., Wang, S., Hu, X., 2013. Transcriptome sequencing of Zhikong scallop (Chlamys farreri) and comparative

U

778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794

Y. Cai et al. / Gene xxx (2014) xxx–xxx

Please cite this article as: Cai, Y., et al., Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.003

795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811