- Email: [email protected]
RNA-seq based on transcriptome reveals differ genetic expressing in Chlamys farreri exposed to carcinogen PAHs
Accepted Manuscript Title: RNA-Seq based on transcriptome reveals differ genetic expressing in Chlamys farreri exposed to carcinogen PAHs Author: Qian Jin Luqing Pan Tong Liu Fengxiao Hu PII: DOI: Reference:
S1382-6689(14)00289-0 http://dx.doi.org/doi:10.1016/j.etap.2014.11.019 ENVTOX 2136
To appear in:
Environmental Toxicology and Pharmacology
Received date: Revised date: Accepted date:
20-8-2014 23-11-2014 26-11-2014
Please cite this article as: Jin, Q., Pan, L., Liu, T., Hu, F.,RNA-Seq based on transcriptome reveals differ genetic expressing in Chlamys farreri exposed to carcinogen PAHs, Environmental Toxicology and Pharmacology (2014), http://dx.doi.org/10.1016/j.etap.2014.11.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Highlights The tested scallops were exposed to a mixture of four carcinogen polycyclic aromatic hydrocarbons (BaP, BaA, BbF and CHR).
high-throughput RNA-sequencing (RNA-seq) technologies.
ip t
209 and 260 genes were identified as significantly up- or down-regulated, respectively, by
The HSP70 mRNA expression was inhibited when scallop exposed to the mixture of the
cr
four PAHs.
CYP450 3A2 to monitor complex marine environment..
us
We suggested CYP450 2P1 mRNA expression as a better candidate biomarker than
an
QO mRNA expression can be a potential marker for prediction of the biological effects of a
Ac ce p
te
d
M
mixture of PAHs on scallop.
1
Page 1 of 24
A manuscript for Environmental Toxicology and Pharmacology
Title: RNA-Seq based on transcriptome reveals differ genetic
Authors: Qian Jin a, Luqing Pan a,*, Tong Liu a, Fengxiao Hu a a
ip t
expressing in Chlamys farreri exposed to carcinogen PAHs Key Laboratory of Mariculture, Ministry of Education, Ocean University of China,
cr
Qingdao 266003, PR China
Tel.: +86 532 82032963; Fax: +86 532 82032963;
an
E-mail address: [email protected] (L.Q. Pan).
us
*Corresponding author: Luqing Pan, Professor, Ph.D.
Address: Fisheries College, Ocean University of China, Yushan Road 5, Qingdao
Ac ce p
te
d
M
266003 China
2
Page 2 of 24
Abstract The effects of a mixture of carcinogen polycyclic aromatic hydrocarbons (BaP, BaA, BbF and CHR) on transcriptional responses in the digestive gland of scallop, Chlamys farreri, were investigated by high-throughput RNA-sequencing (RNA-seq)
ip t
technologies. In total, 209 and 260 genes were identified as significantly up- or
down-regulated, respectively. Functional analysis based on gene ontology (GO)
cr
classification system and the Kyoto Encyclopedia of Genes and Genomes (KEGG)
database revealed that PAHs significantly altered the expression of genes involved in
us
stress response, detoxication, antioxidation which were extensively discussed. In particular, CYP450 2P1 and QO mRNA expression were found to be up-regulated by
an
exposure to PAHs mixture, suggesting that CYP450 2P1 and QO mRNA expression can be a potential marker for prediction of the biological effects of a mixture of PAHs
M
on scallops.
Keywords: Carcinogen PAHs; Digital gene expression; Differential expression;
Ac ce p
te
d
Chlamys farreri
3
Page 3 of 24
1. Introduction Polycyclic Aromatic Hydrocarbons (PAHs) are persistent hydrophobic organic pollutants ubiquitously found in the environment and they originate mainly from
ip t
incomplete combustion of organic materials and fossil fuels. As a kind of global persistent organic pollutants (POPs), PAHs have also been included in the group of
cr
persistent toxic substances (PTS), and 16 unsubstituted PAHs have been listed as “Priority Pollutants” by the United States Environment Protection Agency (US EPA),
us
due to their well-known carcinogenic and mutagenic properties (Chung et al., 2007). Concentrations of PAHs in seawater and sediment are up to 1717.87 ng/L and 1606.89
an
ng/g dry weight, respectively, in the Bohai Sea, suggesting that PAHs pollution may be a serious problem in China (Men et al., 2009).
specific
PAHs,
namely,
M
The European Food Safety Authority (EFSA) introduced a system of four benzo[a]anthracene
(BaA),
benzo[a]pyrene
(BaP),
d
benzo[b]fluoranthene (BbF) and chrysene (CHR), assessing that the sum of the four PAH compounds was the most suitable indicator for PAHs in food (EFSA, 2008). In
te
recent years, PAH pollution has been much aggravated by the development of
Ac ce p
offshore oil and marine traffic, making the concentrations of BaP, BbF BaA and CHR in the marine environment increasingly higher. According to the survey, the above three PAH contents in surface sediment from Zhanjiang Bay, China, were up to 19.90, 36.06 and 39.65 ng/g dry weight, respectively (Huang et al., 2012). BaP is considered as a model substance and many single PAH toxicity studies have been carried out on bivalves to determine the adverse effects of PAHs (Pan et al., 2008, Liu et al., 2014). Yet there are few reports on the toxicity of other kinds of PAHs, such as BaA, BbF and CHR, and little work has been done on co-exposure of the four PAHs. In most case, cytochrome P450 (CYP450), glutathione-S-transferase (GST), antioxidant enzyme activities and DNA damage can be induced or depressed in a dose-dependent manner by PAHs and can therefore be used as biomarkers to effectively assess the biological effects of PAHs (Anzenbacher et al., 2001, Kovacevic et al., 2008, Michel et al., 2013). Among the different types of biomarkers, molecular 4
Page 4 of 24
biomarkers have the potential to provide rapid, cost-effective and sensitive diagnostic tools to assess the physiological impacts of toxicants on organisms. In addition, gene transcription data can provide mechanistic insights into the responses of organisms to toxicants (Connon et al., 2008). With the advancement of molecular biology
ip t
technology, differential gene expression methods are being used more frequently to select better molecular biomarker options (Brown et al., 2006). High-throughput
cr
RNA-sequencing (RNA-seq) technologies, including Roche/454 pyrosequencing, Illumina/Solexa sequencing technology, and Applied bio-systems SOLiD sequencing
us
technology, have led to a revolution in the fileds of transcriptomics and genome characterization in recent years (Martin and Wang, 2011). The Illumina sequencing
an
platform, the HiSeqTM 2000 sequencing system, with relatively low cost, can generate up to two billion paired-end reads with an average length of 100 bp at 99.5%
M
accuracy per run. Digital gene expression (DGE) analysis based on this sequencing platform had been applied to many aquatic organisms, such as Oryzias latipes (Oh et
d
al., 2012), Crassostrea gigas (Zhao et al., 2012) and Chlamys farreri (Fu et al., 2014), in the research of their response to environmental stresses.
te
The scallop Chlamys farreri, one of economic mollusk species, is mainly
Ac ce p
distributed in coastal areas of North China, Korea and Japan, and the scallop is usually cultured in the shallow seas and considered as suitable bio-indicator for biomonitoring studies because of their sessile nature, filter-feeding habits, and pollutant bioconcentration (Goldberg et al., 1978). In this study, aiming to develop a further understanding of the molecular responses of scallops upon exposure to PAHs mixture, Digital Gene Expression (DGE) analysis was conducted by employing Illumina HiSeqTM 2000 sequencing system. Quantitative real-time PCR (qRT-PCR) was performed to verify several selected differentially expressed genes and expound the molecular mechanism coping with PAHs exposure stress. Furthermore, genes associated with defensing PAHs stress identified in this attempt may be potential biomarkers for PAHs pollution monitoring in aquatic environment. 2. Materials and methods 2.1. Animals and PAHs exposure experiment 5
Page 5 of 24
The experimental scallop C. farreri (shell length 53.7 ± 1.5 mm, shell height 57.2 ± 1.8 mm, shell width 18.5 ± 0.8 mm, mass 20.7 ± 1.7 g) were obtained from Shazikou (Yellow Sea, Qingdao, China) on June 12, 2013, and acclimated for one weeks prior to exposure in laboratory condition with aerated seawater at a temperature
ip t
of 22℃ temperature and a salinity of 31‰. During the acclimation period, scallops were fed twice daily with dry powder of spirulina platensis. Seawater was refreshed
cr
every 24 hours.
PAHs used for the experiment were purchased from Dr. Ehrenstorfer GmbH
us
(Augsburg, Germany). A minimal concentration of DMSO (< 0.001%) was used as a vehicle to prevent cellular damage to the scallop. For this experiment, animals were
an
randomly divided into two experimental groups with three replicates: (1) control; (2) 2 μg/L PAHs for a ten days’ exposure experiment. All other conditions were kept the
M
same as those used during acclimation. The artificial mixture of the 4 PAHs was prepared based on the relative concentration ratio of 4 PAHs in the coastal seawater in
d
China (BaP: BbF: BaA: CHR = 5:3:1:1).
2.2. Sampling and total RNA extraction
te
Digestive glands from nine scallops (three scallops per replicate) in the treated
Ac ce p
group and control group were sampled, then were frozen in liquid nitrogen and stored in 1 mL Trizol reagent (Invitrogen, USA) at ˗ 80 °C. Scallops (both the treated
scallops and control ones) for DGE were sampled at 10th day while 0, 3, 6 and 10 d for RT-qPCR.
Total RNA was extracted following the instructions of the Trizol reagent. RNA
degradation and contamination were monitored in 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer
(Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). 2.3 Library preparation and DGE sequencing Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep 6
Page 6 of 24
Kit for Illumina (NEB, USA) following manufacturer’s recommendations. Briefly, the mRNA was purified from a total amount of 3 g RNA per sample using poly-T oligo-attached magnetic beads and then fragmented into small pieces using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction
ip t
Buffer (5X). The fragments were then used for synthesizing the first- and second-strand cDNA. After end repair and single nucleotide A (adenine) addition,
cr
adaptors with a hairpin loop structure were ligated to the DNA fragments. Afterwards, suitable fragments (150-200bp in length) were chosen as templates for polymerase
us
chain reaction (PCR) amplification. 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 95 °C before
an
PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified
M
(AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Library preparations were sequence on an Illumina Hiseq 2000 platform
d
and 100 bp pair-end reads were generated. The clean Illumina sequencing data have been submitted to the GenBank (SRA) with the accession number SRR1303115.
te
2.4 Mapping of digital gene expression (DGE) reads
Ac ce p
The reference unigenes were assembled from transcriptome (SRA accession SRP018044), and further annotated using Non-redundant protein (Nr) database, Nucleotide collection (Nt) database, Swiss-Plot, Gene Ontology database (GO), Clusters of Orthologous Groups (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008; Ye et al., 2006). After removing reads with adaptors, reads containing poly N and low quality reads from raw data, we obtained clean reads. Bowtie v0.12.9 was used to align the high-quality reads to the reference unigene sequences. Numbers of the reads mapped to each unigene were determined using RESM (Li et al., 2011). RPKM (reads per kilobase of exon model per million mapped reads) of each gene was calculated, considering the effects of sequencing depth and gene length for the reads count at the same time. 2.4 Quality control Raw data (raw reads) of fastq format were firstly processed through in-house perl 7
Page 7 of 24
scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated. All the downstream analyses were based on the clean data with high quality.
ip t
2.5 Analysis of differential expression analysis
The differential expression detection of genes across samples was performed
cr
using the DEGSeq R package (1.12.0). The resulting p-values were fitted using the
Benjamini and Hochberg’s approach. Corrected p value of 0.05 and log2 (fold change)
us
of 1 were selected as the threshold for significantly differential expression. 2.6 GO and KEGG
an
Gene Ontology (GO) enrichment analysis of differential expressed genes was conducted by the GOseq package, in which gene length was taken into consideration.
M
A corrected p value less than 0.05 was considered as a threshold for significant enrichment of the gene sets. The topGO software was used to visualize enriched GO
d
categories.
Pathway enrichment analysis can further identify significantly enriched
te
metabolic pathways or signal transduction pathways based on the KEGG database
Ac ce p
resource. KOBAS (2.0) software was utilized to test the statistical enrichment of differential expression genes in KEGG pathways. Pathways with a corrected p value ≤ 0.05 were designated as significantly enriched pathways in DGEs. 2.7 Quantitative real time PCR validation Ten genes were selected for the confirmation of DGE data by Quantitative
real-time PCR (qRT-PCR) and18s was used as a reference gene. 1 µg of total RNA from each sample was used to synthesize the first strand cDNA using the PrimeScriptTM RT Reagent Kit (TaKaRa). The qRT-PCR amplifications were
performed with a Thermo PIKOREAL 96 Real-Time PCR system using SYBR Green Supermix (TaKaRa) according to the manufacturer’s instructions. The qRT-PCR program was 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s, and 52 °C for 20 s for annealing, and then 72 °C for 30 s for extension. The relative expression ratio of a target gene is calculated as described in the 2−ΔΔCt method (Livak and Schmittgen, 8
Page 8 of 24
2001), and expressed in comparison to a reference gene. Statistical analysis was carried out with the statistical software SPSS 19.0. All quantitative data were expressed as means ± S.D. 3. Result
ip t
3.1. Analysis of DGE libraries
Two DGE libraries were composed of digestive glands of the control and
cr
PAHs-exposed using Hiseq 2000 sequencing system (Illumina). Total 12,711,432 and 11,630,310 raw reads were generated in the control and PAHs-exposed DGE libraries,
us
respectively. A total of 12,571,017 and 11,500,871 clean reads were left after removing reads with adaptors, reads containing poly N and low quality reads from
an
raw data. Values of the Q20 percentage and GC percentage of the clean reads in two DGE libraries were 97.40% and 97.46%, 43.87% and 42.92%, respectively (Table 1).
M
So the two DGE libraries were reliable. 3.2. Analysis of mapping
d
We mapped the sequences of the two DGE libraries to the reference transcriptome database of C.farreri, which contains 92999 unigenes. In control and
te
PAHs-treated DGE libraries, 77.07% and 78.60% of the clean reads were mapped to a
Ac ce p
gene in the reference database (Table 1). The number of clean reads mapped to a gene was normalized to RPKM values to
accurately measure the gene expression level (Fig.1). The results showed that mRNAs transcribed from the major types of genes were represented in low RPKM values, and only a small proportion of genes was highly expressed. 3.3. GO analysis and KEGG pathways A total of 469 differentially expressed genes were detected, of which 209
were up-regulated and 260 were down-regulated. To understand the functions of these differentially expressed genes, 151 differentially expressed genes were mapped to terms in GO database and compared with the transcriptome background. The top 30 most highly enriched GO terms were presented in Fig. 2. GO has three ontologies, namely, molecular function, cellular component, and biological process, and genes related to the terms “catalytic activity”, “oxidoreductase activity”, “extracellular 9
Page 9 of 24
space” were dominant in the PAHs-exposed scallop. Interestingly, some genes related to the establishment of meiotic spindle localization, norepinephrine transmembrane transporter activity were all down-regulated in PAHs-exposed digestive gland of C. farreri.
ip t
KEGG assignments were used to classify functional annotations of the identified
genes to further understand the biological functions of the genes. Among the
cr
differentially expressed genes, 613 were assigned to 212 pathways in the KEGG database in 4008 DGE class. 19 terms were significantly enriched (Q value ≤ 0.05). 10
us
significantly enriched KEGG pathways (Corrected p < value) were listed in Table 2. Especially, “Glutathione metabolism” pathway related to detoxification metabolism
an
was significantly enriched. In addition, most of the genes related to pathways “Lysosome”, “Fat digestion and absorption” and “Glutathione metabolism” were
M
up-regulated, while the majority of genes related to pathways “ECM-receptor interaction”, “Betalain biosynthesis” were down-regulated.
d
3.4. Identification and verification of differentially expressed genes The expression profile of 7 transcripts of C. farreri exposed to 0.02, 0.2 and 2
te
μg/L of PAHs was performed by qRT-PCR at 0, 1, 6 and 10 days. The relative mRNA
Ac ce p
expression level of 6 genes including 4 up-regulated, 2 down-regulated genes from DGE libraries was detected by qRT-PCR. These chosen genes list in Table 3. were based on the relation with stress response, detoxication, antioxidation processes that could be important for the exposure event. The relative expression levels of the target genes were calculated with the 2−ΔΔCt
method (Livak and Schmittgen, 2001) (Fig.3). The data revealed that CYP450 3A2, CYP 450 2P1, quinone oxidoreductase (QO) and GST were extremely significantly up-regulated after 1 day exposure at all three concentrations, while HSP70 and Cu/Zn-SOD down-regulated. The HSP70 mRNA expression remained low and significantly different from control values until the end of the experiment (p < 0.01).
But the differences among exposure days were non-significant (p > 0.05). The levels of CYP2P1, CYP3A2 and QO were induced by PAHs at first 1 day then raised gradually by exposure period prolonged. The different between control and treated 10
Page 10 of 24
groups was extremely significant (p < 0.01). GST gene expression in 0.02 and 0.2 μg/L PAHs group peaked at day 6 and remained constant until the end of the experiment, while expression in 2 μg/L PAHs group increased during the whole experiment. The differences of exposed groups were non-significant (p > 0.05). The
the differences among concentrations were remarkable (p < 0.01).
ip t
Cu/Zn-SOD expression PAHs exposure group reached the lowest level at day 10 and
cr
The values of expression were imported to Microsoft Excel and data analyses
were performed using IBM SPSS statistical 19.0 (SPSS Inc.). All data were presented
us
as means ± S.D. (n=3) of the relative mRNA expression. 4. Discussion
an
At present, most researches of differentially expressed genes were focused on single PAH exposure, but less on a mixture of PAHs (Oh et al., 2012, van Delft et al.,
M
2012). High-throughput RNA-sequencing (RNA-seq) technologies, including Roche/454 pyrosequencing, Illumina/Solexa sequencing technology have led to a
d
revolution in the fields of transcriptomics and genome characterization in recent years (Martin and Wang, 2011). Previously, Oh et al. (2012) observed that CYP450 2P1
te
mRNA in liver of medaka was significantly induced by synergistic effects of cPAHs
Ac ce p
mixture, but not by single BaP. On the contrary, that study found that genes encoding malate dehydrogenase, anti-thrombin III, and transferrin were significantly up-regulated in response to exposure to BaP, but not in response to the cPAHs mixture, suggesting that exposure to cPAHs results in a reduction of the effect of BaP. Considering the complex chemical interactions between BaP and other three PAHs, we compared the DGE library with another one conducted by our colleagues (Cai et al., 2014). Under the same environment condition, single BaP exposure (1μg/L)
induced 477 differentially expressed genes in digestive gland of Chlamys farreri at day 10. GST mRNA was down-regulated under single BaP exposure with 1.0169 fold change, while GST was up-regulated under PAHs mixture exposure with 1.1036 fold change. CYP450 3A2 was up-regulated in both single and mix PAH exposure, but with different fold changes, 1.2754 and 1.6721. Besides that, mRNA expression of HSP70, CYP450 2P1 and quinone oxidoreductase were not significantly altered due 11
Page 11 of 24
to single BaP, but to PAHs mixture. In comparison with 469 differentially expressed genes identified in our attempt, the result indicated that impact of PAHs mixture on gene expression could not be identified as synergy or antagonism respectively. Considering vast contaminations in ocean ecosystem, researches about effect of PAHs
ip t
mixture in gene expression are meaningful, but needed further efforts. Furthermore, as presented in Fig.3, qRT-PCR results demonstrated the reliability and accuracy of DGE
cr
data.
Molecular biomarkers can be used to investigate the impact or the physiological
us
response of organism to aquatic pollution. Differentially expressed genes in the digestive gland of C.farreri exposed to PAHs mixture under experimental conditions
closely related with detoxification (Table 3).
an
were investigated. Some genes were identified from the DGE library, of which they
M
Heat shock proteins (HSPs) have been identified in several bivalve species in response to various chemical, physical and pathogenic stresses and appear to represent
d
a general marker of non-specific stress (Qiu et al., 2007). According to the apparent molecular mass, HSPs are classified into several families, including HSP90 (83–90
te
kDa), HSP70 (66–78 kDa), HSP60, HSP47, and small heat shock protein (sHSP;
Ac ce p
15–43 kDa) (Morimoto, 1993). All organisms from bacteria to mammals exposed to environmental stressors respond by synthetizing heat-shock protein (HSP70) (Schlesinger et al., 1982). HSP70 gene expression is implicated in protein repair, transport and protection from oxidative stress (Contardo-Jara et al., 2010). In zebra mussels, an exposure to BaP induced a persistent increase in HSP70 mRNA level after 28 days (Châtel et al., 2012). Similarly, in a field study, Copat et al. (2012) found that in Engraulis encrasicolus tissues collected from the Catania Gulf containing 16 PAHs, HSP70 protein was weakly expressed. Interestingly, our result showed that gene HSP70 expression was remarkably inhibited after one day exposure, but approximately invariant at day 6 and 10 in the digestive gland of scallop, even the highest mRNA relative expression was 0.6 fold compared to control samples. It was induced that gene HSP70 stabilization is likely a strategy for PAHs mixture tolerance in C.farreri, which accordance with previous research in Tar Pond killifish (Christine 12
Page 12 of 24
et al., 2009). However, in the present study, HSP70 mRNA expression was not correlated with the amount of PAHs present in the different stations measured previously in the laboratory (Michel et al., 2013), which was parallel to our data. This may suggest that other parameters such as temperature and pH might also affect
ip t
HSP70 mRNA level induction, as demonstrated in other field studies (Hamer et al.,
2004). The fact is that HSP proteins participate in several important processes such as
cr
folding of the nascent polypeptide chain and protection of proteins from denaturation or aggregation, the suppression of HSPs by PAHs mixture may result in a reduction in
us
protein stability and consequent accumulation of abnormal proteins, leading to changes in cell metabolism and apoptosis.
an
Previous studies have indicated the digestive gland was the main target organ for detoxification of environmental pollutant and identification of related biomarkers.
M
Among differentially expressed genes, we identified several genes related to detoxification. CYP450 is a major family of enzymes involved in the primary or
d
phase I metabolism of xenobiotics, including drugs, carcinogens and pesticides (Anzenbacher and Anzenbacherova, 2001). We clearly observed that PAHs mixture
te
provoked differential expression of cytochrome P450 genes, including CYP2P1,
Ac ce p
CYP3A2. Oh et al. (2012) discovered that CYP450 2P1 mRNA was significantly up-regulated in response to exposure to cPAHs, but not BaP alone. In addition, the gene CYP3A2 was only linked to the benzo[a]pyrene metabolism (Lupp et al., 2001). Our qRT-PCR results revealed that CYP2P1 mRNA expression (4.6 fold change) was more severely changeable than CYP3A2 (1.3 fold change) when scallops coped with 10 day PAHs mixture stress. It may suggest that CYP450 2P1 mRNA can be induced by synergistic effects of cPAHs mixture, while the synergistic effects less worked on CYP450 3A2. Moreover, CYP450 2P1 could be regarded a better candidate biomarker than CYP450 3A2 to monitor complex marine environment. It is proved that there were many kinds of BaP phase I metabolites in mammals, such as BaP epoxides, phenols, dihydrodiols, tetrahydrotetrol, quinones and diolepoxides. Furthermore, BaP appears to be metabolised predominantly to quinones in mussel rather than diols as in vertebrates (Sjölin and Livingstone, 1997), and 13
Page 13 of 24
relative to mechanisms of toxicity, quinones undergo redox cycling with production of O2•− making oxidative stress a prominent feature of PAH toxicity in mussel. These processes cause interference with electron transport, either inducing enzyme systems such as cytochrome P450 that mediate oxidation reactions or depleting protective
ip t
antioxidants. Fortunately, enzymatic (e.g., SOD, GST and GPx) and non-enzymatic
(e.g., GSH, ascorbate, carotenoids, and free proline) antioxidant defense systems in
cr
living organisms can provide cellular protection against oxidative stress. Otherwise, as the pollutants are metabolized, reactive molecules are produced and interact with
us
the genetic material and cause DNA damage as a toxic terminal of the induction of a cascade of cellular events. DNA integrity reduction represents the genotoxicity which
an
may have a long-term effect on the sustainability of a particular population. In the present study, quinone oxidoreductase (QO) mRNA expression was up-regulated
M
remarkable, suggesting that scallops were able to mount a successful adaptive response towards the threats posed by PAHs. Because QO is a principal phase II
d
enzyme that plays an important role in the complete reduction and detoxification of highly reactive quinines. In addition to metabolizing foreign compounds, QO also has
te
antioxidant functions, such as protecting against lipid peroxidationand defending
Ac ce p
against intracellular oxidative stress by scavenging superoxide (Siegel et al., 2004, Talalay and Dinkova-Kostova, 2004). Thus, our results indicated that PAHs mixture exposure may disturb the redox balance of the tested tissues. To further examine the functional distribution of these differentially expressed
genes and classify the functional annotations of them, we applied KEGG analysis on these identified genes. Higher corrected p-value take part in glutathione metabolism
(ko: 00480) were found in our database. Glutathione (GSH) is reported to act as an effective antioxidant in marine animal system and as a reactant in conjugation with electrophilic substances (Kovacevic et al., 2008). Glutathione S-transferase (GST) is a dominated antioxidant enzyme, conjugating electrophilic compounds with GSH. Previous field studies have shown a positive correlation between PAHs burden and GST activity (Gowland et al., 2002). Although GSTs mRNA expression levels could be induced in BaP treated groups (Liu et al., 2014), the GST gene was not 14
Page 14 of 24
differentially regulated as expected in the DGE library. PAHs exposure experiment was carried out to ascertain how GST worked during the 10 days. In the present research, variation of GST mRNA expression was slight in all treated columns, but the difference with the control group was significant (p < 0.05). The result inferred GST
ip t
expression was induced by individual aryl hydrocarbon receptor agonists (BaP), but weakened by PAHs mixture containing non-agonists (Garner et al., 2012). What’s
cr
more, the differences among concentration gradients were positive in our result. Xu et al. (2010) also observed a good linear correlation between the BaP exposed
us
concentrations and the GSTs expression levels at 10 d in the digestive gland of scallops (p < 0.01). So GSTs have been recognized as biomarkers of exposure to
an
oxidative stress-inducing chemical pollutants including PAHs mixture (Boutet et al., 2004).
M
Another key gene belong to the antioxidant defense system is superoxide dismutase (SOD), which can eliminate reactive oxygen species (ROS) produced by
d
phase I reactions. ROS can be responsible for lipid peroxidation processes and eventually result in DNA damage. Wang et al. (2011) reported that Cu/Zn-SOD
te
mRNA expression levels in the digestive gland of R. philippinarum were induced by
Ac ce p
0.01 mg/L BaP. However, in the present study, Cu/Zn-SOD mRNA expression was down-regulated obvious even in the 0.02μg/L PAHs group (p < 0.05), indicating suppression of Cu/Zn-SOD mRNA expression caused by PAHs mixture. But Liu et al. (2014) hold the opinion that Mn-SOD appears to be a better biomarker for BaP exposure than Cu/Zn-SOD. While, in the present study, all three treatments inhibited the SOD mRNA expression significantly at the whole exposure time, and the expression decreased along with mixture concentration increased, indicating significant dose-time effect. Thus, Cu/Zn-SOD could reflect the PAHs pollution at low concentration. 5. Conclusion Our study examined the transcriptome response of C.farreri to PAHs mixture stress using digital gene expression analysis. The results demonstrated that PAHs mixture significantly altered expression of genes implicated in stress response, 15
Page 15 of 24
detoxication, antioxidation and so on. As to assess multi-PAHs pollution based on scallops C.farreri, CYP450 2P1 could offer more alert information than CYP3A2. Besides, regarding QO as PAHs pollution biomarker applied in scallop was rarely reported. Moreover, PAHs mixture inhibits the expression of a key gene involved in
ip t
glutathione metabolism and may disturb the detoxication system of C.farreri, which would put forward new evidence about making GST as biomarker in bivalve. The
cr
current results provide a strong basis for future research on molecular mechanism response to PAHs exposure to bivalve. More validation experiments are also needed
us
to confirm if these differentially expressed genes identified in this study can serve as
an
potential biomarkers for PAHs pollution.
Acknowledgements
M
This work was supported by the National Natural Sciences Foundation of China (30972237) and the State Ocean Administration Special Public Project of China
d
(201105013). We thank the staff at the Laboratory of Environmental Physiology of Aquatic Animal for the help with sampling and taking care of the scallops.
te
Conflict of interest statement
Ac ce p
The authors have nothing to disclose. References
Anzenbacher, P., Anzenbacherova, E., 2001. Cytochromes P450 and metabolism of xenobiotics. CMLS, Cell. Mol. Life Sci. 58, 737-747.
Boutet, I., Tanguy, A., Moraga, D., 2004. Characterisation and expression of four mRNA sequences encoding glutathione S-transferases pi, mu, omega and sigma classes in the Pacific oyster Crassostrea gigas exposed to hydrocarbons and pesticides. Mar. Biol., 146, 53–64 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, 128-135.
Cai, Y.F., Pan, L.Q., Hu, F.X., Jin, Q., Liu, T., 2014. Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene. 16
Page 16 of 24
GENE. 551, 261-270. Châtel, A., Faucet-Marquis, V., Perret, M., Gourlay-Francé, C., Uher, E., Pfohl-Leszkowicz, A., Vincent-Hubert, F., 2012. Genotoxicity assessment and detoxification induction in Dreissena polymorpha exposed to benzo[a]pyrene
ip t
MUTAGENESIS. 27, 703–711
Christine Paetzold, S., Ross, N. W., Richards, R. C., Jones, M., Hellou, J., & Bard, S.
cr
M., 2009. Up-regulation of hepatic ABCC2, ABCG2, CYP1A1 and GST in
multixenobiotic-resistant killifish (Fundulus heteroclitus) from the Sydney Tar
us
Ponds, Nova Scotia, Canada. Mar. Environ. Res. 68, 37-47.
Chung, M., Hu, R., Cheung, K., Wong, M., 2007. Pollutants in Hong Kong soils:
an
polycyclic aromatic hydrocarbons. CHEMOSPHERE. 67, 464–473. Connon, R., Hooper, H.L., Sibly, R.M., Lim, F.L., Heckmann, L.H., Moore, D.J.,
M
Watanabe, H., Soetaert, A., Cook, K., Maund, S.J., Hutchinson, T.H., Moggs, J., Coen, W.D., Iguchi, T., Callaghan, A., 2008. Linking molecular and population
42, 2181-2188.
d
stress responses in Daphnia magna exposed to cadmium. Environ. Sci. Technol.
te
Contardo-Jara, V., Pflugmacher, S., Nutzmann, G., Kloas, W., Wiegand C., 2010. The
Ac ce p
beta receptor blocker metoprolol alters detoxification processes in the non-target organism Dreissena polymorpha. Environ. Pollut., 158, 2059–2066
Copat, C., Brundo, M.V., Arena, G., Grasso, A., Oliveri Conti, G., Ledda, C., Fallico, R., Sciacca, S., Ferrante, M., 2012. Seasonal variation of bioaccumulation in Engraulis encrasicolus (Linneaus, 1758) and related biomarkers of exposure.
Ecotox. Envieon. Safe. 86, 31-37.
EFSA, 2008. Scientific Opinion of the Panel on Contaminants in the Food Chain on a request from the European Commission on Polycyclic Aromatic Hydrocarbons in Food. EFSA Journal, 724, 1–114 Fu, X., Sun, Y., Wang, J., Xing, Q., Zou, J., Li, R.,Wang, Z., Wang, S., Hu, X., Zhang, L., Bao, Z., 2014. Sequencing-based gene network analysis provides a core set of gene resource for understanding thermal adaptation in Zhikong scallop Chlamys farreri. Mol. Ecol. Resour. 14, 184-198. 17
Page 17 of 24
Garner, L. V., & Di Giulio, R. T., 2012. Glutathione transferase pi class 2 (GSTp2) protects against the cardiac deformities caused by exposure to PAHs but not PCB-126 in zebrafish embryos. Comp. Biochem. Phys. C. 155, 573-579. Goldberg, E.D., Bowen, V.T., Farrington, J.W., Harvey, G., Martin, J.H., Parker, P.L.,
ip t
Risebrough, R.W., Robertson, W., Schneider, E., Gamble, E., 1978. The mussel watch. Enviro.Cnserv. 5, 101-125.
cr
Gowland, B.T.G., Mcintosh, A.D., Davies, I.M., Moffat, C.F., Webster, L., 2002.
Implications from a field study regarding the relationship between polycyclic
us
aromatic hydrocarbons and glutathione S-transferase activity in mussels. Mar. Environ. Res. 54, 231-235.
an
Hamer, B., Hamer, D.P., Muller, W.E., Batel, R., 2004. Stress-70 proteins in marine mussel Mytilus galloprovincialis as biomarkers of environmental pollution: a
M
field study. Environ. Int. 30, 873–882.
Huang, W., Wang, Z., Yan, W., 2012. Distribution and sources of polycyclic aromatic
d
hydrocarbons (PAHs) in sediments from Zhanjiang Bay and Leizhou Bay, South China. Mar. Pollut. Bull. 64, 1962–1969.
te
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., Katayama, T.,
Ac ce p
Kawashima, S., Okuda, S., Tokimatsu, T.,Yamanishi, Y., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, 480-484.
Kovacevic, T.B., Borokovic, S.S., Pavlovic, S.Z., Despotovic, S.G., Saicic, Z.S., 2008. Glutathione as a suitable biomarker in hepatopancreas, gills and muscles of three fresh water cryfish species. Arch. Biol. Sci. 60, 59–66.
Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC BIOINFORMATICS. 12, 323.
Liu, D., Pan, L., Cai, Y., Li, Z., & Miao, J., 2014. Response of detoxification gene mRNA expression and selection of molecular biomarkers in the clam Ruditapes philippinarum exposed to benzo[a]pyrene. Environ. Pollut. 189, 1-8. Livak, K.J., Schmittgen, T.D. 2001. Analysis of relative gene expression data using real-time
quantitative
PCR
and
the
2(T)(−delta
delta
C)
method.
METHODS .25,402-408. 18
Page 18 of 24
Lupp, A., Danz, M., Müller, D., 2001. Morphology and cytochrome P450 isoforms expression in precision-cut rat liver slices. TOXICOLOGY. 161, 53-66. Martin, J. A., Wang, Z., 2011. Next-generation transcriptome assembly. Nat. Rev. Genet. 12, 671-682..
ip t
Men, B., He, M., Tan, L., Lin, C., Quan, X., 2009. Distributions of polycyclic
aromatic hydrocarbons in the Daliao River Estuary of Liaodong Bay, Bohai Sea
cr
(China). Mar. Pollut. Bull. 58, 818-826.
Michel, C., Bourgeault, A., Gourlay-Francé, C., Palais, F., Geffard, A.,
us
Vincent-Hubert, F., 2013. Seasonal and PAH impact on DNA strand-break levels in gills of transplanted zebra mussels. Ecotoxicol. Environ. Saf. 92, 18–26.
an
Morimoto, R.I., 1993. Cell in stress: transcriptional activation of heat shock genes. SCIENCE. 259, 1409–1410, at doi: http://dx.doi.org/10.1126/science 8451637.
M
Oh, J.H., Moon, H.B., Choe, E.S., 2012. Alterations in differentially expressed genes by exposure to a mixture of carcinogenic polycyclic aromatic hydrocarbons in
d
the liver of Oryzias latipes. Environ. Toxicol. Phar. 33, 403-407. Pan, L., Ren, J., Liu, J., 2005. Effects of benzo[k]fluoranthene exposure on the
te
biomarkers of scallop Chlamys farreri. Comp. Biochem. Physiol. C: Toxicol.
Ac ce p
Pharmacol. 141, 248–256.
Qiu, L.M., Song, L.S., Yu, Y.D., Xu, W., Ni, D.J., Zhang, Q.C., 2007. Identification and characterization of a myeloid differentiation factor 88 (MyD88) cDNA from Zhikong scallop Chlamys farreri. Fish Shellfish Immun. 23, 614–623.
Schlesinger, M.J., Ashburner, M., Tissieres, A., 1982. Heat-Shock from Bacteria to Men 1982 Cold Spring Harbor Laboratory, New York, p. 286
Siegel, D., Gustafson, D.L., Dehn, D.L., Han, J.Y., Boonchoong, P., Berliner, L.J., Ross, D., 2004. NAD(P)H: quinone oxidoreductase 1: role as a superoxide scavenger. Mol. Pharmocol. 65, 1238-1247. Sjölin, A. M., Livingstone, D. R., 1997. Redox cycling of aromatic hydrocarbon quinones catalysed by digestive gland microsomes of the common mussel (Mytilus edulis L.). Aquat. Toxicol. 38, 83-99. Talalay, P., Dinkova-Kostova, A.T., 2004. Role of nicotinamide quinone 19
Page 19 of 24
oxidoreductase 1 (NQO1) in protection against toxicity of electrophiles and reactive oxygen intermediates. Method Enzymol. 382, 355-364. Wang, L., Pan, L.Q., Liu, N., Liu, D., Xu, C.Q., Miao, J.J., 2011. Biomarkers and bio-accumulation of clam Ruditapes philippinarum in response to combined
ip t
cadmium and benzo[a]pyrene exposure. Food Chem. Toxicol. 49, 3407-3417.
Xu, C., Pan, L., Liu, N., Wang, L., & Miao, J., 2010. Cloning, characterization and
cr
tissue distribution of a pi-class glutathione S-transferase from clam (Venerupis
philippinarum): Response to benzo[a]pyrene exposure. Comp. Biochem. Phys. C.
us
152, 160-166.
Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., Wang, J., Li, S., Li, R.,
an
Bolund, L.,Wang, J., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34, 293- 297.
M
Zhao, X.L., Yu, H., Kong, L.F., Li, Q., 2012. Transcriptomic responses to salinity
Ac ce p
te
d
stress in the pacific oyster Crassostrea gigas. PloS one. 7, e46244.
20
Page 20 of 24
Table 1. Control
PAHs-exposed
12711432
11630310
Clean reads
12571017
11500871
0.03
0.03
Q20 (%)
97.46
97.40
Total reads
12571017
11500871
Total mapped (%)
77.07
78.60
an
us
cr
Error rate (%)
ip t
Category Raw reads
Summary for the C.farreri DGE in control and PAHs-exposed libraries
M
Table 2. The 10 most significantly enriched KEGG pathway (Corrected p value < 0.05) Pathway
Up-/down-regulated sequences
Corrected p value
ko04974
Protein digestion and absorption
7/9
2.99E-05
ko04512
ECM-receptor interaction
5/12
0.0015
ko04640
Hematopoietic cell lineage
5/3
0.0025
ko04142
Lysosome
15/2
0.0040
te
d
Pathway ID
Tyrosine metabolism
3/6
0.0040
ko04975
Fat digestion and absorption
7/0
0.0040
ko04614
Renin-angiotensin system
4/1
0.0040
ko00965
Betalain biosynthesis
1/3
0.0040
ko04145
Phagosome
9/8
0.0050
ko00480
Glutathione metabolism
5/2
0.0236
Ac ce p
ko00350
Table 3. Genes selected for Real-time qRT-PCR. Gene number
Description
log2.Fold_changea
CL6151.Contig1_zhikongshanbei
Hsp70,putative [Ixodes scapularis]
-2.4552
CL8118.Contig1_zhikongshanbei
Cytochrome P450 2P1[Fundulus heteroclitus]
1.2745
Unigene64457_zhikongshanbei
Cytochrome P450 3A2 [Rattus norvegicus]
1.6127
CL4845.Contig2_zhikongshanbei
sulfide:quinone oxidoreductase [Arenicola marina]
2.6718
CL10249.Contig1_zhikongshanbei
Glutathione S-transferase
1.1036
CL1826.Contig1_zhikongshanbei
Superoxide dismutase [Cu-Zn] [Brugia pahangi]
-3.2591
a
Up- or down-regulated in DGE of PAHs exposed clams C.farreri (10 d).
21
Page 21 of 24
us
cr
ip t
Fig. 1
Ac ce p
te
d
M
an
Fig.2
22
Page 22 of 24
te
d
M
an
us
cr
ip t
Fig.3
Fig.1 Length distribution of clean reads mapped.
Ac ce p
Fig.2 Histogram presentation of Gene Ontology classification. The results indicated the most enriched GO terms (PAHs exposed vs Control). Terms were summarized in three main categories: biological process, cellular component and molecular function. * indicated the enriched GO term. Fig.3 Relative expression by Quantitative real-time PCR of 6 selected genes. The Ct values of 0d were used as baseline with the value 1. The relative expression levels of the genes were calculated using 2
−ΔΔCt
model. ** Significantly different from the control group (p < 0.01).
23
Page 23 of 24
Ac ce p
te
d
M
an
us
cr
ip t
Conflict of interest statement The authors have nothing to disclose.
24
Page 24 of 24
S1382-6689(14)00289-0 http://dx.doi.org/doi:10.1016/j.etap.2014.11.019 ENVTOX 2136
To appear in:
Environmental Toxicology and Pharmacology
Received date: Revised date: Accepted date:
20-8-2014 23-11-2014 26-11-2014
Please cite this article as: Jin, Q., Pan, L., Liu, T., Hu, F.,RNA-Seq based on transcriptome reveals differ genetic expressing in Chlamys farreri exposed to carcinogen PAHs, Environmental Toxicology and Pharmacology (2014), http://dx.doi.org/10.1016/j.etap.2014.11.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Highlights The tested scallops were exposed to a mixture of four carcinogen polycyclic aromatic hydrocarbons (BaP, BaA, BbF and CHR).
high-throughput RNA-sequencing (RNA-seq) technologies.
ip t
209 and 260 genes were identified as significantly up- or down-regulated, respectively, by
The HSP70 mRNA expression was inhibited when scallop exposed to the mixture of the
cr
four PAHs.
CYP450 3A2 to monitor complex marine environment..
us
We suggested CYP450 2P1 mRNA expression as a better candidate biomarker than
an
QO mRNA expression can be a potential marker for prediction of the biological effects of a
Ac ce p
te
d
M
mixture of PAHs on scallop.
1
Page 1 of 24
A manuscript for Environmental Toxicology and Pharmacology
Title: RNA-Seq based on transcriptome reveals differ genetic
Authors: Qian Jin a, Luqing Pan a,*, Tong Liu a, Fengxiao Hu a a
ip t
expressing in Chlamys farreri exposed to carcinogen PAHs Key Laboratory of Mariculture, Ministry of Education, Ocean University of China,
cr
Qingdao 266003, PR China
Tel.: +86 532 82032963; Fax: +86 532 82032963;
an
E-mail address: [email protected] (L.Q. Pan).
us
*Corresponding author: Luqing Pan, Professor, Ph.D.
Address: Fisheries College, Ocean University of China, Yushan Road 5, Qingdao
Ac ce p
te
d
M
266003 China
2
Page 2 of 24
Abstract The effects of a mixture of carcinogen polycyclic aromatic hydrocarbons (BaP, BaA, BbF and CHR) on transcriptional responses in the digestive gland of scallop, Chlamys farreri, were investigated by high-throughput RNA-sequencing (RNA-seq)
ip t
technologies. In total, 209 and 260 genes were identified as significantly up- or
down-regulated, respectively. Functional analysis based on gene ontology (GO)
cr
classification system and the Kyoto Encyclopedia of Genes and Genomes (KEGG)
database revealed that PAHs significantly altered the expression of genes involved in
us
stress response, detoxication, antioxidation which were extensively discussed. In particular, CYP450 2P1 and QO mRNA expression were found to be up-regulated by
an
exposure to PAHs mixture, suggesting that CYP450 2P1 and QO mRNA expression can be a potential marker for prediction of the biological effects of a mixture of PAHs
M
on scallops.
Keywords: Carcinogen PAHs; Digital gene expression; Differential expression;
Ac ce p
te
d
Chlamys farreri
3
Page 3 of 24
1. Introduction Polycyclic Aromatic Hydrocarbons (PAHs) are persistent hydrophobic organic pollutants ubiquitously found in the environment and they originate mainly from
ip t
incomplete combustion of organic materials and fossil fuels. As a kind of global persistent organic pollutants (POPs), PAHs have also been included in the group of
cr
persistent toxic substances (PTS), and 16 unsubstituted PAHs have been listed as “Priority Pollutants” by the United States Environment Protection Agency (US EPA),
us
due to their well-known carcinogenic and mutagenic properties (Chung et al., 2007). Concentrations of PAHs in seawater and sediment are up to 1717.87 ng/L and 1606.89
an
ng/g dry weight, respectively, in the Bohai Sea, suggesting that PAHs pollution may be a serious problem in China (Men et al., 2009).
specific
PAHs,
namely,
M
The European Food Safety Authority (EFSA) introduced a system of four benzo[a]anthracene
(BaA),
benzo[a]pyrene
(BaP),
d
benzo[b]fluoranthene (BbF) and chrysene (CHR), assessing that the sum of the four PAH compounds was the most suitable indicator for PAHs in food (EFSA, 2008). In
te
recent years, PAH pollution has been much aggravated by the development of
Ac ce p
offshore oil and marine traffic, making the concentrations of BaP, BbF BaA and CHR in the marine environment increasingly higher. According to the survey, the above three PAH contents in surface sediment from Zhanjiang Bay, China, were up to 19.90, 36.06 and 39.65 ng/g dry weight, respectively (Huang et al., 2012). BaP is considered as a model substance and many single PAH toxicity studies have been carried out on bivalves to determine the adverse effects of PAHs (Pan et al., 2008, Liu et al., 2014). Yet there are few reports on the toxicity of other kinds of PAHs, such as BaA, BbF and CHR, and little work has been done on co-exposure of the four PAHs. In most case, cytochrome P450 (CYP450), glutathione-S-transferase (GST), antioxidant enzyme activities and DNA damage can be induced or depressed in a dose-dependent manner by PAHs and can therefore be used as biomarkers to effectively assess the biological effects of PAHs (Anzenbacher et al., 2001, Kovacevic et al., 2008, Michel et al., 2013). Among the different types of biomarkers, molecular 4
Page 4 of 24
biomarkers have the potential to provide rapid, cost-effective and sensitive diagnostic tools to assess the physiological impacts of toxicants on organisms. In addition, gene transcription data can provide mechanistic insights into the responses of organisms to toxicants (Connon et al., 2008). With the advancement of molecular biology
ip t
technology, differential gene expression methods are being used more frequently to select better molecular biomarker options (Brown et al., 2006). High-throughput
cr
RNA-sequencing (RNA-seq) technologies, including Roche/454 pyrosequencing, Illumina/Solexa sequencing technology, and Applied bio-systems SOLiD sequencing
us
technology, have led to a revolution in the fileds of transcriptomics and genome characterization in recent years (Martin and Wang, 2011). The Illumina sequencing
an
platform, the HiSeqTM 2000 sequencing system, with relatively low cost, can generate up to two billion paired-end reads with an average length of 100 bp at 99.5%
M
accuracy per run. Digital gene expression (DGE) analysis based on this sequencing platform had been applied to many aquatic organisms, such as Oryzias latipes (Oh et
d
al., 2012), Crassostrea gigas (Zhao et al., 2012) and Chlamys farreri (Fu et al., 2014), in the research of their response to environmental stresses.
te
The scallop Chlamys farreri, one of economic mollusk species, is mainly
Ac ce p
distributed in coastal areas of North China, Korea and Japan, and the scallop is usually cultured in the shallow seas and considered as suitable bio-indicator for biomonitoring studies because of their sessile nature, filter-feeding habits, and pollutant bioconcentration (Goldberg et al., 1978). In this study, aiming to develop a further understanding of the molecular responses of scallops upon exposure to PAHs mixture, Digital Gene Expression (DGE) analysis was conducted by employing Illumina HiSeqTM 2000 sequencing system. Quantitative real-time PCR (qRT-PCR) was performed to verify several selected differentially expressed genes and expound the molecular mechanism coping with PAHs exposure stress. Furthermore, genes associated with defensing PAHs stress identified in this attempt may be potential biomarkers for PAHs pollution monitoring in aquatic environment. 2. Materials and methods 2.1. Animals and PAHs exposure experiment 5
Page 5 of 24
The experimental scallop C. farreri (shell length 53.7 ± 1.5 mm, shell height 57.2 ± 1.8 mm, shell width 18.5 ± 0.8 mm, mass 20.7 ± 1.7 g) were obtained from Shazikou (Yellow Sea, Qingdao, China) on June 12, 2013, and acclimated for one weeks prior to exposure in laboratory condition with aerated seawater at a temperature
ip t
of 22℃ temperature and a salinity of 31‰. During the acclimation period, scallops were fed twice daily with dry powder of spirulina platensis. Seawater was refreshed
cr
every 24 hours.
PAHs used for the experiment were purchased from Dr. Ehrenstorfer GmbH
us
(Augsburg, Germany). A minimal concentration of DMSO (< 0.001%) was used as a vehicle to prevent cellular damage to the scallop. For this experiment, animals were
an
randomly divided into two experimental groups with three replicates: (1) control; (2) 2 μg/L PAHs for a ten days’ exposure experiment. All other conditions were kept the
M
same as those used during acclimation. The artificial mixture of the 4 PAHs was prepared based on the relative concentration ratio of 4 PAHs in the coastal seawater in
d
China (BaP: BbF: BaA: CHR = 5:3:1:1).
2.2. Sampling and total RNA extraction
te
Digestive glands from nine scallops (three scallops per replicate) in the treated
Ac ce p
group and control group were sampled, then were frozen in liquid nitrogen and stored in 1 mL Trizol reagent (Invitrogen, USA) at ˗ 80 °C. Scallops (both the treated
scallops and control ones) for DGE were sampled at 10th day while 0, 3, 6 and 10 d for RT-qPCR.
Total RNA was extracted following the instructions of the Trizol reagent. RNA
degradation and contamination were monitored in 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer
(Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). 2.3 Library preparation and DGE sequencing Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep 6
Page 6 of 24
Kit for Illumina (NEB, USA) following manufacturer’s recommendations. Briefly, the mRNA was purified from a total amount of 3 g RNA per sample using poly-T oligo-attached magnetic beads and then fragmented into small pieces using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction
ip t
Buffer (5X). The fragments were then used for synthesizing the first- and second-strand cDNA. After end repair and single nucleotide A (adenine) addition,
cr
adaptors with a hairpin loop structure were ligated to the DNA fragments. Afterwards, suitable fragments (150-200bp in length) were chosen as templates for polymerase
us
chain reaction (PCR) amplification. 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 95 °C before
an
PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified
M
(AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Library preparations were sequence on an Illumina Hiseq 2000 platform
d
and 100 bp pair-end reads were generated. The clean Illumina sequencing data have been submitted to the GenBank (SRA) with the accession number SRR1303115.
te
2.4 Mapping of digital gene expression (DGE) reads
Ac ce p
The reference unigenes were assembled from transcriptome (SRA accession SRP018044), and further annotated using Non-redundant protein (Nr) database, Nucleotide collection (Nt) database, Swiss-Plot, Gene Ontology database (GO), Clusters of Orthologous Groups (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008; Ye et al., 2006). After removing reads with adaptors, reads containing poly N and low quality reads from raw data, we obtained clean reads. Bowtie v0.12.9 was used to align the high-quality reads to the reference unigene sequences. Numbers of the reads mapped to each unigene were determined using RESM (Li et al., 2011). RPKM (reads per kilobase of exon model per million mapped reads) of each gene was calculated, considering the effects of sequencing depth and gene length for the reads count at the same time. 2.4 Quality control Raw data (raw reads) of fastq format were firstly processed through in-house perl 7
Page 7 of 24
scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated. All the downstream analyses were based on the clean data with high quality.
ip t
2.5 Analysis of differential expression analysis
The differential expression detection of genes across samples was performed
cr
using the DEGSeq R package (1.12.0). The resulting p-values were fitted using the
Benjamini and Hochberg’s approach. Corrected p value of 0.05 and log2 (fold change)
us
of 1 were selected as the threshold for significantly differential expression. 2.6 GO and KEGG
an
Gene Ontology (GO) enrichment analysis of differential expressed genes was conducted by the GOseq package, in which gene length was taken into consideration.
M
A corrected p value less than 0.05 was considered as a threshold for significant enrichment of the gene sets. The topGO software was used to visualize enriched GO
d
categories.
Pathway enrichment analysis can further identify significantly enriched
te
metabolic pathways or signal transduction pathways based on the KEGG database
Ac ce p
resource. KOBAS (2.0) software was utilized to test the statistical enrichment of differential expression genes in KEGG pathways. Pathways with a corrected p value ≤ 0.05 were designated as significantly enriched pathways in DGEs. 2.7 Quantitative real time PCR validation Ten genes were selected for the confirmation of DGE data by Quantitative
real-time PCR (qRT-PCR) and18s was used as a reference gene. 1 µg of total RNA from each sample was used to synthesize the first strand cDNA using the PrimeScriptTM RT Reagent Kit (TaKaRa). The qRT-PCR amplifications were
performed with a Thermo PIKOREAL 96 Real-Time PCR system using SYBR Green Supermix (TaKaRa) according to the manufacturer’s instructions. The qRT-PCR program was 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s, and 52 °C for 20 s for annealing, and then 72 °C for 30 s for extension. The relative expression ratio of a target gene is calculated as described in the 2−ΔΔCt method (Livak and Schmittgen, 8
Page 8 of 24
2001), and expressed in comparison to a reference gene. Statistical analysis was carried out with the statistical software SPSS 19.0. All quantitative data were expressed as means ± S.D. 3. Result
ip t
3.1. Analysis of DGE libraries
Two DGE libraries were composed of digestive glands of the control and
cr
PAHs-exposed using Hiseq 2000 sequencing system (Illumina). Total 12,711,432 and 11,630,310 raw reads were generated in the control and PAHs-exposed DGE libraries,
us
respectively. A total of 12,571,017 and 11,500,871 clean reads were left after removing reads with adaptors, reads containing poly N and low quality reads from
an
raw data. Values of the Q20 percentage and GC percentage of the clean reads in two DGE libraries were 97.40% and 97.46%, 43.87% and 42.92%, respectively (Table 1).
M
So the two DGE libraries were reliable. 3.2. Analysis of mapping
d
We mapped the sequences of the two DGE libraries to the reference transcriptome database of C.farreri, which contains 92999 unigenes. In control and
te
PAHs-treated DGE libraries, 77.07% and 78.60% of the clean reads were mapped to a
Ac ce p
gene in the reference database (Table 1). The number of clean reads mapped to a gene was normalized to RPKM values to
accurately measure the gene expression level (Fig.1). The results showed that mRNAs transcribed from the major types of genes were represented in low RPKM values, and only a small proportion of genes was highly expressed. 3.3. GO analysis and KEGG pathways A total of 469 differentially expressed genes were detected, of which 209
were up-regulated and 260 were down-regulated. To understand the functions of these differentially expressed genes, 151 differentially expressed genes were mapped to terms in GO database and compared with the transcriptome background. The top 30 most highly enriched GO terms were presented in Fig. 2. GO has three ontologies, namely, molecular function, cellular component, and biological process, and genes related to the terms “catalytic activity”, “oxidoreductase activity”, “extracellular 9
Page 9 of 24
space” were dominant in the PAHs-exposed scallop. Interestingly, some genes related to the establishment of meiotic spindle localization, norepinephrine transmembrane transporter activity were all down-regulated in PAHs-exposed digestive gland of C. farreri.
ip t
KEGG assignments were used to classify functional annotations of the identified
genes to further understand the biological functions of the genes. Among the
cr
differentially expressed genes, 613 were assigned to 212 pathways in the KEGG database in 4008 DGE class. 19 terms were significantly enriched (Q value ≤ 0.05). 10
us
significantly enriched KEGG pathways (Corrected p < value) were listed in Table 2. Especially, “Glutathione metabolism” pathway related to detoxification metabolism
an
was significantly enriched. In addition, most of the genes related to pathways “Lysosome”, “Fat digestion and absorption” and “Glutathione metabolism” were
M
up-regulated, while the majority of genes related to pathways “ECM-receptor interaction”, “Betalain biosynthesis” were down-regulated.
d
3.4. Identification and verification of differentially expressed genes The expression profile of 7 transcripts of C. farreri exposed to 0.02, 0.2 and 2
te
μg/L of PAHs was performed by qRT-PCR at 0, 1, 6 and 10 days. The relative mRNA
Ac ce p
expression level of 6 genes including 4 up-regulated, 2 down-regulated genes from DGE libraries was detected by qRT-PCR. These chosen genes list in Table 3. were based on the relation with stress response, detoxication, antioxidation processes that could be important for the exposure event. The relative expression levels of the target genes were calculated with the 2−ΔΔCt
method (Livak and Schmittgen, 2001) (Fig.3). The data revealed that CYP450 3A2, CYP 450 2P1, quinone oxidoreductase (QO) and GST were extremely significantly up-regulated after 1 day exposure at all three concentrations, while HSP70 and Cu/Zn-SOD down-regulated. The HSP70 mRNA expression remained low and significantly different from control values until the end of the experiment (p < 0.01).
But the differences among exposure days were non-significant (p > 0.05). The levels of CYP2P1, CYP3A2 and QO were induced by PAHs at first 1 day then raised gradually by exposure period prolonged. The different between control and treated 10
Page 10 of 24
groups was extremely significant (p < 0.01). GST gene expression in 0.02 and 0.2 μg/L PAHs group peaked at day 6 and remained constant until the end of the experiment, while expression in 2 μg/L PAHs group increased during the whole experiment. The differences of exposed groups were non-significant (p > 0.05). The
the differences among concentrations were remarkable (p < 0.01).
ip t
Cu/Zn-SOD expression PAHs exposure group reached the lowest level at day 10 and
cr
The values of expression were imported to Microsoft Excel and data analyses
were performed using IBM SPSS statistical 19.0 (SPSS Inc.). All data were presented
us
as means ± S.D. (n=3) of the relative mRNA expression. 4. Discussion
an
At present, most researches of differentially expressed genes were focused on single PAH exposure, but less on a mixture of PAHs (Oh et al., 2012, van Delft et al.,
M
2012). High-throughput RNA-sequencing (RNA-seq) technologies, including Roche/454 pyrosequencing, Illumina/Solexa sequencing technology have led to a
d
revolution in the fields of transcriptomics and genome characterization in recent years (Martin and Wang, 2011). Previously, Oh et al. (2012) observed that CYP450 2P1
te
mRNA in liver of medaka was significantly induced by synergistic effects of cPAHs
Ac ce p
mixture, but not by single BaP. On the contrary, that study found that genes encoding malate dehydrogenase, anti-thrombin III, and transferrin were significantly up-regulated in response to exposure to BaP, but not in response to the cPAHs mixture, suggesting that exposure to cPAHs results in a reduction of the effect of BaP. Considering the complex chemical interactions between BaP and other three PAHs, we compared the DGE library with another one conducted by our colleagues (Cai et al., 2014). Under the same environment condition, single BaP exposure (1μg/L)
induced 477 differentially expressed genes in digestive gland of Chlamys farreri at day 10. GST mRNA was down-regulated under single BaP exposure with 1.0169 fold change, while GST was up-regulated under PAHs mixture exposure with 1.1036 fold change. CYP450 3A2 was up-regulated in both single and mix PAH exposure, but with different fold changes, 1.2754 and 1.6721. Besides that, mRNA expression of HSP70, CYP450 2P1 and quinone oxidoreductase were not significantly altered due 11
Page 11 of 24
to single BaP, but to PAHs mixture. In comparison with 469 differentially expressed genes identified in our attempt, the result indicated that impact of PAHs mixture on gene expression could not be identified as synergy or antagonism respectively. Considering vast contaminations in ocean ecosystem, researches about effect of PAHs
ip t
mixture in gene expression are meaningful, but needed further efforts. Furthermore, as presented in Fig.3, qRT-PCR results demonstrated the reliability and accuracy of DGE
cr
data.
Molecular biomarkers can be used to investigate the impact or the physiological
us
response of organism to aquatic pollution. Differentially expressed genes in the digestive gland of C.farreri exposed to PAHs mixture under experimental conditions
closely related with detoxification (Table 3).
an
were investigated. Some genes were identified from the DGE library, of which they
M
Heat shock proteins (HSPs) have been identified in several bivalve species in response to various chemical, physical and pathogenic stresses and appear to represent
d
a general marker of non-specific stress (Qiu et al., 2007). According to the apparent molecular mass, HSPs are classified into several families, including HSP90 (83–90
te
kDa), HSP70 (66–78 kDa), HSP60, HSP47, and small heat shock protein (sHSP;
Ac ce p
15–43 kDa) (Morimoto, 1993). All organisms from bacteria to mammals exposed to environmental stressors respond by synthetizing heat-shock protein (HSP70) (Schlesinger et al., 1982). HSP70 gene expression is implicated in protein repair, transport and protection from oxidative stress (Contardo-Jara et al., 2010). In zebra mussels, an exposure to BaP induced a persistent increase in HSP70 mRNA level after 28 days (Châtel et al., 2012). Similarly, in a field study, Copat et al. (2012) found that in Engraulis encrasicolus tissues collected from the Catania Gulf containing 16 PAHs, HSP70 protein was weakly expressed. Interestingly, our result showed that gene HSP70 expression was remarkably inhibited after one day exposure, but approximately invariant at day 6 and 10 in the digestive gland of scallop, even the highest mRNA relative expression was 0.6 fold compared to control samples. It was induced that gene HSP70 stabilization is likely a strategy for PAHs mixture tolerance in C.farreri, which accordance with previous research in Tar Pond killifish (Christine 12
Page 12 of 24
et al., 2009). However, in the present study, HSP70 mRNA expression was not correlated with the amount of PAHs present in the different stations measured previously in the laboratory (Michel et al., 2013), which was parallel to our data. This may suggest that other parameters such as temperature and pH might also affect
ip t
HSP70 mRNA level induction, as demonstrated in other field studies (Hamer et al.,
2004). The fact is that HSP proteins participate in several important processes such as
cr
folding of the nascent polypeptide chain and protection of proteins from denaturation or aggregation, the suppression of HSPs by PAHs mixture may result in a reduction in
us
protein stability and consequent accumulation of abnormal proteins, leading to changes in cell metabolism and apoptosis.
an
Previous studies have indicated the digestive gland was the main target organ for detoxification of environmental pollutant and identification of related biomarkers.
M
Among differentially expressed genes, we identified several genes related to detoxification. CYP450 is a major family of enzymes involved in the primary or
d
phase I metabolism of xenobiotics, including drugs, carcinogens and pesticides (Anzenbacher and Anzenbacherova, 2001). We clearly observed that PAHs mixture
te
provoked differential expression of cytochrome P450 genes, including CYP2P1,
Ac ce p
CYP3A2. Oh et al. (2012) discovered that CYP450 2P1 mRNA was significantly up-regulated in response to exposure to cPAHs, but not BaP alone. In addition, the gene CYP3A2 was only linked to the benzo[a]pyrene metabolism (Lupp et al., 2001). Our qRT-PCR results revealed that CYP2P1 mRNA expression (4.6 fold change) was more severely changeable than CYP3A2 (1.3 fold change) when scallops coped with 10 day PAHs mixture stress. It may suggest that CYP450 2P1 mRNA can be induced by synergistic effects of cPAHs mixture, while the synergistic effects less worked on CYP450 3A2. Moreover, CYP450 2P1 could be regarded a better candidate biomarker than CYP450 3A2 to monitor complex marine environment. It is proved that there were many kinds of BaP phase I metabolites in mammals, such as BaP epoxides, phenols, dihydrodiols, tetrahydrotetrol, quinones and diolepoxides. Furthermore, BaP appears to be metabolised predominantly to quinones in mussel rather than diols as in vertebrates (Sjölin and Livingstone, 1997), and 13
Page 13 of 24
relative to mechanisms of toxicity, quinones undergo redox cycling with production of O2•− making oxidative stress a prominent feature of PAH toxicity in mussel. These processes cause interference with electron transport, either inducing enzyme systems such as cytochrome P450 that mediate oxidation reactions or depleting protective
ip t
antioxidants. Fortunately, enzymatic (e.g., SOD, GST and GPx) and non-enzymatic
(e.g., GSH, ascorbate, carotenoids, and free proline) antioxidant defense systems in
cr
living organisms can provide cellular protection against oxidative stress. Otherwise, as the pollutants are metabolized, reactive molecules are produced and interact with
us
the genetic material and cause DNA damage as a toxic terminal of the induction of a cascade of cellular events. DNA integrity reduction represents the genotoxicity which
an
may have a long-term effect on the sustainability of a particular population. In the present study, quinone oxidoreductase (QO) mRNA expression was up-regulated
M
remarkable, suggesting that scallops were able to mount a successful adaptive response towards the threats posed by PAHs. Because QO is a principal phase II
d
enzyme that plays an important role in the complete reduction and detoxification of highly reactive quinines. In addition to metabolizing foreign compounds, QO also has
te
antioxidant functions, such as protecting against lipid peroxidationand defending
Ac ce p
against intracellular oxidative stress by scavenging superoxide (Siegel et al., 2004, Talalay and Dinkova-Kostova, 2004). Thus, our results indicated that PAHs mixture exposure may disturb the redox balance of the tested tissues. To further examine the functional distribution of these differentially expressed
genes and classify the functional annotations of them, we applied KEGG analysis on these identified genes. Higher corrected p-value take part in glutathione metabolism
(ko: 00480) were found in our database. Glutathione (GSH) is reported to act as an effective antioxidant in marine animal system and as a reactant in conjugation with electrophilic substances (Kovacevic et al., 2008). Glutathione S-transferase (GST) is a dominated antioxidant enzyme, conjugating electrophilic compounds with GSH. Previous field studies have shown a positive correlation between PAHs burden and GST activity (Gowland et al., 2002). Although GSTs mRNA expression levels could be induced in BaP treated groups (Liu et al., 2014), the GST gene was not 14
Page 14 of 24
differentially regulated as expected in the DGE library. PAHs exposure experiment was carried out to ascertain how GST worked during the 10 days. In the present research, variation of GST mRNA expression was slight in all treated columns, but the difference with the control group was significant (p < 0.05). The result inferred GST
ip t
expression was induced by individual aryl hydrocarbon receptor agonists (BaP), but weakened by PAHs mixture containing non-agonists (Garner et al., 2012). What’s
cr
more, the differences among concentration gradients were positive in our result. Xu et al. (2010) also observed a good linear correlation between the BaP exposed
us
concentrations and the GSTs expression levels at 10 d in the digestive gland of scallops (p < 0.01). So GSTs have been recognized as biomarkers of exposure to
an
oxidative stress-inducing chemical pollutants including PAHs mixture (Boutet et al., 2004).
M
Another key gene belong to the antioxidant defense system is superoxide dismutase (SOD), which can eliminate reactive oxygen species (ROS) produced by
d
phase I reactions. ROS can be responsible for lipid peroxidation processes and eventually result in DNA damage. Wang et al. (2011) reported that Cu/Zn-SOD
te
mRNA expression levels in the digestive gland of R. philippinarum were induced by
Ac ce p
0.01 mg/L BaP. However, in the present study, Cu/Zn-SOD mRNA expression was down-regulated obvious even in the 0.02μg/L PAHs group (p < 0.05), indicating suppression of Cu/Zn-SOD mRNA expression caused by PAHs mixture. But Liu et al. (2014) hold the opinion that Mn-SOD appears to be a better biomarker for BaP exposure than Cu/Zn-SOD. While, in the present study, all three treatments inhibited the SOD mRNA expression significantly at the whole exposure time, and the expression decreased along with mixture concentration increased, indicating significant dose-time effect. Thus, Cu/Zn-SOD could reflect the PAHs pollution at low concentration. 5. Conclusion Our study examined the transcriptome response of C.farreri to PAHs mixture stress using digital gene expression analysis. The results demonstrated that PAHs mixture significantly altered expression of genes implicated in stress response, 15
Page 15 of 24
detoxication, antioxidation and so on. As to assess multi-PAHs pollution based on scallops C.farreri, CYP450 2P1 could offer more alert information than CYP3A2. Besides, regarding QO as PAHs pollution biomarker applied in scallop was rarely reported. Moreover, PAHs mixture inhibits the expression of a key gene involved in
ip t
glutathione metabolism and may disturb the detoxication system of C.farreri, which would put forward new evidence about making GST as biomarker in bivalve. The
cr
current results provide a strong basis for future research on molecular mechanism response to PAHs exposure to bivalve. More validation experiments are also needed
us
to confirm if these differentially expressed genes identified in this study can serve as
an
potential biomarkers for PAHs pollution.
Acknowledgements
M
This work was supported by the National Natural Sciences Foundation of China (30972237) and the State Ocean Administration Special Public Project of China
d
(201105013). We thank the staff at the Laboratory of Environmental Physiology of Aquatic Animal for the help with sampling and taking care of the scallops.
te
Conflict of interest statement
Ac ce p
The authors have nothing to disclose. References
Anzenbacher, P., Anzenbacherova, E., 2001. Cytochromes P450 and metabolism of xenobiotics. CMLS, Cell. Mol. Life Sci. 58, 737-747.
Boutet, I., Tanguy, A., Moraga, D., 2004. Characterisation and expression of four mRNA sequences encoding glutathione S-transferases pi, mu, omega and sigma classes in the Pacific oyster Crassostrea gigas exposed to hydrocarbons and pesticides. Mar. Biol., 146, 53–64 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, 128-135.
Cai, Y.F., Pan, L.Q., Hu, F.X., Jin, Q., Liu, T., 2014. Deep sequencing-based transcriptome profiling analysis of Chlamys farreri exposed to benzo[a]pyrene. 16
Page 16 of 24
GENE. 551, 261-270. Châtel, A., Faucet-Marquis, V., Perret, M., Gourlay-Francé, C., Uher, E., Pfohl-Leszkowicz, A., Vincent-Hubert, F., 2012. Genotoxicity assessment and detoxification induction in Dreissena polymorpha exposed to benzo[a]pyrene
ip t
MUTAGENESIS. 27, 703–711
Christine Paetzold, S., Ross, N. W., Richards, R. C., Jones, M., Hellou, J., & Bard, S.
cr
M., 2009. Up-regulation of hepatic ABCC2, ABCG2, CYP1A1 and GST in
multixenobiotic-resistant killifish (Fundulus heteroclitus) from the Sydney Tar
us
Ponds, Nova Scotia, Canada. Mar. Environ. Res. 68, 37-47.
Chung, M., Hu, R., Cheung, K., Wong, M., 2007. Pollutants in Hong Kong soils:
an
polycyclic aromatic hydrocarbons. CHEMOSPHERE. 67, 464–473. Connon, R., Hooper, H.L., Sibly, R.M., Lim, F.L., Heckmann, L.H., Moore, D.J.,
M
Watanabe, H., Soetaert, A., Cook, K., Maund, S.J., Hutchinson, T.H., Moggs, J., Coen, W.D., Iguchi, T., Callaghan, A., 2008. Linking molecular and population
42, 2181-2188.
d
stress responses in Daphnia magna exposed to cadmium. Environ. Sci. Technol.
te
Contardo-Jara, V., Pflugmacher, S., Nutzmann, G., Kloas, W., Wiegand C., 2010. The
Ac ce p
beta receptor blocker metoprolol alters detoxification processes in the non-target organism Dreissena polymorpha. Environ. Pollut., 158, 2059–2066
Copat, C., Brundo, M.V., Arena, G., Grasso, A., Oliveri Conti, G., Ledda, C., Fallico, R., Sciacca, S., Ferrante, M., 2012. Seasonal variation of bioaccumulation in Engraulis encrasicolus (Linneaus, 1758) and related biomarkers of exposure.
Ecotox. Envieon. Safe. 86, 31-37.
EFSA, 2008. Scientific Opinion of the Panel on Contaminants in the Food Chain on a request from the European Commission on Polycyclic Aromatic Hydrocarbons in Food. EFSA Journal, 724, 1–114 Fu, X., Sun, Y., Wang, J., Xing, Q., Zou, J., Li, R.,Wang, Z., Wang, S., Hu, X., Zhang, L., Bao, Z., 2014. Sequencing-based gene network analysis provides a core set of gene resource for understanding thermal adaptation in Zhikong scallop Chlamys farreri. Mol. Ecol. Resour. 14, 184-198. 17
Page 17 of 24
Garner, L. V., & Di Giulio, R. T., 2012. Glutathione transferase pi class 2 (GSTp2) protects against the cardiac deformities caused by exposure to PAHs but not PCB-126 in zebrafish embryos. Comp. Biochem. Phys. C. 155, 573-579. Goldberg, E.D., Bowen, V.T., Farrington, J.W., Harvey, G., Martin, J.H., Parker, P.L.,
ip t
Risebrough, R.W., Robertson, W., Schneider, E., Gamble, E., 1978. The mussel watch. Enviro.Cnserv. 5, 101-125.
cr
Gowland, B.T.G., Mcintosh, A.D., Davies, I.M., Moffat, C.F., Webster, L., 2002.
Implications from a field study regarding the relationship between polycyclic
us
aromatic hydrocarbons and glutathione S-transferase activity in mussels. Mar. Environ. Res. 54, 231-235.
an
Hamer, B., Hamer, D.P., Muller, W.E., Batel, R., 2004. Stress-70 proteins in marine mussel Mytilus galloprovincialis as biomarkers of environmental pollution: a
M
field study. Environ. Int. 30, 873–882.
Huang, W., Wang, Z., Yan, W., 2012. Distribution and sources of polycyclic aromatic
d
hydrocarbons (PAHs) in sediments from Zhanjiang Bay and Leizhou Bay, South China. Mar. Pollut. Bull. 64, 1962–1969.
te
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., Katayama, T.,
Ac ce p
Kawashima, S., Okuda, S., Tokimatsu, T.,Yamanishi, Y., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, 480-484.
Kovacevic, T.B., Borokovic, S.S., Pavlovic, S.Z., Despotovic, S.G., Saicic, Z.S., 2008. Glutathione as a suitable biomarker in hepatopancreas, gills and muscles of three fresh water cryfish species. Arch. Biol. Sci. 60, 59–66.
Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC BIOINFORMATICS. 12, 323.
Liu, D., Pan, L., Cai, Y., Li, Z., & Miao, J., 2014. Response of detoxification gene mRNA expression and selection of molecular biomarkers in the clam Ruditapes philippinarum exposed to benzo[a]pyrene. Environ. Pollut. 189, 1-8. Livak, K.J., Schmittgen, T.D. 2001. Analysis of relative gene expression data using real-time
quantitative
PCR
and
the
2(T)(−delta
delta
C)
method.
METHODS .25,402-408. 18
Page 18 of 24
Lupp, A., Danz, M., Müller, D., 2001. Morphology and cytochrome P450 isoforms expression in precision-cut rat liver slices. TOXICOLOGY. 161, 53-66. Martin, J. A., Wang, Z., 2011. Next-generation transcriptome assembly. Nat. Rev. Genet. 12, 671-682..
ip t
Men, B., He, M., Tan, L., Lin, C., Quan, X., 2009. Distributions of polycyclic
aromatic hydrocarbons in the Daliao River Estuary of Liaodong Bay, Bohai Sea
cr
(China). Mar. Pollut. Bull. 58, 818-826.
Michel, C., Bourgeault, A., Gourlay-Francé, C., Palais, F., Geffard, A.,
us
Vincent-Hubert, F., 2013. Seasonal and PAH impact on DNA strand-break levels in gills of transplanted zebra mussels. Ecotoxicol. Environ. Saf. 92, 18–26.
an
Morimoto, R.I., 1993. Cell in stress: transcriptional activation of heat shock genes. SCIENCE. 259, 1409–1410, at doi: http://dx.doi.org/10.1126/science 8451637.
M
Oh, J.H., Moon, H.B., Choe, E.S., 2012. Alterations in differentially expressed genes by exposure to a mixture of carcinogenic polycyclic aromatic hydrocarbons in
d
the liver of Oryzias latipes. Environ. Toxicol. Phar. 33, 403-407. Pan, L., Ren, J., Liu, J., 2005. Effects of benzo[k]fluoranthene exposure on the
te
biomarkers of scallop Chlamys farreri. Comp. Biochem. Physiol. C: Toxicol.
Ac ce p
Pharmacol. 141, 248–256.
Qiu, L.M., Song, L.S., Yu, Y.D., Xu, W., Ni, D.J., Zhang, Q.C., 2007. Identification and characterization of a myeloid differentiation factor 88 (MyD88) cDNA from Zhikong scallop Chlamys farreri. Fish Shellfish Immun. 23, 614–623.
Schlesinger, M.J., Ashburner, M., Tissieres, A., 1982. Heat-Shock from Bacteria to Men 1982 Cold Spring Harbor Laboratory, New York, p. 286
Siegel, D., Gustafson, D.L., Dehn, D.L., Han, J.Y., Boonchoong, P., Berliner, L.J., Ross, D., 2004. NAD(P)H: quinone oxidoreductase 1: role as a superoxide scavenger. Mol. Pharmocol. 65, 1238-1247. Sjölin, A. M., Livingstone, D. R., 1997. Redox cycling of aromatic hydrocarbon quinones catalysed by digestive gland microsomes of the common mussel (Mytilus edulis L.). Aquat. Toxicol. 38, 83-99. Talalay, P., Dinkova-Kostova, A.T., 2004. Role of nicotinamide quinone 19
Page 19 of 24
oxidoreductase 1 (NQO1) in protection against toxicity of electrophiles and reactive oxygen intermediates. Method Enzymol. 382, 355-364. Wang, L., Pan, L.Q., Liu, N., Liu, D., Xu, C.Q., Miao, J.J., 2011. Biomarkers and bio-accumulation of clam Ruditapes philippinarum in response to combined
ip t
cadmium and benzo[a]pyrene exposure. Food Chem. Toxicol. 49, 3407-3417.
Xu, C., Pan, L., Liu, N., Wang, L., & Miao, J., 2010. Cloning, characterization and
cr
tissue distribution of a pi-class glutathione S-transferase from clam (Venerupis
philippinarum): Response to benzo[a]pyrene exposure. Comp. Biochem. Phys. C.
us
152, 160-166.
Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., Wang, J., Li, S., Li, R.,
an
Bolund, L.,Wang, J., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34, 293- 297.
M
Zhao, X.L., Yu, H., Kong, L.F., Li, Q., 2012. Transcriptomic responses to salinity
Ac ce p
te
d
stress in the pacific oyster Crassostrea gigas. PloS one. 7, e46244.
20
Page 20 of 24
Table 1. Control
PAHs-exposed
12711432
11630310
Clean reads
12571017
11500871
0.03
0.03
Q20 (%)
97.46
97.40
Total reads
12571017
11500871
Total mapped (%)
77.07
78.60
an
us
cr
Error rate (%)
ip t
Category Raw reads
Summary for the C.farreri DGE in control and PAHs-exposed libraries
M
Table 2. The 10 most significantly enriched KEGG pathway (Corrected p value < 0.05) Pathway
Up-/down-regulated sequences
Corrected p value
ko04974
Protein digestion and absorption
7/9
2.99E-05
ko04512
ECM-receptor interaction
5/12
0.0015
ko04640
Hematopoietic cell lineage
5/3
0.0025
ko04142
Lysosome
15/2
0.0040
te
d
Pathway ID
Tyrosine metabolism
3/6
0.0040
ko04975
Fat digestion and absorption
7/0
0.0040
ko04614
Renin-angiotensin system
4/1
0.0040
ko00965
Betalain biosynthesis
1/3
0.0040
ko04145
Phagosome
9/8
0.0050
ko00480
Glutathione metabolism
5/2
0.0236
Ac ce p
ko00350
Table 3. Genes selected for Real-time qRT-PCR. Gene number
Description
log2.Fold_changea
CL6151.Contig1_zhikongshanbei
Hsp70,putative [Ixodes scapularis]
-2.4552
CL8118.Contig1_zhikongshanbei
Cytochrome P450 2P1[Fundulus heteroclitus]
1.2745
Unigene64457_zhikongshanbei
Cytochrome P450 3A2 [Rattus norvegicus]
1.6127
CL4845.Contig2_zhikongshanbei
sulfide:quinone oxidoreductase [Arenicola marina]
2.6718
CL10249.Contig1_zhikongshanbei
Glutathione S-transferase
1.1036
CL1826.Contig1_zhikongshanbei
Superoxide dismutase [Cu-Zn] [Brugia pahangi]
-3.2591
a
Up- or down-regulated in DGE of PAHs exposed clams C.farreri (10 d).
21
Page 21 of 24
us
cr
ip t
Fig. 1
Ac ce p
te
d
M
an
Fig.2
22
Page 22 of 24
te
d
M
an
us
cr
ip t
Fig.3
Fig.1 Length distribution of clean reads mapped.
Ac ce p
Fig.2 Histogram presentation of Gene Ontology classification. The results indicated the most enriched GO terms (PAHs exposed vs Control). Terms were summarized in three main categories: biological process, cellular component and molecular function. * indicated the enriched GO term. Fig.3 Relative expression by Quantitative real-time PCR of 6 selected genes. The Ct values of 0d were used as baseline with the value 1. The relative expression levels of the genes were calculated using 2
−ΔΔCt
model. ** Significantly different from the control group (p < 0.01).
23
Page 23 of 24
Ac ce p
te
d
M
an
us
cr
ip t
Conflict of interest statement The authors have nothing to disclose.
24
Page 24 of 24