- Email: [email protected]
RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta
Accepted Manuscript RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta Zhenhui Wang, Yina Shao, Chenghua Li, Weiwei Zhang, Xuemei Duan, Xuelin Zhao, Qiongfen Qiu, Chunhua Jin PII:
S1050-4648(16)30535-6
DOI:
10.1016/j.fsi.2016.08.051
Reference:
YFSIM 4152
To appear in:
Fish and Shellfish Immunology
Received Date: 15 June 2016 Revised Date:
14 August 2016
Accepted Date: 24 August 2016
Please cite this article as: Wang Z, Shao Y, Li C, Zhang W, Duan X, Zhao X, Qiu Q, Jin C, RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta, Fish and Shellfish Immunology (2016), doi: 10.1016/j.fsi.2016.08.051. 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.
ACCEPTED MANUSCRIPT
RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam
RI PT
Sinonovacula constricta
Zhenhui Wang, Yina Shao, Chenghua Li*, Weiwei Zhang, Xuemei Duan, Xuelin
SC
Zhao, Qiongfen Qiu, Chunhua Jin
M AN U
School of Marine Sciences, Ningbo University, Ningbo 315211, PR China
* Corresponding author:
818 Fenghua Road,
EP
Ningbo University,
TE D
Chenghua Li
AC C
Ningbo, Zhejiang Province 315211, P. R. China
Tel: 86-574-87608368
Fax: 86-574-87608368,
Email: [email protected] Notes The authors declare no competing financial interest.
ACCEPTED MANUSCRIPT
RNA-seq analysis revealed ROS-mediated related genes
Sinonovacula constricta
RI PT
involved in cadmium detoxification in the razor clam
SC
Abstract: The razor clam Sinonovacula constricta is a eukaryotic benthic intertidal bivalve species that is tolerant to different heavy metals, such as cadmium ion (Cd2+).
M AN U
However, the mechanism by which S. constricta responds to Cd2+-induced stress remains unclear. In this study, eight transcriptome libraries were constructed and sequenced using razor clams exposed to Cd2+ for 12 and 48 h. A total of 18,330 unigenes with an average length of 500 bp were annotated. Among these 18,330
TE D
unigenes, 582 and 649 displayed differential expression profiles at 12 and 48 h in the gill, respectively. The corresponding differential unigenes in the hepatopancreas were 1056 and 382. Gene Ontology annotation revealed that these unigenes were highly
EP
pronounced in metabolic process, cellular process, binding, and catalytic activity. Notably, ROS production-related genes, such as heat shock proteins 32,
AC C
metallothionein, and glutathione, were synchronously enriched in all experimental samples with induced expression profiles, which was also validated by qPCR. Our results highlighted the relation between immune regulation and Cd2+-induced stress in razor clam and provided new insights into the molecular mechanisms of heavy toxicology. Keywords: Sinonovacula constricta, Cadmium, RNA-seq, Heat shock protein32, Metallothionein, Glutathione.
ACCEPTED MANUSCRIPT
1. Introduction The razor clam Sinonovacula constricta is a mollusk that thrives in the lower-tomid intertidal zones along the coast of the West Pacific Ocean. Razor clam is a
RI PT
commercially important benthic shellfish and one of the four major clam species in Chinese aquaculture production. In 2015, the cultured razor clam yield was approximately 720,084 tons, which accounted for 30% of the total mudflat shellfish
SC
production in China [1]. Many efforts have been exerted to understand the ecological habits [2], nutritional components [3], cultural techniques [4], and genetic structure
M AN U
and biochemical characteristics [5] of this species. Hence, an increasing number of razor clam genes have been discovered and deposited in GenBank. However, the number of identified genes in razor clam is fewer than those in other bivalves, such as oyster, scallop, and clam [6,7]. Moreover, most sequences are containing
TE D
microsatellite sequences or annotated as mitochondrial genes [8,9]. Few functional genes associated with growth and immune regulation have been well documented; these genes include insulin-like growth factor, growth and differentiation factor 8,
EP
selenium-binding protein [10], and small heat shock proteins [11]. Razor clam genes specific to heavy mental exposure were scarcely investigated.
AC C
With the rapid development of coastal industries, heavy metal pollution of the
ocean is worsening [12]. Cadmium ion (Cd2+) poisons multi-organs and has a long half-life. The element has been listed as the 6th most poisonous substance by the Agency for Toxic Substances and Disease Registry (ATSDR). The great influence on food safety of increased heavy metal contents in shellfish has become a concern with the improvement of health standards. Copper or cadmium exposure decreases the feeding activity of Ruditapes philippinarum by prolonging shell closure [13]. Haberkorn et al. [14] verified that cadmium exposure of Crassostrea gigas modifies
ACCEPTED MANUSCRIPT the lipid composition of the digestive gland; this modification indicates alterations in digestion, resulting in reduced toxin accumulation. Cadmium exposure of oysters may also increase phenoloxidase activity and circulating hemocytes, and decrease phagocytosis, to compensate for the increase in their mortality [15]. In mammals, Cd-
RI PT
induced stress increases the expression level of heat shock protein 32 (HSP32) as an alternative factor of metallothionein (MT) through ROS augmentation [16]. However, this process remains to be comprehensively investigated in a specific species. In
SC
addition, preliminary groundwork and foundation preparation work must be conducted to provide a breakthrough point for the specific and in-depth research of
M AN U
complex heavy metal detoxification mechanisms of S. constricta. RNA-seq technologies are efficient and reliable for transcriptome sequencing [17]. Recent advances in RNA-seq have allowed us to generate an unprecedented global view of the transcriptome and provided us with an efficient method of exploring the whole
TE D
transcriptional landscape [17]. For S. constricta, effective molecular markers and certain key genes associated with growth have been identified and validated [19], providing an efficient and reliable method for regulatory gene screening in non-model
EP
organisms. The present study is the first to analyze comparatively the transcriptional profiles of razor clam in Cd2+-exposed gill and hepatopancreas at two observation
AC C
time points, then construct cDNA library and further to mine related genes on regulation and detoxification for cadmium in S. constricta, finally in order to provide basic research direction for molecular pathways mechanism study in the future.
2. Materials and Methods 2.1 Sample collection
ACCEPTED MANUSCRIPT Sixty healthy S. constricta (shell length: 5–6 cm, shell width: 2 2.5 cm, weight: 10 ± 1 g) were obtained from a commercial shellfish farm (Ninghai, Zhejiang, China) in October 2014. All S. constricta individuals used for the study were randomly
RI PT
distributed in six tanks and maintained in a 10 L open-circuit filtered seawater with a salinity level of 20 at 16 °C with aeration [20]. After accommodation for 3 days, three tanks were exposed to 5 µg/L Cd2+ by adding 50 µL of stock solution of CdCl2 (1 g/L).
SC
The other three tanks were used as controls. Identical tissue samples of the hepatopancreas and the gill were collected at 12 and 48 h after the six tanks were
M AN U
exposed to cadmium. Afterward, tissue samples with the same exposure time from three parallel tanks were mixed up to 100 mg and immediately stored with 1 mL of RNAlater (Life technology, Waltham, MA, USA) for later use. Finally, eight groups of samples for constructing cDNA libraries and subsequent real-time quantitative PCR were obtained.
TE D
2.2 RNA extraction and quantification The saved gill and hepatopancreas samples were grinded by adding liquid nitrogen, and each sample from the eight groups was lysed in 1 mL of TRIzol reagent
EP
(Invitrogen, Carlsbad, CA, USA) for RNA extraction in accordance with the
AC C
manufacturer’s instructions. Then, the total RNA concentration was measured on a NanoDrop (Thermo Scientific, Waltham, MA, USA), and RNA integrity [21] was checked using an Agilent 2100 BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). Total RNA was treated with Dnase I (Ambion, Thermo Scientific, Waltham, MA, USA) and Poly-A RNA was extracted from each total RNA sample using MicroPoly(A)-PuristTM Kit (Ambion, Thermo Scientific, Waltham, MA, USA) in accordance with the manufacturers’ protocols. Equal quantities of high-quality mRNA from each material were used for cDNA synthesis.
ACCEPTED MANUSCRIPT
2.3 cDNA library construction and sequencing A single-end fragment library was constructed following the SOLiD Total RNA-Seq Kit protocol (Life Technologies, PN4452437). The mRNA was fragmented by RNase III and purified using the RiboMinus Concentration Module (Invitrogen,
RI PT
Carlsbad, CA, USA). The RNA fragment was linked with the adaptors using the hybridization master mix (SOLiD Total RNA-Seq Kit) before it was subsequently subjected to reverse transcription. The purified cDNA was size-selected after
SC
electrophoresis with the Novex TBE-Urea Gel (Invitrogen, Carlsbad, CA, USA) at 180 V for 20 min as described by Wang et al [22]. The gel block containing 150–250
M AN U
nt cDNA was precisely excised and used as amplification template. The PCR reactions were initiated at 95 °C for 5 min, followed by 16 cycles of 95 °C for 30 s, 62 °C for 3 s, and 72 °C for 30 s in a thermal cycler. All of the components used in the amplification were from the SOLiD Total RNA-Seq Kit. The yield and size
TE D
distribution of the PCR products were verified using Agilent 2100 Bioanalyzer. Emulsion PCR and bead enrichment were performed using the SOLiD EZ BeadTM system (Life Technologies, Waltham, MA, USA). Workflow analysis was
EP
performed to verify the quality and density of the template beads, and the enriched beads for each sample were deposited on the sequencing slide. Finally, the libraries
AC C
were sequenced using the SOLiD 4 sequencing platform, and color-space reads were obtained. The raw sequencing reads were submitted to NCBI Short Read Archive under the accession number of SRP080845.
2.4 Sequence data analysis and functional annotation After the initial images were obtained, the data were transformed into sequence data. The raw reads were cleaned by removing the adaptor sequences and ambiguous or low-quality reads. De novo transcriptome assembly was carried out using the shortread assembly program Trinity consisting of three software modules: Inchworm,
ACCEPTED MANUSCRIPT Chrysalis, and Butterfly, as follows [23]. First, Inchworm assembles reads into the unique sequences of transcripts. Inchworm uses a greedy k-mer (k=25) based approach for fast and efficient transcript assembly, recovering only a single (best) representative for a set of alternative variants that share k-mers (due to alternative
RI PT
splicing, gene duplication, or allelic variation). Then, Chrysalis clusters related contigs that correspond to portions of alternatively-spliced transcripts or otherwise unique portions of paralogous genes. Chrysalis then constructs a de Bruijn graph for
SC
each cluster of related contigs, each reflecting the complexity of overlaps between variants. Finally, Butterfly analyzes the paths taken by reads and read pairings in the
M AN U
context of the corresponding de Bruijn graph and reports all plausible transcript sequences, resolving alternatively spliced isoforms and transcripts derived from paralogous genes [24]. The predicted protein-coding sequences were searched on different databases, including Swiss-Prot(http://www.uniprot.org/), NCBI non-
TE D
redundant (Nr) protein (ftp://ftp.ncbi.nih.gov/blast/db/), Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/), euKaryotic Orthologous Groups (KOG) (http://www.ncbi.nlm.nih.gov/KOG/), Pfam (http://pfam.xfam.org/), and
EP
Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), by using Blastx with an E-value<1e-5. Gene names were assigned to each protein
AC C
sequence on the basis of the best BLAST hit from all BLAST results. GO enrichment analysis was implemented through one-tailed Fisher’s exact test
with a filter value set to 0.001. The differentially expressed genes (DEGs) were selected as the test set, whereas all identified transcripts were used as the reference set. The significantly overenriched GO in the test set was reported, and the DEGs were mapped to the KEGG database using the built-in function of Blast2GO for pathway retrieval. All mapping and retrieving steps were performed using the default setting.
2.5 Gene expression difference analysis
ACCEPTED MANUSCRIPT We counted the number of reads that can be mapped to the assembly contigs from the two tissues obtained at 12 and 48 h, respectively, and transformed them into FPKM (Fragments Per Kilo bases per Million reads) [25]. Then, we calculated the abundance of differential expression of each gene on a contig from the gill and observed
at
different
times
by
using
TBCC
RI PT
hepatopancreas
and an absolute value of log2(fold_change) >1.
2.6 Real-time quantitative PCR
SC
(Tophat+Bowtie+Cufflinks+Cuffdiff) [26] with a false-discovery rate (FDR) <0.001
Animals challenge experiment was performed according to section 2.1. Total
M AN U
RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and the cDNA was synthesized using a PrimeScript RT-PCR Kit (TaKaRa, Shiga, Japan). A real-time quantitative PCR experiment associated with Cd stress such as HSP32, GSH and MT, using SYBR Green detection chemistry (TaKaRa, Shiga, Japan), was
TE D
performed on a Rotor-Gene 6000 real-time quantitative PCR detection system. The primer information for the real-time quantitative PCR is shown in Table 3. Actin of S. constricta (ScActin) served as internal controls to normalize the targets for
EP
quantification. Each reaction was performed in a final volume of 20 µL, which contained 2 µL of the cDNA, 1 µL of each primer (10 µM), 6 µL of RNase-free water
AC C
and 10 µL of the SYBR Green PCR Master Mix (TaKara, Shiga, Japan). The amplification profile was as follows: denaturation at 94 °C for 5 min, followed by 40 cycles of 94 °C for 15 s, 60 °C for 30 s and 70 °C for 30 s. At the end of the PCR cycles, melting curve analyses were performed. The Ct value is defined as the fractional cycle number at which the fluorescence passes the fixed threshold. Each sample was analyzed in quintuplicate.
ACCEPTED MANUSCRIPT
3. Result 3.1 Sequencing and assembly Eight libraries, which include the four control groups (PC library) and four treated groups (PT library), were sequenced using the SOLiD 4 sequencing platform.
RI PT
The characteristics of these libraries are summarized in Table 1. A total of 7,205,607 and 10,301,319 raw reads in the gill were produced from the control groups at 12 and 48 h, respectively, compared with the 9,861,573 and 9,578,956 raw nucleotide reads
SC
from the treated groups. In the hepatopancreas, a total of 7,311,534 and 8,406,059 raw nucleotide reads were generated for the control groups, whereas 22,115,620 and
M AN U
13,666,784 raw nucleotide reads were obtained from the treated groups at the corresponding exposure times. After removing low-quality reads, 45.68% and 42.19% clean reads of the control groups were obtained from the gill, and 42.48% and 45.68% were obtained from the treated groups at 12 and 48 h. The corresponding clean data
TE D
for the hepatopancreas were 49.49%, 41.36%, 52.57%, and 49.18%, respectively. Through assembly and elimination of redundancy, 127,005 contigs were generated from the eight libraries, in which 33,134 contigs were longer than 500 bp. Genes with
EP
different expression patterns were identified between the PC (normal) and PT (stimulation for 12 or 48 h) libraries. All P-values were adjusted with an FDR
AC C
correction for multiple testing by the Benjamini–Hochberg method [27]. Comparison of the reads obtained from the eight libraries revealed that 5,516 unique tags were present between the libraries (absolute value of log2(fold_change) >1) at P < 0.001. The counts of DEGs in the eight libraries are displayed in Fig. 1.
3.2 Functional annotation The contigs were annotated on the basis of the identity of the translation frame, and 18,330 contigs were predicted to play roles in mapping protein-coding sequences. Among the annotated unigenes, 10,984 with a length of more than 300 bp and 7,346
ACCEPTED MANUSCRIPT with a length of more than 1,000 bp were identified. The statistics of annotated genes is displayed in Table 2. The putative protein-coding sequences were searched using BLAST [28] on the Swiss-Prot [29], NCBI Nr protein [30], Gene Ontology (GO) [31], COG 32], KOG [33], and KEGG [34] (E-value<1e-5) databases. Annotated putative
RI PT
proteins were acquired by matching their biological functions to other species, whereas some sequences without annotations possibly represented the poorly conserved region in S. constricta.
SC
As mentioned above, GO analysis was carried out on the putative proteins. Of the 23,152 putative proteins, 7,871 unique sequences on the average were assigned to one
M AN U
or more GO terms, and 21,430 GO assignments were obtained in total. All GO assignments fell under the broad categories for all three major GO functional fields, namely, Biological Processes, Cellular Components, and Molecular Function (Fig. 2). In addition to GO terms, 15,621 unigenes were mapped to 262 KEGG pathways based
TE D
on KO terms for assignments. The well-represented pathways were focal adhesion (234), PI3K-Akt signaling pathway (221), and Rap1 signaling pathway (217). Some immune response-related pathways, such as NF-kappa B signaling pathway (67) and
EP
Toll-like receptor signaling pathway (67), were also identified. Assignments of COG and KOG, which was only used to compare homologous genes in eukaryotes, were
AC C
used to predict and classify the possible functions of these unigenes. On the basis of sequence homology, 4,798 unigenes (26.18%) and 11,244 unigenes (61.34%) were annotated and divided into 25 specific categories (Figs. 3 and 4) for COG and KOG analyses, respectively. In the COG analysis, the general function prediction, which involved 1,617 unigenes (25.21%), was the largest; followed by the categories replication, recombination, and repair (6.56%); post-translational modification, protein turnover, and chaperones (6.47%); and transcription (6.28%). By contrast, the signal transduction mechanisms, which had the maximum value among all the
ACCEPTED MANUSCRIPT functional classification, contained 2,043 unigenes (16.04%), followed by general function prediction (14.97%); post-translational modification, protein turnover, and chaperones (9.63%); and intracellular trafficking, secretion, and vesicular transport (7.52%) for the KOG analysis.
RI PT
3.3 Analysis of differential expression genes between PC and PT libraries Biological functions, which can characterize the responses of S. constricta to heavy metals such as Cd2+, were focused saliently on the sequence libraries.
SC
Comparing their occurrences in the PC and PT libraries, we established a group of figures of transcripts that indicated different expressions of Cd2+ to stress at different
M AN U
time periods in the two types of tissues considered (Fig. 5). Cellular and singleorganism processes both contained the most number of downregulated genes after 12 h (71, 68) and 48 h (55,55). However, the processes with the third largest number of downregulated genes in the gill were metabolic process at 12 h (53) and localization
TE D
process at 48 h (33). Meanwhile, the biological processes with the most number of upregulated genes in the gill at 12 and 48 h were cellular process (72,17), singleorganism process (69,14), and metabolic process (61,14). In the hepatopancreas, the
EP
processes containing the most DEGs at 12 h were cellular, metabolic, and singleorganism processes, and these processes were consistent between downregulated
AC C
(35,32,27) and upregulated genes (98,88,81). Meanwhile, the results were identical with the observations at 12 h but different from the biological processes (i.e., cellular, single-organism, and metabolic processes) observed at 48 h. As for the molecular function, the largest number of processes for the upregulated or downregulated genes appeared during the binding process for each observation time of the sample tissues. For the cellular component, more than 90% of the DEGs were involved in the constitution of the cytoplasmic matrix components; however, only a small portion was involved on the cellular membrane or organelle membrane. Moreover, cellular
ACCEPTED MANUSCRIPT processes were the most prominent among all biological processes in both the hepatopancreas and the gill.
3.4 ROS-mediated genes in response to Cd stress Some ROS-mediated genes in response to Cd stress were discovered through the
RI PT
analysis of DEGs (Table 4). The key DEGs were significantly involved in cellular processes. In particular, HSP32, MT, and glutathione (GSH) were differentially expressed in the two tissues. The expression of the HSP32 gene (highest value of
SC
log2(fold_change) =10.14) was higher in the gill than in the hepatopancreas (the highest value of log2(fold_change) =8.3). Meanwhile, the expression levels of MT
M AN U
(highest value of log2(fold_change) = 5.6) and GSH (highest value of log2(fold_change) =3.4) were higher in the hepatopancreas than in the gill (highest value of log2(fold_change) =1.3 and 0.8, respectively).
3.5 Validation of differential expressed genes by qRT-PCR
TE D
Three key genes associated with Cd detoxification containing HSP32, GSH and MT, which are three of the most significant genes in response to Cd stress of S. constricta, were selected to perform qRT-PCR and compare with the RNA-seq data
EP
by gene’s transcript level. The results showed that the expression patterns of HSP32, GSH and MT all showed induced expression profiles trend by qRT-PCR (Fig. 6),
AC C
which was totally consistent with RNA-seq analysis.
4. Discussion
ROS accumulated in the organism because of the replacement of zinc with Cd2+
in Zn/Cu superoxide dismutase after Cd2+-induced stress. Logically, the key genes related to ROS mediation were significantly expressed during the decrease and elimination of ROS to alleviate cell damage. In fact, recent research has confirmed that the synthesis level of the key gene HSP32 increases via ROS-activated AP-1 [14].
ACCEPTED MANUSCRIPT Meanwhile, the expression levels of genes synergistic with HSP32, such as MT and GSH, significantly increase along with HSP32. In the present study, S. constricta transcriptomes were extracted from the gill and the hepatopancreas at different time periods (12 and 48 h) under Cd2+-induced
RI PT
stress and normal conditions. The upregulated and downregulated genes significantly differed between the control and treated groups at 12 and 48 h. The identified DEGs included HSP32, MT, and GSH. The upregulated level of these genes after Cd2+-
SC
induced stress indicated their important functions in the regulation and detoxification of cadmium, thereby providing a basis for explaining the putative molecular
M AN U
mechanism of regulation and detoxification of cadmium in S. constricta. Cellular processes were the most prominent among all the biological processes in both the gill and the hepatopancreas. The molecular mechanisms behind the response to Cd2+-induced stress of the hepatopancreas is possibly similar to that of the
TE D
gill and other tissues. Thus, the gill and the hepatopancreas are important in the modulation and clearance of cadmium. Tetsuya et al. [14] confirmed that the expression of the HSP32 gene is more augmented in the hepatopancreas and kidney of
EP
mice exposed to cadmium as compared with other organs; this finding explains why the hepatopancreas is not a target organ of cadmium carcinogenesis. In the present
AC C
study, the key DEGs (HSP32, MT, and GSH) between the two tissues were involved in cellular processes and indicated that the hepatopancreas may be more effective than the gill in regulating and eliminating cadmium in S. constricta. Functional annotation of DEGs between the PT and PC libraries revealed biological processes that may characterize the responses to Cd2+-induced stress. Several key genes, such as HSP32 (highest value of log2(fold_change) = 10.14 in the gill at 48 h), MT (highest value of log2(fold_change) = 5.6 in the hepatopancreas at 48 h), and GSH (highest value of log2(fold_change) = 3.4 in the hepatopancreas at 48 h),
ACCEPTED MANUSCRIPT which are related to cellular processes, including metabolic processes and molecular transducer activities, are implicated in the regulation of Cd2+-induced stress response. Moreover, cellular processes related to immune system functions, biological adhesion, and anti-apoplectic processes were all highly represented.
RI PT
Similar to most marine mollusks [35,36], S. constricta can resist and endure the increasing Cd2+ concentration stress in oceans [37]. Heat shock proteins (HSPs), such as HSP32, HSP90, HSP70, and small heat shock proteins (sHSPs), may be the direct
SC
regulators and cause of change in S. constricta transcriptome data on the dominant expression of the genes in the hepatopancreas. Results showed that the expression of
M AN U
HSP32 was the highest among all HSP families with different gene expression patterns. Recent research has shown that physical and chemical stresses, such as heavy metal, acids, and oxidizers [38–40], can activate HSP expression in addition to thermal stimulation [41]. Anguo Zhang et al. [10] found that the expression levels of
TE D
sHSPs significantly increase in hemocytes after the PbCl2 exposure of S. constricta, indicating that HSPs are involved in mediating exposure to environmental pollutants. Tetsuya Abe et al. [14,42] found that the survival rate of control mice after the
EP
subcutaneous injection of CdCl2 is higher than that of HSP32 gene-deficient transgenic mice. This finding indicates that HSPs are obviously important molecules
AC C
for resisting external heavy metal stimulation caused by the transfer of an organism from a normal to a heavy metal laden environment. Recent studies have also indicated that the synthesis of HSPs is promoted by ROS-activated AP-1; the ROS are induced by Cd, which replaced the zinc in Zn/Cu superoxide dismutase that binds with the HSP32 promoter [43,44]. In the present study, Zn/Cu superoxide dismutase (highest absolute value of log2(fold_change) = 2.5) and HSPs, especially HSP32, were significantly upregulated. This result suggests that HSP32 protects against ROS induced by Cd2+ via Zn/Cu superoxide dismutase in S. constricta. Previous research
ACCEPTED MANUSCRIPT confirmed that Na+-K+-ATP enzyme regulates the permeability and transport of metal ions across the cellular membrane [45]. Simultaneously, HSPs can also regulate Na+K+-ATP enzyme activity [40], which was upregulated in the present study. These findings suggest that the HSP32 gene promotes intracellular resistance against ROS
RI PT
caused by Cd2+ exposure and suppresses the cellular transference of Cd2+ via the downregulation of the Na+-K+-ATP enzyme gene, resulting in the stable resistance function and low accumulation of Cd2+ in S. constricta.
SC
Aside from genes participating in Cd2+ resistance channels, genes related to detoxification responses to Cd2+ were also found. To our knowledge, the timely
M AN U
removal of ROS that damage proteins and nucleic acids in cells maintains the normal functioning of the organism [46]. MT functions in the detoxification of Cd2+ by decreasing ROS production [47]. MT is a class of metallothione with a low molecular weight, high cysteine content, and wide abundance in an organism; this type of
TE D
metallothione can be stimulated by heavy metals. The upstream coding sequence of MT exists in various transcriptional regulatory elements, such as metal response elements (MREs), GREs, IL-6, and MLTF/ARE [48]. MREs, which play a positive
EP
regulatory role for MT, were studied more than other transcriptional regulatory elements. Tetsuya Abe et al. [14] demonstrated that the transcription of MT genes
AC C
induced by metals such as zinc and cadmium is triggered through the activation of metal transcription factor-I (MTF-I) that binds to MREs. They also found that the expression of the MT gene is extensively induced by Cd2+ in mice [14]. Snyder. R.D found that a high expression of MT in the hepatopancreas may suppress the induction of ROS by Cd2+ in rats [49]. The present results indicate that MT, as well as MREs (highest value of log2(fold_change) = 3.2) and MTF-I (highest value of log2(fold_change) = 2.8), was up-regulated, especially in the hepatopancreas of the treated groups of S. constricta. This result is consistent with previous functional
ACCEPTED MANUSCRIPT research that revealed a feasible detoxification pathway involved in MT for Cd2+ [50– 53]. Recent studies have also indicated that MT has a higher induction efficiency in the hepatopancreas of C. gigas [54] and kidney of rats after Cd exposure, demonstrating the specificity of MT gene expression in different organs induced by
RI PT
Cd2+ [55]. In the present study, MT was highly expressed in the hepatopancreas (value of log2(fold_change) = 5.6) and gill (value of log2(fold_change) = 3.8) at 48 h. We inferred that the production of ROS increases with increasing Cd2+ concentration. The
SC
expression of MT must be regulated via MTF-I and MREs to decrease ROS and thus protect organisms against damage from oxidative stress.
M AN U
GSH plays an important role in protection against the toxicity of Cd2+ to target organs [56]. For detoxification, GSH integrates with a GSH active group that contains a -SH in cysteine residues and bind to Cd2+ via GSH transferases [57]. Therefore, GSH, especially that in rat hepatopancreas cells, participates in bio-transformation by
TE D
transforming harmful poisons in organs into harmless substances and excreting them outside the body [58]. GSH also suppresses microsome-mediated interactions by acting as a “scavenger” of reactive intermediates [59]. In the present study, RNA-seq
EP
data revealed that the expression of GSH significantly differed in the two tissues (value of log2(fold_change) = 3.4). Considering both recent and present results, we
AC C
inferred that GSH induced by Cd2+ promoted the excretion of Cd2+. In summary, HSP32 modulated and resisted the damage in organs of ROS initially induced by Cd. Afterward, the proteins of MT and GSH induced by ROS through some pathways directly eliminated ROS and promoted the detoxification of Cd2+ in S. constricta. Immune responses, such as apoptosis, programmed cell death, or cell suicide, were also observed. Among the DEGs, apoptosis suppressors ITG-α and ITG-β directly inhibit the activity of caspase [60]. In addition, ITG-β is involved in signal transduction and cell adhesion regulation [61]. These downregulated genes, including
ACCEPTED MANUSCRIPT caspase, trigger apoptosis or are activated in response to other triggers of apoptosis. We also found APAF1-interacting protein, a gene that encodes for a cytoplasmic protein that forms one of the central hubs in the apoptosis regulatory network containing (from the N terminal) a caspase recruitment domain, an ATPase domain
RI PT
(NB-ARC), few short helical domains, and several copies of the WD40 repeat domain. Upon binding cytochrome C and dATP, this protein forms an oligomeric apoptosome [62]. The apoptosome binds and cleaves caspase 9 protein, releasing its mature and
SC
activated form [63]. In the present study, the expression of genes related to apoptosis focused on both the inhibition and promotion of apoptosis. Basing from the present
M AN U
data, we conclude that the balanced promotion and inhibition of apoptotic factors ensure the stability of S. constricta in adapting to an environment laden with heavy
TE D
metals, such as Cd and Pb.
Acknowledgments
This work was financially supported through the Public welfare project in Zhejiang
EP
Province (2015C32004), Natural Science Foundation of Ningbo (2015C10009, 2015C10062), Collaborative Innovation Center for Zhejiang Marine High-Efficiency
AC C
and Healthy Aquaculture, and the K. C. Wong Magna Fund at Ningbo University.
References
[1] Yuan XC, Zhao WW. China Fishery Statistical Yearbook. Fisheries Agency of China Agriculture Ministry. 2015; 56. [2] Wang SK, Jin BS, Qin HM, Sheng Q, Wu JH. Trophic Dynamics of Filter Feeding Bivalves in the Yangtze Estuarine Intertidal Marsh: Stable Isotope and Fatty Acid Analyses. Plos One. 2015; 10(8): e0135604.
ACCEPTED MANUSCRIPT [3] Jin S, Xue CB, Wang GL, Lu TX, Wang YN, Chen YE. Bacterial flora composition and its dynamics in tidal-flat Sinonovacula constricta aquaculture area. Chinese journal of applied ecology. 2005; 16(7): 1380-2. [4] Lu ZB, Du Q, Fang MJ, Qian XM, Cai QH, Xu CY. Carrying capacity of shellfish
16(5):961-6.
RI PT
culture in Dadeng Island sea area of Xiamen. 2005; Ying Yong Sheng Tai Xue Bao.
[5] Niu DH, Wang F, Xie SM, Sun FY, Wang Z, Peng MX et al. Developmental Analysis
and
Identification
of
Genes
Involved
SC
Transcriptome
in
Larval
Metamorphosis of the Razor Clam, Sinonovacula constricta. Marine Biotechnology.
M AN U
2016; 18(2): 168-75.
[6] Feng BB, Dong LL, Niu DH, Meng SS, Zhang B, Liu DB, et al. Identification of Immune Genes of the Agamaki Clam (Sinonovacula constricta) by Sequencing and Bioinformatic Analysis of ESTs. Marine Biotechnology. 2010; 12(3): 182-91.
TE D
[7] Niu DH, Xie SM, Bai ZY, Wang L, Jin K, Li JL. Identification, expression, and responses to bacterial challenge of the cathepsin C gene from the razor clam
241-5.
EP
Sinonovacula constricta. Developmental and Comparative Immunology. 2014; 46(2):
[8] Ma HT, Jiang HB, Liu XQ, Wu XP, Wei XM. Polymorphic microsatellite loci for
AC C
the razor clam, Sinonovacula constricta. Genetics and Molecular Research. 2015; 14(1): 145-8.
[9] Yuan Y, Li Q, Yu H, Kong LF. The Complete Mitochondrial Genomes of Six Heterodont Bivalves (Tellinoidea and Solenoidea): Variable Gene Arrangements and Phylogenetic Implications. PLoS One. 2012; 7(2): e32353. [10] Lu YL, Zhang AG, Li CH, Zhang P, Su XR, Li Y, et al. The link between selenium binding protein from Sinonovacula constricta and environmental pollutions exposure. Fish and Shellfish Immunology. 2013; 35(2): 271-7.
ACCEPTED MANUSCRIPT [11] Zhang AG, Lu YL, Li CH, Zhang P, Su XR, Li Y et al. A small heat shock protein (sHSP) from Sinonovacula constricta against heavy metals stresses. Fish and Shellfish Immunology. 2013; 34(6): 1605-10. [12] Natalie D, Jill MY, Chiu, Beverly HK, Pob, Anna CK et al. Heavy metal
RI PT
contamination along the China coastline: A comprehensive study using Artificial Mussels and native mussels. Journal of Environmental Management. 2016; 180(9): 238-46.
SC
[13] Tran D, Fournier E, Durrieu G, Massabuau JC. Copper detection in the Asiatic
clam Corbicila fluminea: optimum valve closure response. Toxicology. 2004; 66(3):
M AN U
333-43.
[14] Haberkorn H, Lambert C, Le Goic N, Moal J, Palacios E, Lassus P et al. Effects of Alexandrium minutum exposure upon physiological and hematological variables of diploid and triploid oysters, Crassostrea gigas. Aquatic Toxicology. 2010; 2(97): 96108.
[15] Lin W, Rice MA, Chien PK. The differential effects of three heavy metals on
TE D
particle filtration and amino acid uptake in the pacific oyster, Crassostrea gigas. Journal of Shellfish Research. 12: 111.
[16] AbeT, Yamamoto O,Gotoh S, Yan Y, Todaka Y, Higashi K. Cadmium-Induced
EP
mRNA Expression of Hsp32 Is Augmented in Metallothionein-I and -II Knock-out
AC C
Mice. Archives of Biochemistry and Biophysics. 2000; 1(382): 81-8. [17] Zhang PJ, Li CH, Zhu L, Su XR, Li Y, Jin CH, et al. De Novo Assembly of the Sea Cucumber Apostichopus japonicus Hemocytes Transcriptome to Identify miRNA Targets Associated with Skin Ulceration Syndrome. PLoS One. 2013; 8(9): e73506. [18] Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary toolfor transcriptomics. Nature Review Genetics. 2009; 10(1): 57-63. [19] Niu DH, Wang L, Sun FY, Liu ZJ, Li JL. Development of Molecular Resources for an Intertidal Clam, Sinonovacula constricta, Using 454 Transcriptome Sequencing.
ACCEPTED MANUSCRIPT PLoS One. 2013; 8(7): e67456. [20] Moreira R, Pereiro P, Canchaya C, Posada D, Figueras A, Novoa B. RNA-Seq in Mytilus galloprovincialis: comparative transcriptomics and expression profiles among different tissues. BMC Genomics. 2015; 16(1): 728.
RI PT
[21] Schroeder A, Mueller O, Stocker S, Salowsky R, Leiber M, Gassmann M, et al. The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Molecular Biology. 2006; 7:3.
SC
[22] Wang, L. et al. The promotion of cytoskeleton integration and redox in the haemocyte of shrimp Litopenaeus vannamei after the successive stimulation of
M AN U
recombinant VP28. Developmental and Comparative Immunology. 2014; 45: 123-32. [23] Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Fulllength transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology. 2011; 29(7): 644-52.
TE D
[24] Rice P, Longden I, Bleasby A. EMBOSS: the European Molecular Biology Open Software Suite. Trends in Genetics. 2000; 16(6): 276-7. [25] Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and
621-8.
EP
quantifying mammalian transcriptomes by RNA-Seq. Nature Methods. 2008; 5(7):
AC C
[26] Wang LK, Feng ZX, Wang X, Wang XW, Zhang XG. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010; 26(1):136-8.
[27] Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2005; 57: 289-300. [28] Altschul SF, Madden TL, Schaffer AA, et al. Gapped BLAST and PSI BLAST: A New Generation of Protein Database Search Programs. Nucleic Acids Research.
ACCEPTED MANUSCRIPT 1997; 25(17): 3389-402. [29] Apweiler R, Bairoch A, Wu CH, et al. UniProt: the Universal Protein knowledgebase. Nucleic Acids Research. 2004; 32(Database issue): D115-9. [30] Deng YY, Li JQ, Wu SF, et al. Integrated nr Database in Protein Annotation
RI PT
System and Its Localization. Computer Engineering. 2006; 32(5): 71-4.
[31] Ashburner M, Ball CA, Blake JA, et al. Gene ontology: a tool for the unification of biology. Nature genetics. 2000; 25(1): 25-9.
SC
[32] Tatusov RL, Galperin MY, Natale DA. The COG database: a tool for genome scale annalysis of protein functions and evolution. Nucleic Acids Research. 2000;
M AN U
28(1): 33-6.
[33] Koonin EV, Fedorova ND, Jackson JD, et al. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biology. 2004; 5(2): R7.
TE D
[34] Kanehisa M, Goto S, Kawashima S, et al. The KEGG resource for deciphering the genome. Nucleic Acids Research. 2004; 32(Database issue): D277-D280. [35] Berger VJ, Kharazova AD. Mechanisms of salinity adaptation in marine
EP
molluscs. Hydrobiologia. 1997; 355: 115-26. [36] Lockwood BL, Somero GN. Transcriptomic responses to salinity stress in
AC C
invasive and native blue mussels (genus Mytilus). Molecular Ecology. 2011; 20(3): 517-29.
[37] Bergman HL, Dorward-King EJ. Reassessment of metals criteria for aquatic life protection: priorities for research and implementation. Pensacola: SETAC Press. 1997; 31-56. [38] Trautinger F, Kindås-Mügge I, Barlan B, Neuner P, Knobler RM. 72-kD heat shock protein is a mediator of resistance to ultraviolet B light. Journal of Investigative
ACCEPTED MANUSCRIPT Dermatology. 1995; 105(2): 160-2. [39] Lee S, Prochaska DJ, Fang F, Barnum SR. A 16.6-kilodalton protein in the Cyanobacterium synechocystis sp. PCC 6803 plays a role in the heat shock response. Current Microbiology. 1998; 37(6): 403-7.
RI PT
[40] Köhler HR, Triebskorn R, Stöcker W, Kloetzel PM, Alberti G. The 70 kD heat shock protein (HSP 70) in soil invertebrates: a possible tool for monitoring environmental toxicants. Archives of environmental contamination and toxicology.
SC
1992; 22(4): 334-8.
[41] Joshi CP, Nguyen HT. 5' untranslated leader sequences of eukaryotic mRNAs
M AN U
encoding heat shock induced proteins. Nucleic Acids Research. 1995; 23(4): 541-9. [42] Chen Q, Osteryoung K, Vierling E. A 21-kDa chloroplast heat shock protein assembles into high molecular weight complexes in vivo and in Organelle. The Journal of Biological Chemistry. 1994; 269(18): 1216-23.
TE D
[43] Kusumoto T, Maehara Y, Sakaguchi Y, Emi Y, Kohnoe S, Sugimachi K. 1Hexylcarbamoyl-5-fluorouracil alters the expression of heat shock protein in HeLa cells. Anticancer Drugs. 1991; 2(1): 45-8.
EP
[44] Riordan M, Sreedharan R, Wang S, Thulin G, Mann A, Stankewich M, et al. HSP70 binding modulates detachment of Na-K-ATPase following energy deprivation
AC C
in renal epithelial cells. American Journal of Physiology Renal Physiology. 2005; 288(6): F1236-42.
[45] Kroese LJ, Scheffer PG. 8-hydroxy-2'-deoxyguanosine and cardiovascular disease: a systematic review. Current Atherosclerosis Reports. 2014; 16(11): 452. [46] Adams TK, Saydam N, Steiner F, Schaffner W, Freedman JH. Activation of gene expression by metal-responsive signal transduction pathways. Environ Health Perspect. 2002; 110 Suppl 5: 813-7. [47] Dorian C, Gattone VH, Klaassen CD. Accumulation and degradation of the
ACCEPTED MANUSCRIPT protein moiety of cadmium-metallothionein (CdMT) in the mouse kidney. Toxicology and applied pharmacology. 1992; 117(2): 242-8. [48] Dalton T, Fu K, Enders GC, Palmiter RD, Andrews GK. Analysis of the effects of overexpression of metallothionein-I in transgenic mice on the reproductive
RI PT
toxicology of cadmium. Environmental Health Perspectives. 1996; 104(1): 68-76.
[49] Snyder RD. A review and investigation into the mechanistic basis of the genotoxicity of antihistamines. Mutation Research. 1998; 411(3): 235-48.
SC
[50] Abdel-Mageed AB, Agrawal KC. Activation of nuclear factor kappaB: potential role in metallothionein-mediated mitogenic response. Cancer Research. 1998; 58(11):
M AN U
2335-8.
[51] Lahiri A, Abraham C. Activation of pattern recognition receptors up-regulates metallothioneins, thereby increasing intracellular accumulation of zinc, autophagy, and bacterial clearance by macrophages. Gastroenterology. 147(4): 835-46.
TE D
[52] Galdes A, Hill HA, Kägi JH, Vasak M, Bremner I, Young BW. The 1H N.M.R. spectra of metallothioneins. Experientia Suppl. 1979; 34:241-8. [53] Tsujikawa K, Suzuki N, Shimaoka T, Iida K, Kohama Y, Otaki N, et
EP
al.Abnormal accumulation of copper-metallothionein in the hepatopancreas and kidney of Long-Evans rats with a cinnamon-like coat color (LEC rats).Biological &
AC C
pharmaceutical bulletin. 1994; 17(5): 591-5. [54] Zhao XL, Yu H, Kong LF, Li Q. Transcriptomic responses to salinity stress in the Pacific oyster Crassostrea gigas. PLoS One. 2012; 7(9): e46244. [55] Liu J, LiuY, MichalskaAE, ChooKH, KlaasenCD. Metallothionein plays less of a protective role in cadmium-metallothionein-induced nephrotoxicity than in cadmium chloride-induced hepatotoxicity. The Journal of Pharmacology. 1996; 276(3): 121623.
ACCEPTED MANUSCRIPT [56] Gazdag Z, Fujs S, Koszegi B, Kálmán N, Papp G, Emri T, et al. The abc1-/coq8respiratory-deficient mutant of Schizosaccharomyces pombe suffers from glutathione underproduction and hyperaccumulates Cd2+. Folia Microbiologica. 2011; 56(4): 353-
RI PT
9. [57] Dourado DF, Fernandes PA, Mannervik B, Ramos MJ. Glutathione transferase: new model for glutathione activation. Chemistry. 2008; 14(31): 9591-8.
SC
[58] Tan KH, Meyer DJ, Gillies N, Ketterer B. Detoxification of DNA hydroperoxide by glutathione transferases and the purification and characterization of glutathione
841-5.
M AN U
transferases of the rat hepatopancreas nucleus. Biochemical Journal. 1988; 254(3):
[59] Mazzullo M, Colacci A, Grilli S, Prodi G, Arfellini G. 1,1,2-Trichloroethane: evidence of genotoxicity from short-term tests. Japanese journal of cancer research : Gann. 1986; 77(6): 532-9.
TE D
[60] Murphy BJ, Andrews GK, Bittel D, Discher DJ, McCue J, Green CJ, et al. Activation of metallothionein gene expression by hypoxia involves metal response elements and metal transcription factor-1. Cancer Research. 1999; 59(6): 1315-22.
EP
[61] Eckelman BP, Salvesen GS, Scott FL. Human inhibitor of apoptosis proteins:
AC C
why XIAP is the black sheep of the family.EMBO Reports. 2006; 7(10): 988-94. [62] Richter BWM, Duckett CS. The IAP proteins: caspases inhibitors and beyond. Science STKE. 2000; 44: 1-4. [63] Ko DC, Gamazon ER, Shukla KP, Pfuetzner RA, Whittington D, Holden TD, et al. Functional genetic screen of human diversity reveals that a methionine salvage enzyme regulates inflammatory cell death. Proceedings of the National Academy of Sciences of the United States of America. 2012; 109(35): E2343-52.
ACCEPTED MANUSCRIPT
Figure legends
RI PT
Fig 1. Summary statistics of the counts of the different expression genes in gill (A) and hepatopancreas (B).
Fig 2. Classification of the annotated amino-acid sequences. Amino-acid sequences were grouped into different functional sub-categories: cellular component, molecular
SC
function and biological process.
Fig 3. COG function classification consensus sequence.
M AN U
Fig 4. KOG function classification consensus sequence.
Fig 5. Summary statistics of the counts of genes involving in different processes according to GO analysis. A: counts of genes down-regulated in gill 12 h vs control; B: counts of genes up-regulated in gill 12 h vs control; C: counts of genes downregulated in gill 48 h vs control; D: counts of genes up-regulated in gill 48 h vs control; E: counts of genes down-regulated in hepatopancreas 12 h vs control; F:
TE D
counts of genes up-regulated in hepatopancreas 12 h vs control; G: counts of genes down-regulated in hepatopancreas 48 h vs control; H: counts of genes up-regulated in hepatopancreas 48 h vs control.
Fig 6. Validation of relative expression level of GSH, MT and HSP32 by qRT-PCR.
EP
A: Relative expression level of GSH in gill; B: Relative expression level of MT in gill; C: Relative expression level of HSP32 in gill; D: Relative expression level of GSH in
AC C
hepatopancreas; E: Relative expression levelof MT in hepatopancreas; F: Relative expression level of HSP32 in hepatopancreas.
ACCEPTED MANUSCRIPT
Table 1. Summary statistics of the transcriptome sequencing for S. constricta from
Gill controlled
(12h)
(48h)
Gill treated (48h)
Hepatopancreas
Hepatopan-
Hepatopancreas
controlled
creas treated
Hepatopancreas controlled
(12h)
(12h)
(48h)
(48h)
SC
(12h)
Gill treated
M AN U
Gill controlled
RI PT
control group (PC) and treated group (PT)
treated
7,205,607
9,861,573 10,301,319
9,578,956
7,311,534
22,115,620
8,406,059
13,666,784
clean reads
3,291,760
4,189,228
4,346,584
4,375,850
3,618,758
11,665,491
3,476,846
6,720,966
clean reads (%)
45.68%
42.48%
42.19%
45.68%
49.49%
52.57%
41.36%
49.18%
mapped reads
4,616,632
6,685,160
mapped reads (%)
64.07%
67.79%
TE D
Raw sequencing reads
6,091,258
4,921,393
15,107,180
5,786,731
9,613,215
65.29%
63.95%
67.31%
68.31%
68.84%
70.34%
AC C
EP
6,725,731
ACCEPTED MANUSCRIPT Table 2. Summary statistics of the number of unigenes mapping all annotated databases
Unigene
≥300nt
≥1000nt
COG_Annotation
4,798
2,436
2,362
GO_Annotation
6,173
3,539
KEGG_Annotation
9,314
5,325
KOG_Annotation
11,244
6,408
Pfam_Annotation
12,971
7,107
Swissprot_Annotation
8,964
nr_Annotation
17,770
All_Annotated
18,330
SC
M AN U
TE D EP AC C
RI PT
Annotated databases
2,634 3,989 4,836
5,864
4,848
4,116
10,585
7,185
10,984
7,346
ACCEPTED MANUSCRIPT Table 3. Real-time quantitative PCR primer information in this study
Names
Sequences(5′-3′)
Application
ScHSP32
GGCATCTGGATTGGAACTATGATAG
Real-time PCR
CTTACCCCAACTGCATAAATGTTTC GTGTATGAGTCCGTGGTGGTGTGTA GAATGCAGGTCCTCTCGTATCTCAT ScMT
ATGAGGCATGTATAGACAGACAGTG
ScActin
CGGGATCCATGGCTGAGACAATG
EP
TE D
M AN U
CCCAAGCTTCTAGCTAAGGAGTTTCTGGTC
AC C
Real-time PCR
SC
GTCTCATTAAGGTCTCCATGGTAAC
Real-time PCR
RI PT
ScGSH
Real-time PCR
ACCEPTED MANUSCRIPT Table 4. Detailed annotations of differentially expressed genes associated with cadmium detoxification
Gene name
Highest absolute value log2(fold_change)
P-value
UN000905
HSP32 △
10.14
0.0001036
UN014045
MT △
5.6
0.0002495
Expression tissue/time
Gene Ontology
Up/Down
gill/48 h
F: catalytic activity (GO:0003824); P: metabolic process (GO:0008152); F: metal ion binding (GO:0046872)
up
P: response to metal ion (GO:0010038)
up
F: catalytic activity (GO:0003824)
up
gill/48 h
F: ATP binding (GO:0005524); P: response to stress (GO:0006950)
up
gill/24 h
P: response to stress (GO:0006950); F: unfolded protein binding (GO:0051082)
down
P: intracellular protein transport (GO:0006886); F: protein transporter activity (GO:0008565);P:vesic le-mediated transport (GO:0016192)
up
F: metal ion binding (GO:0046872)
down
F:sodium-potassiumexchani-ng ATPase activity (GO:0005391); C: integral component of membrane (GO:0016021); F: metal ion binding (GO:0046872)
down
SC
Hepatopancreas/48 h
UN106383
HSP70
2.76
HSP90
-1.78
0.000252 6
0.0003581
0.0004574
AP-1
2.34
AC C
UN068669
UN017282
Zn/Cu superoxide dismutase
Na+-K+ATPase
Hepatopa-
ncreas/48 h
Hepatopa-
0.0003513 ncreas/24 h
EP
UN066133
3.4
M AN U
UN081663
GSH △
TE D
UN005612
Hepatopa-2.50
0.0005741 ncreas/24 h
-1.91
RI PT
Transcript ID
0.0004097
gill/48 h
ACCEPTED MANUSCRIPT
MREs
3.2
0.0001533 ncreas/48 h
HepatopaUN017157
MTF-1
2.8
0.0003507
M AN U TE D EP AC C
ion
F: metal ion binding (GO:0046872); F: organic cyclic compound binding (GO:0097159); F: heterocyclic compound binding (GO:1901363)
SC
ncreas/48 h
F: manganese binding (GO:0030145)
RI PT
HepatopaUN009810
up
up
ACCEPTED MANUSCRIPT
Figure 1
SC
RI PT
A
AC C
EP
TE D
M AN U
B
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 2
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 3
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 4
ACCEPTED MANUSCRIPT
Figure 5
SC
RI PT
A
M AN U
B
AC C
C
EP
TE D
E
ACCEPTED MANUSCRIPT
SC
RI PT
D
EP AC C
F
TE D
M AN U
E
ACCEPTED MANUSCRIPT
SC
RI PT
G
AC C
EP
TE D
M AN U
H
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 6
ACCEPTED MANUSCRIPT
Highlights: Eight cDNA libraries were constructed and sequenced in Cd2+ exposed Sinonovacula constricta.
RI PT
Total of 18,330 unigenes with an average length of 500 bp were annotated. Different expression unigenes were identified in gill and liver at two examined time points.
SC
ROS production-related genes were synchronously enriched in all experimental
AC C
EP
TE D
M AN U
samples
S1050-4648(16)30535-6
DOI:
10.1016/j.fsi.2016.08.051
Reference:
YFSIM 4152
To appear in:
Fish and Shellfish Immunology
Received Date: 15 June 2016 Revised Date:
14 August 2016
Accepted Date: 24 August 2016
Please cite this article as: Wang Z, Shao Y, Li C, Zhang W, Duan X, Zhao X, Qiu Q, Jin C, RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam Sinonovacula constricta, Fish and Shellfish Immunology (2016), doi: 10.1016/j.fsi.2016.08.051. 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.
ACCEPTED MANUSCRIPT
RNA-seq analysis revealed ROS-mediated related genes involved in cadmium detoxification in the razor clam
RI PT
Sinonovacula constricta
Zhenhui Wang, Yina Shao, Chenghua Li*, Weiwei Zhang, Xuemei Duan, Xuelin
SC
Zhao, Qiongfen Qiu, Chunhua Jin
M AN U
School of Marine Sciences, Ningbo University, Ningbo 315211, PR China
* Corresponding author:
818 Fenghua Road,
EP
Ningbo University,
TE D
Chenghua Li
AC C
Ningbo, Zhejiang Province 315211, P. R. China
Tel: 86-574-87608368
Fax: 86-574-87608368,
Email: [email protected] Notes The authors declare no competing financial interest.
ACCEPTED MANUSCRIPT
RNA-seq analysis revealed ROS-mediated related genes
Sinonovacula constricta
RI PT
involved in cadmium detoxification in the razor clam
SC
Abstract: The razor clam Sinonovacula constricta is a eukaryotic benthic intertidal bivalve species that is tolerant to different heavy metals, such as cadmium ion (Cd2+).
M AN U
However, the mechanism by which S. constricta responds to Cd2+-induced stress remains unclear. In this study, eight transcriptome libraries were constructed and sequenced using razor clams exposed to Cd2+ for 12 and 48 h. A total of 18,330 unigenes with an average length of 500 bp were annotated. Among these 18,330
TE D
unigenes, 582 and 649 displayed differential expression profiles at 12 and 48 h in the gill, respectively. The corresponding differential unigenes in the hepatopancreas were 1056 and 382. Gene Ontology annotation revealed that these unigenes were highly
EP
pronounced in metabolic process, cellular process, binding, and catalytic activity. Notably, ROS production-related genes, such as heat shock proteins 32,
AC C
metallothionein, and glutathione, were synchronously enriched in all experimental samples with induced expression profiles, which was also validated by qPCR. Our results highlighted the relation between immune regulation and Cd2+-induced stress in razor clam and provided new insights into the molecular mechanisms of heavy toxicology. Keywords: Sinonovacula constricta, Cadmium, RNA-seq, Heat shock protein32, Metallothionein, Glutathione.
ACCEPTED MANUSCRIPT
1. Introduction The razor clam Sinonovacula constricta is a mollusk that thrives in the lower-tomid intertidal zones along the coast of the West Pacific Ocean. Razor clam is a
RI PT
commercially important benthic shellfish and one of the four major clam species in Chinese aquaculture production. In 2015, the cultured razor clam yield was approximately 720,084 tons, which accounted for 30% of the total mudflat shellfish
SC
production in China [1]. Many efforts have been exerted to understand the ecological habits [2], nutritional components [3], cultural techniques [4], and genetic structure
M AN U
and biochemical characteristics [5] of this species. Hence, an increasing number of razor clam genes have been discovered and deposited in GenBank. However, the number of identified genes in razor clam is fewer than those in other bivalves, such as oyster, scallop, and clam [6,7]. Moreover, most sequences are containing
TE D
microsatellite sequences or annotated as mitochondrial genes [8,9]. Few functional genes associated with growth and immune regulation have been well documented; these genes include insulin-like growth factor, growth and differentiation factor 8,
EP
selenium-binding protein [10], and small heat shock proteins [11]. Razor clam genes specific to heavy mental exposure were scarcely investigated.
AC C
With the rapid development of coastal industries, heavy metal pollution of the
ocean is worsening [12]. Cadmium ion (Cd2+) poisons multi-organs and has a long half-life. The element has been listed as the 6th most poisonous substance by the Agency for Toxic Substances and Disease Registry (ATSDR). The great influence on food safety of increased heavy metal contents in shellfish has become a concern with the improvement of health standards. Copper or cadmium exposure decreases the feeding activity of Ruditapes philippinarum by prolonging shell closure [13]. Haberkorn et al. [14] verified that cadmium exposure of Crassostrea gigas modifies
ACCEPTED MANUSCRIPT the lipid composition of the digestive gland; this modification indicates alterations in digestion, resulting in reduced toxin accumulation. Cadmium exposure of oysters may also increase phenoloxidase activity and circulating hemocytes, and decrease phagocytosis, to compensate for the increase in their mortality [15]. In mammals, Cd-
RI PT
induced stress increases the expression level of heat shock protein 32 (HSP32) as an alternative factor of metallothionein (MT) through ROS augmentation [16]. However, this process remains to be comprehensively investigated in a specific species. In
SC
addition, preliminary groundwork and foundation preparation work must be conducted to provide a breakthrough point for the specific and in-depth research of
M AN U
complex heavy metal detoxification mechanisms of S. constricta. RNA-seq technologies are efficient and reliable for transcriptome sequencing [17]. Recent advances in RNA-seq have allowed us to generate an unprecedented global view of the transcriptome and provided us with an efficient method of exploring the whole
TE D
transcriptional landscape [17]. For S. constricta, effective molecular markers and certain key genes associated with growth have been identified and validated [19], providing an efficient and reliable method for regulatory gene screening in non-model
EP
organisms. The present study is the first to analyze comparatively the transcriptional profiles of razor clam in Cd2+-exposed gill and hepatopancreas at two observation
AC C
time points, then construct cDNA library and further to mine related genes on regulation and detoxification for cadmium in S. constricta, finally in order to provide basic research direction for molecular pathways mechanism study in the future.
2. Materials and Methods 2.1 Sample collection
ACCEPTED MANUSCRIPT Sixty healthy S. constricta (shell length: 5–6 cm, shell width: 2 2.5 cm, weight: 10 ± 1 g) were obtained from a commercial shellfish farm (Ninghai, Zhejiang, China) in October 2014. All S. constricta individuals used for the study were randomly
RI PT
distributed in six tanks and maintained in a 10 L open-circuit filtered seawater with a salinity level of 20 at 16 °C with aeration [20]. After accommodation for 3 days, three tanks were exposed to 5 µg/L Cd2+ by adding 50 µL of stock solution of CdCl2 (1 g/L).
SC
The other three tanks were used as controls. Identical tissue samples of the hepatopancreas and the gill were collected at 12 and 48 h after the six tanks were
M AN U
exposed to cadmium. Afterward, tissue samples with the same exposure time from three parallel tanks were mixed up to 100 mg and immediately stored with 1 mL of RNAlater (Life technology, Waltham, MA, USA) for later use. Finally, eight groups of samples for constructing cDNA libraries and subsequent real-time quantitative PCR were obtained.
TE D
2.2 RNA extraction and quantification The saved gill and hepatopancreas samples were grinded by adding liquid nitrogen, and each sample from the eight groups was lysed in 1 mL of TRIzol reagent
EP
(Invitrogen, Carlsbad, CA, USA) for RNA extraction in accordance with the
AC C
manufacturer’s instructions. Then, the total RNA concentration was measured on a NanoDrop (Thermo Scientific, Waltham, MA, USA), and RNA integrity [21] was checked using an Agilent 2100 BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). Total RNA was treated with Dnase I (Ambion, Thermo Scientific, Waltham, MA, USA) and Poly-A RNA was extracted from each total RNA sample using MicroPoly(A)-PuristTM Kit (Ambion, Thermo Scientific, Waltham, MA, USA) in accordance with the manufacturers’ protocols. Equal quantities of high-quality mRNA from each material were used for cDNA synthesis.
ACCEPTED MANUSCRIPT
2.3 cDNA library construction and sequencing A single-end fragment library was constructed following the SOLiD Total RNA-Seq Kit protocol (Life Technologies, PN4452437). The mRNA was fragmented by RNase III and purified using the RiboMinus Concentration Module (Invitrogen,
RI PT
Carlsbad, CA, USA). The RNA fragment was linked with the adaptors using the hybridization master mix (SOLiD Total RNA-Seq Kit) before it was subsequently subjected to reverse transcription. The purified cDNA was size-selected after
SC
electrophoresis with the Novex TBE-Urea Gel (Invitrogen, Carlsbad, CA, USA) at 180 V for 20 min as described by Wang et al [22]. The gel block containing 150–250
M AN U
nt cDNA was precisely excised and used as amplification template. The PCR reactions were initiated at 95 °C for 5 min, followed by 16 cycles of 95 °C for 30 s, 62 °C for 3 s, and 72 °C for 30 s in a thermal cycler. All of the components used in the amplification were from the SOLiD Total RNA-Seq Kit. The yield and size
TE D
distribution of the PCR products were verified using Agilent 2100 Bioanalyzer. Emulsion PCR and bead enrichment were performed using the SOLiD EZ BeadTM system (Life Technologies, Waltham, MA, USA). Workflow analysis was
EP
performed to verify the quality and density of the template beads, and the enriched beads for each sample were deposited on the sequencing slide. Finally, the libraries
AC C
were sequenced using the SOLiD 4 sequencing platform, and color-space reads were obtained. The raw sequencing reads were submitted to NCBI Short Read Archive under the accession number of SRP080845.
2.4 Sequence data analysis and functional annotation After the initial images were obtained, the data were transformed into sequence data. The raw reads were cleaned by removing the adaptor sequences and ambiguous or low-quality reads. De novo transcriptome assembly was carried out using the shortread assembly program Trinity consisting of three software modules: Inchworm,
ACCEPTED MANUSCRIPT Chrysalis, and Butterfly, as follows [23]. First, Inchworm assembles reads into the unique sequences of transcripts. Inchworm uses a greedy k-mer (k=25) based approach for fast and efficient transcript assembly, recovering only a single (best) representative for a set of alternative variants that share k-mers (due to alternative
RI PT
splicing, gene duplication, or allelic variation). Then, Chrysalis clusters related contigs that correspond to portions of alternatively-spliced transcripts or otherwise unique portions of paralogous genes. Chrysalis then constructs a de Bruijn graph for
SC
each cluster of related contigs, each reflecting the complexity of overlaps between variants. Finally, Butterfly analyzes the paths taken by reads and read pairings in the
M AN U
context of the corresponding de Bruijn graph and reports all plausible transcript sequences, resolving alternatively spliced isoforms and transcripts derived from paralogous genes [24]. The predicted protein-coding sequences were searched on different databases, including Swiss-Prot(http://www.uniprot.org/), NCBI non-
TE D
redundant (Nr) protein (ftp://ftp.ncbi.nih.gov/blast/db/), Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/), euKaryotic Orthologous Groups (KOG) (http://www.ncbi.nlm.nih.gov/KOG/), Pfam (http://pfam.xfam.org/), and
EP
Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), by using Blastx with an E-value<1e-5. Gene names were assigned to each protein
AC C
sequence on the basis of the best BLAST hit from all BLAST results. GO enrichment analysis was implemented through one-tailed Fisher’s exact test
with a filter value set to 0.001. The differentially expressed genes (DEGs) were selected as the test set, whereas all identified transcripts were used as the reference set. The significantly overenriched GO in the test set was reported, and the DEGs were mapped to the KEGG database using the built-in function of Blast2GO for pathway retrieval. All mapping and retrieving steps were performed using the default setting.
2.5 Gene expression difference analysis
ACCEPTED MANUSCRIPT We counted the number of reads that can be mapped to the assembly contigs from the two tissues obtained at 12 and 48 h, respectively, and transformed them into FPKM (Fragments Per Kilo bases per Million reads) [25]. Then, we calculated the abundance of differential expression of each gene on a contig from the gill and observed
at
different
times
by
using
TBCC
RI PT
hepatopancreas
and an absolute value of log2(fold_change) >1.
2.6 Real-time quantitative PCR
SC
(Tophat+Bowtie+Cufflinks+Cuffdiff) [26] with a false-discovery rate (FDR) <0.001
Animals challenge experiment was performed according to section 2.1. Total
M AN U
RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and the cDNA was synthesized using a PrimeScript RT-PCR Kit (TaKaRa, Shiga, Japan). A real-time quantitative PCR experiment associated with Cd stress such as HSP32, GSH and MT, using SYBR Green detection chemistry (TaKaRa, Shiga, Japan), was
TE D
performed on a Rotor-Gene 6000 real-time quantitative PCR detection system. The primer information for the real-time quantitative PCR is shown in Table 3. Actin of S. constricta (ScActin) served as internal controls to normalize the targets for
EP
quantification. Each reaction was performed in a final volume of 20 µL, which contained 2 µL of the cDNA, 1 µL of each primer (10 µM), 6 µL of RNase-free water
AC C
and 10 µL of the SYBR Green PCR Master Mix (TaKara, Shiga, Japan). The amplification profile was as follows: denaturation at 94 °C for 5 min, followed by 40 cycles of 94 °C for 15 s, 60 °C for 30 s and 70 °C for 30 s. At the end of the PCR cycles, melting curve analyses were performed. The Ct value is defined as the fractional cycle number at which the fluorescence passes the fixed threshold. Each sample was analyzed in quintuplicate.
ACCEPTED MANUSCRIPT
3. Result 3.1 Sequencing and assembly Eight libraries, which include the four control groups (PC library) and four treated groups (PT library), were sequenced using the SOLiD 4 sequencing platform.
RI PT
The characteristics of these libraries are summarized in Table 1. A total of 7,205,607 and 10,301,319 raw reads in the gill were produced from the control groups at 12 and 48 h, respectively, compared with the 9,861,573 and 9,578,956 raw nucleotide reads
SC
from the treated groups. In the hepatopancreas, a total of 7,311,534 and 8,406,059 raw nucleotide reads were generated for the control groups, whereas 22,115,620 and
M AN U
13,666,784 raw nucleotide reads were obtained from the treated groups at the corresponding exposure times. After removing low-quality reads, 45.68% and 42.19% clean reads of the control groups were obtained from the gill, and 42.48% and 45.68% were obtained from the treated groups at 12 and 48 h. The corresponding clean data
TE D
for the hepatopancreas were 49.49%, 41.36%, 52.57%, and 49.18%, respectively. Through assembly and elimination of redundancy, 127,005 contigs were generated from the eight libraries, in which 33,134 contigs were longer than 500 bp. Genes with
EP
different expression patterns were identified between the PC (normal) and PT (stimulation for 12 or 48 h) libraries. All P-values were adjusted with an FDR
AC C
correction for multiple testing by the Benjamini–Hochberg method [27]. Comparison of the reads obtained from the eight libraries revealed that 5,516 unique tags were present between the libraries (absolute value of log2(fold_change) >1) at P < 0.001. The counts of DEGs in the eight libraries are displayed in Fig. 1.
3.2 Functional annotation The contigs were annotated on the basis of the identity of the translation frame, and 18,330 contigs were predicted to play roles in mapping protein-coding sequences. Among the annotated unigenes, 10,984 with a length of more than 300 bp and 7,346
ACCEPTED MANUSCRIPT with a length of more than 1,000 bp were identified. The statistics of annotated genes is displayed in Table 2. The putative protein-coding sequences were searched using BLAST [28] on the Swiss-Prot [29], NCBI Nr protein [30], Gene Ontology (GO) [31], COG 32], KOG [33], and KEGG [34] (E-value<1e-5) databases. Annotated putative
RI PT
proteins were acquired by matching their biological functions to other species, whereas some sequences without annotations possibly represented the poorly conserved region in S. constricta.
SC
As mentioned above, GO analysis was carried out on the putative proteins. Of the 23,152 putative proteins, 7,871 unique sequences on the average were assigned to one
M AN U
or more GO terms, and 21,430 GO assignments were obtained in total. All GO assignments fell under the broad categories for all three major GO functional fields, namely, Biological Processes, Cellular Components, and Molecular Function (Fig. 2). In addition to GO terms, 15,621 unigenes were mapped to 262 KEGG pathways based
TE D
on KO terms for assignments. The well-represented pathways were focal adhesion (234), PI3K-Akt signaling pathway (221), and Rap1 signaling pathway (217). Some immune response-related pathways, such as NF-kappa B signaling pathway (67) and
EP
Toll-like receptor signaling pathway (67), were also identified. Assignments of COG and KOG, which was only used to compare homologous genes in eukaryotes, were
AC C
used to predict and classify the possible functions of these unigenes. On the basis of sequence homology, 4,798 unigenes (26.18%) and 11,244 unigenes (61.34%) were annotated and divided into 25 specific categories (Figs. 3 and 4) for COG and KOG analyses, respectively. In the COG analysis, the general function prediction, which involved 1,617 unigenes (25.21%), was the largest; followed by the categories replication, recombination, and repair (6.56%); post-translational modification, protein turnover, and chaperones (6.47%); and transcription (6.28%). By contrast, the signal transduction mechanisms, which had the maximum value among all the
ACCEPTED MANUSCRIPT functional classification, contained 2,043 unigenes (16.04%), followed by general function prediction (14.97%); post-translational modification, protein turnover, and chaperones (9.63%); and intracellular trafficking, secretion, and vesicular transport (7.52%) for the KOG analysis.
RI PT
3.3 Analysis of differential expression genes between PC and PT libraries Biological functions, which can characterize the responses of S. constricta to heavy metals such as Cd2+, were focused saliently on the sequence libraries.
SC
Comparing their occurrences in the PC and PT libraries, we established a group of figures of transcripts that indicated different expressions of Cd2+ to stress at different
M AN U
time periods in the two types of tissues considered (Fig. 5). Cellular and singleorganism processes both contained the most number of downregulated genes after 12 h (71, 68) and 48 h (55,55). However, the processes with the third largest number of downregulated genes in the gill were metabolic process at 12 h (53) and localization
TE D
process at 48 h (33). Meanwhile, the biological processes with the most number of upregulated genes in the gill at 12 and 48 h were cellular process (72,17), singleorganism process (69,14), and metabolic process (61,14). In the hepatopancreas, the
EP
processes containing the most DEGs at 12 h were cellular, metabolic, and singleorganism processes, and these processes were consistent between downregulated
AC C
(35,32,27) and upregulated genes (98,88,81). Meanwhile, the results were identical with the observations at 12 h but different from the biological processes (i.e., cellular, single-organism, and metabolic processes) observed at 48 h. As for the molecular function, the largest number of processes for the upregulated or downregulated genes appeared during the binding process for each observation time of the sample tissues. For the cellular component, more than 90% of the DEGs were involved in the constitution of the cytoplasmic matrix components; however, only a small portion was involved on the cellular membrane or organelle membrane. Moreover, cellular
ACCEPTED MANUSCRIPT processes were the most prominent among all biological processes in both the hepatopancreas and the gill.
3.4 ROS-mediated genes in response to Cd stress Some ROS-mediated genes in response to Cd stress were discovered through the
RI PT
analysis of DEGs (Table 4). The key DEGs were significantly involved in cellular processes. In particular, HSP32, MT, and glutathione (GSH) were differentially expressed in the two tissues. The expression of the HSP32 gene (highest value of
SC
log2(fold_change) =10.14) was higher in the gill than in the hepatopancreas (the highest value of log2(fold_change) =8.3). Meanwhile, the expression levels of MT
M AN U
(highest value of log2(fold_change) = 5.6) and GSH (highest value of log2(fold_change) =3.4) were higher in the hepatopancreas than in the gill (highest value of log2(fold_change) =1.3 and 0.8, respectively).
3.5 Validation of differential expressed genes by qRT-PCR
TE D
Three key genes associated with Cd detoxification containing HSP32, GSH and MT, which are three of the most significant genes in response to Cd stress of S. constricta, were selected to perform qRT-PCR and compare with the RNA-seq data
EP
by gene’s transcript level. The results showed that the expression patterns of HSP32, GSH and MT all showed induced expression profiles trend by qRT-PCR (Fig. 6),
AC C
which was totally consistent with RNA-seq analysis.
4. Discussion
ROS accumulated in the organism because of the replacement of zinc with Cd2+
in Zn/Cu superoxide dismutase after Cd2+-induced stress. Logically, the key genes related to ROS mediation were significantly expressed during the decrease and elimination of ROS to alleviate cell damage. In fact, recent research has confirmed that the synthesis level of the key gene HSP32 increases via ROS-activated AP-1 [14].
ACCEPTED MANUSCRIPT Meanwhile, the expression levels of genes synergistic with HSP32, such as MT and GSH, significantly increase along with HSP32. In the present study, S. constricta transcriptomes were extracted from the gill and the hepatopancreas at different time periods (12 and 48 h) under Cd2+-induced
RI PT
stress and normal conditions. The upregulated and downregulated genes significantly differed between the control and treated groups at 12 and 48 h. The identified DEGs included HSP32, MT, and GSH. The upregulated level of these genes after Cd2+-
SC
induced stress indicated their important functions in the regulation and detoxification of cadmium, thereby providing a basis for explaining the putative molecular
M AN U
mechanism of regulation and detoxification of cadmium in S. constricta. Cellular processes were the most prominent among all the biological processes in both the gill and the hepatopancreas. The molecular mechanisms behind the response to Cd2+-induced stress of the hepatopancreas is possibly similar to that of the
TE D
gill and other tissues. Thus, the gill and the hepatopancreas are important in the modulation and clearance of cadmium. Tetsuya et al. [14] confirmed that the expression of the HSP32 gene is more augmented in the hepatopancreas and kidney of
EP
mice exposed to cadmium as compared with other organs; this finding explains why the hepatopancreas is not a target organ of cadmium carcinogenesis. In the present
AC C
study, the key DEGs (HSP32, MT, and GSH) between the two tissues were involved in cellular processes and indicated that the hepatopancreas may be more effective than the gill in regulating and eliminating cadmium in S. constricta. Functional annotation of DEGs between the PT and PC libraries revealed biological processes that may characterize the responses to Cd2+-induced stress. Several key genes, such as HSP32 (highest value of log2(fold_change) = 10.14 in the gill at 48 h), MT (highest value of log2(fold_change) = 5.6 in the hepatopancreas at 48 h), and GSH (highest value of log2(fold_change) = 3.4 in the hepatopancreas at 48 h),
ACCEPTED MANUSCRIPT which are related to cellular processes, including metabolic processes and molecular transducer activities, are implicated in the regulation of Cd2+-induced stress response. Moreover, cellular processes related to immune system functions, biological adhesion, and anti-apoplectic processes were all highly represented.
RI PT
Similar to most marine mollusks [35,36], S. constricta can resist and endure the increasing Cd2+ concentration stress in oceans [37]. Heat shock proteins (HSPs), such as HSP32, HSP90, HSP70, and small heat shock proteins (sHSPs), may be the direct
SC
regulators and cause of change in S. constricta transcriptome data on the dominant expression of the genes in the hepatopancreas. Results showed that the expression of
M AN U
HSP32 was the highest among all HSP families with different gene expression patterns. Recent research has shown that physical and chemical stresses, such as heavy metal, acids, and oxidizers [38–40], can activate HSP expression in addition to thermal stimulation [41]. Anguo Zhang et al. [10] found that the expression levels of
TE D
sHSPs significantly increase in hemocytes after the PbCl2 exposure of S. constricta, indicating that HSPs are involved in mediating exposure to environmental pollutants. Tetsuya Abe et al. [14,42] found that the survival rate of control mice after the
EP
subcutaneous injection of CdCl2 is higher than that of HSP32 gene-deficient transgenic mice. This finding indicates that HSPs are obviously important molecules
AC C
for resisting external heavy metal stimulation caused by the transfer of an organism from a normal to a heavy metal laden environment. Recent studies have also indicated that the synthesis of HSPs is promoted by ROS-activated AP-1; the ROS are induced by Cd, which replaced the zinc in Zn/Cu superoxide dismutase that binds with the HSP32 promoter [43,44]. In the present study, Zn/Cu superoxide dismutase (highest absolute value of log2(fold_change) = 2.5) and HSPs, especially HSP32, were significantly upregulated. This result suggests that HSP32 protects against ROS induced by Cd2+ via Zn/Cu superoxide dismutase in S. constricta. Previous research
ACCEPTED MANUSCRIPT confirmed that Na+-K+-ATP enzyme regulates the permeability and transport of metal ions across the cellular membrane [45]. Simultaneously, HSPs can also regulate Na+K+-ATP enzyme activity [40], which was upregulated in the present study. These findings suggest that the HSP32 gene promotes intracellular resistance against ROS
RI PT
caused by Cd2+ exposure and suppresses the cellular transference of Cd2+ via the downregulation of the Na+-K+-ATP enzyme gene, resulting in the stable resistance function and low accumulation of Cd2+ in S. constricta.
SC
Aside from genes participating in Cd2+ resistance channels, genes related to detoxification responses to Cd2+ were also found. To our knowledge, the timely
M AN U
removal of ROS that damage proteins and nucleic acids in cells maintains the normal functioning of the organism [46]. MT functions in the detoxification of Cd2+ by decreasing ROS production [47]. MT is a class of metallothione with a low molecular weight, high cysteine content, and wide abundance in an organism; this type of
TE D
metallothione can be stimulated by heavy metals. The upstream coding sequence of MT exists in various transcriptional regulatory elements, such as metal response elements (MREs), GREs, IL-6, and MLTF/ARE [48]. MREs, which play a positive
EP
regulatory role for MT, were studied more than other transcriptional regulatory elements. Tetsuya Abe et al. [14] demonstrated that the transcription of MT genes
AC C
induced by metals such as zinc and cadmium is triggered through the activation of metal transcription factor-I (MTF-I) that binds to MREs. They also found that the expression of the MT gene is extensively induced by Cd2+ in mice [14]. Snyder. R.D found that a high expression of MT in the hepatopancreas may suppress the induction of ROS by Cd2+ in rats [49]. The present results indicate that MT, as well as MREs (highest value of log2(fold_change) = 3.2) and MTF-I (highest value of log2(fold_change) = 2.8), was up-regulated, especially in the hepatopancreas of the treated groups of S. constricta. This result is consistent with previous functional
ACCEPTED MANUSCRIPT research that revealed a feasible detoxification pathway involved in MT for Cd2+ [50– 53]. Recent studies have also indicated that MT has a higher induction efficiency in the hepatopancreas of C. gigas [54] and kidney of rats after Cd exposure, demonstrating the specificity of MT gene expression in different organs induced by
RI PT
Cd2+ [55]. In the present study, MT was highly expressed in the hepatopancreas (value of log2(fold_change) = 5.6) and gill (value of log2(fold_change) = 3.8) at 48 h. We inferred that the production of ROS increases with increasing Cd2+ concentration. The
SC
expression of MT must be regulated via MTF-I and MREs to decrease ROS and thus protect organisms against damage from oxidative stress.
M AN U
GSH plays an important role in protection against the toxicity of Cd2+ to target organs [56]. For detoxification, GSH integrates with a GSH active group that contains a -SH in cysteine residues and bind to Cd2+ via GSH transferases [57]. Therefore, GSH, especially that in rat hepatopancreas cells, participates in bio-transformation by
TE D
transforming harmful poisons in organs into harmless substances and excreting them outside the body [58]. GSH also suppresses microsome-mediated interactions by acting as a “scavenger” of reactive intermediates [59]. In the present study, RNA-seq
EP
data revealed that the expression of GSH significantly differed in the two tissues (value of log2(fold_change) = 3.4). Considering both recent and present results, we
AC C
inferred that GSH induced by Cd2+ promoted the excretion of Cd2+. In summary, HSP32 modulated and resisted the damage in organs of ROS initially induced by Cd. Afterward, the proteins of MT and GSH induced by ROS through some pathways directly eliminated ROS and promoted the detoxification of Cd2+ in S. constricta. Immune responses, such as apoptosis, programmed cell death, or cell suicide, were also observed. Among the DEGs, apoptosis suppressors ITG-α and ITG-β directly inhibit the activity of caspase [60]. In addition, ITG-β is involved in signal transduction and cell adhesion regulation [61]. These downregulated genes, including
ACCEPTED MANUSCRIPT caspase, trigger apoptosis or are activated in response to other triggers of apoptosis. We also found APAF1-interacting protein, a gene that encodes for a cytoplasmic protein that forms one of the central hubs in the apoptosis regulatory network containing (from the N terminal) a caspase recruitment domain, an ATPase domain
RI PT
(NB-ARC), few short helical domains, and several copies of the WD40 repeat domain. Upon binding cytochrome C and dATP, this protein forms an oligomeric apoptosome [62]. The apoptosome binds and cleaves caspase 9 protein, releasing its mature and
SC
activated form [63]. In the present study, the expression of genes related to apoptosis focused on both the inhibition and promotion of apoptosis. Basing from the present
M AN U
data, we conclude that the balanced promotion and inhibition of apoptotic factors ensure the stability of S. constricta in adapting to an environment laden with heavy
TE D
metals, such as Cd and Pb.
Acknowledgments
This work was financially supported through the Public welfare project in Zhejiang
EP
Province (2015C32004), Natural Science Foundation of Ningbo (2015C10009, 2015C10062), Collaborative Innovation Center for Zhejiang Marine High-Efficiency
AC C
and Healthy Aquaculture, and the K. C. Wong Magna Fund at Ningbo University.
References
[1] Yuan XC, Zhao WW. China Fishery Statistical Yearbook. Fisheries Agency of China Agriculture Ministry. 2015; 56. [2] Wang SK, Jin BS, Qin HM, Sheng Q, Wu JH. Trophic Dynamics of Filter Feeding Bivalves in the Yangtze Estuarine Intertidal Marsh: Stable Isotope and Fatty Acid Analyses. Plos One. 2015; 10(8): e0135604.
ACCEPTED MANUSCRIPT [3] Jin S, Xue CB, Wang GL, Lu TX, Wang YN, Chen YE. Bacterial flora composition and its dynamics in tidal-flat Sinonovacula constricta aquaculture area. Chinese journal of applied ecology. 2005; 16(7): 1380-2. [4] Lu ZB, Du Q, Fang MJ, Qian XM, Cai QH, Xu CY. Carrying capacity of shellfish
16(5):961-6.
RI PT
culture in Dadeng Island sea area of Xiamen. 2005; Ying Yong Sheng Tai Xue Bao.
[5] Niu DH, Wang F, Xie SM, Sun FY, Wang Z, Peng MX et al. Developmental Analysis
and
Identification
of
Genes
Involved
SC
Transcriptome
in
Larval
Metamorphosis of the Razor Clam, Sinonovacula constricta. Marine Biotechnology.
M AN U
2016; 18(2): 168-75.
[6] Feng BB, Dong LL, Niu DH, Meng SS, Zhang B, Liu DB, et al. Identification of Immune Genes of the Agamaki Clam (Sinonovacula constricta) by Sequencing and Bioinformatic Analysis of ESTs. Marine Biotechnology. 2010; 12(3): 182-91.
TE D
[7] Niu DH, Xie SM, Bai ZY, Wang L, Jin K, Li JL. Identification, expression, and responses to bacterial challenge of the cathepsin C gene from the razor clam
241-5.
EP
Sinonovacula constricta. Developmental and Comparative Immunology. 2014; 46(2):
[8] Ma HT, Jiang HB, Liu XQ, Wu XP, Wei XM. Polymorphic microsatellite loci for
AC C
the razor clam, Sinonovacula constricta. Genetics and Molecular Research. 2015; 14(1): 145-8.
[9] Yuan Y, Li Q, Yu H, Kong LF. The Complete Mitochondrial Genomes of Six Heterodont Bivalves (Tellinoidea and Solenoidea): Variable Gene Arrangements and Phylogenetic Implications. PLoS One. 2012; 7(2): e32353. [10] Lu YL, Zhang AG, Li CH, Zhang P, Su XR, Li Y, et al. The link between selenium binding protein from Sinonovacula constricta and environmental pollutions exposure. Fish and Shellfish Immunology. 2013; 35(2): 271-7.
ACCEPTED MANUSCRIPT [11] Zhang AG, Lu YL, Li CH, Zhang P, Su XR, Li Y et al. A small heat shock protein (sHSP) from Sinonovacula constricta against heavy metals stresses. Fish and Shellfish Immunology. 2013; 34(6): 1605-10. [12] Natalie D, Jill MY, Chiu, Beverly HK, Pob, Anna CK et al. Heavy metal
RI PT
contamination along the China coastline: A comprehensive study using Artificial Mussels and native mussels. Journal of Environmental Management. 2016; 180(9): 238-46.
SC
[13] Tran D, Fournier E, Durrieu G, Massabuau JC. Copper detection in the Asiatic
clam Corbicila fluminea: optimum valve closure response. Toxicology. 2004; 66(3):
M AN U
333-43.
[14] Haberkorn H, Lambert C, Le Goic N, Moal J, Palacios E, Lassus P et al. Effects of Alexandrium minutum exposure upon physiological and hematological variables of diploid and triploid oysters, Crassostrea gigas. Aquatic Toxicology. 2010; 2(97): 96108.
[15] Lin W, Rice MA, Chien PK. The differential effects of three heavy metals on
TE D
particle filtration and amino acid uptake in the pacific oyster, Crassostrea gigas. Journal of Shellfish Research. 12: 111.
[16] AbeT, Yamamoto O,Gotoh S, Yan Y, Todaka Y, Higashi K. Cadmium-Induced
EP
mRNA Expression of Hsp32 Is Augmented in Metallothionein-I and -II Knock-out
AC C
Mice. Archives of Biochemistry and Biophysics. 2000; 1(382): 81-8. [17] Zhang PJ, Li CH, Zhu L, Su XR, Li Y, Jin CH, et al. De Novo Assembly of the Sea Cucumber Apostichopus japonicus Hemocytes Transcriptome to Identify miRNA Targets Associated with Skin Ulceration Syndrome. PLoS One. 2013; 8(9): e73506. [18] Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary toolfor transcriptomics. Nature Review Genetics. 2009; 10(1): 57-63. [19] Niu DH, Wang L, Sun FY, Liu ZJ, Li JL. Development of Molecular Resources for an Intertidal Clam, Sinonovacula constricta, Using 454 Transcriptome Sequencing.
ACCEPTED MANUSCRIPT PLoS One. 2013; 8(7): e67456. [20] Moreira R, Pereiro P, Canchaya C, Posada D, Figueras A, Novoa B. RNA-Seq in Mytilus galloprovincialis: comparative transcriptomics and expression profiles among different tissues. BMC Genomics. 2015; 16(1): 728.
RI PT
[21] Schroeder A, Mueller O, Stocker S, Salowsky R, Leiber M, Gassmann M, et al. The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Molecular Biology. 2006; 7:3.
SC
[22] Wang, L. et al. The promotion of cytoskeleton integration and redox in the haemocyte of shrimp Litopenaeus vannamei after the successive stimulation of
M AN U
recombinant VP28. Developmental and Comparative Immunology. 2014; 45: 123-32. [23] Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Fulllength transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology. 2011; 29(7): 644-52.
TE D
[24] Rice P, Longden I, Bleasby A. EMBOSS: the European Molecular Biology Open Software Suite. Trends in Genetics. 2000; 16(6): 276-7. [25] Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and
621-8.
EP
quantifying mammalian transcriptomes by RNA-Seq. Nature Methods. 2008; 5(7):
AC C
[26] Wang LK, Feng ZX, Wang X, Wang XW, Zhang XG. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010; 26(1):136-8.
[27] Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2005; 57: 289-300. [28] Altschul SF, Madden TL, Schaffer AA, et al. Gapped BLAST and PSI BLAST: A New Generation of Protein Database Search Programs. Nucleic Acids Research.
ACCEPTED MANUSCRIPT 1997; 25(17): 3389-402. [29] Apweiler R, Bairoch A, Wu CH, et al. UniProt: the Universal Protein knowledgebase. Nucleic Acids Research. 2004; 32(Database issue): D115-9. [30] Deng YY, Li JQ, Wu SF, et al. Integrated nr Database in Protein Annotation
RI PT
System and Its Localization. Computer Engineering. 2006; 32(5): 71-4.
[31] Ashburner M, Ball CA, Blake JA, et al. Gene ontology: a tool for the unification of biology. Nature genetics. 2000; 25(1): 25-9.
SC
[32] Tatusov RL, Galperin MY, Natale DA. The COG database: a tool for genome scale annalysis of protein functions and evolution. Nucleic Acids Research. 2000;
M AN U
28(1): 33-6.
[33] Koonin EV, Fedorova ND, Jackson JD, et al. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biology. 2004; 5(2): R7.
TE D
[34] Kanehisa M, Goto S, Kawashima S, et al. The KEGG resource for deciphering the genome. Nucleic Acids Research. 2004; 32(Database issue): D277-D280. [35] Berger VJ, Kharazova AD. Mechanisms of salinity adaptation in marine
EP
molluscs. Hydrobiologia. 1997; 355: 115-26. [36] Lockwood BL, Somero GN. Transcriptomic responses to salinity stress in
AC C
invasive and native blue mussels (genus Mytilus). Molecular Ecology. 2011; 20(3): 517-29.
[37] Bergman HL, Dorward-King EJ. Reassessment of metals criteria for aquatic life protection: priorities for research and implementation. Pensacola: SETAC Press. 1997; 31-56. [38] Trautinger F, Kindås-Mügge I, Barlan B, Neuner P, Knobler RM. 72-kD heat shock protein is a mediator of resistance to ultraviolet B light. Journal of Investigative
ACCEPTED MANUSCRIPT Dermatology. 1995; 105(2): 160-2. [39] Lee S, Prochaska DJ, Fang F, Barnum SR. A 16.6-kilodalton protein in the Cyanobacterium synechocystis sp. PCC 6803 plays a role in the heat shock response. Current Microbiology. 1998; 37(6): 403-7.
RI PT
[40] Köhler HR, Triebskorn R, Stöcker W, Kloetzel PM, Alberti G. The 70 kD heat shock protein (HSP 70) in soil invertebrates: a possible tool for monitoring environmental toxicants. Archives of environmental contamination and toxicology.
SC
1992; 22(4): 334-8.
[41] Joshi CP, Nguyen HT. 5' untranslated leader sequences of eukaryotic mRNAs
M AN U
encoding heat shock induced proteins. Nucleic Acids Research. 1995; 23(4): 541-9. [42] Chen Q, Osteryoung K, Vierling E. A 21-kDa chloroplast heat shock protein assembles into high molecular weight complexes in vivo and in Organelle. The Journal of Biological Chemistry. 1994; 269(18): 1216-23.
TE D
[43] Kusumoto T, Maehara Y, Sakaguchi Y, Emi Y, Kohnoe S, Sugimachi K. 1Hexylcarbamoyl-5-fluorouracil alters the expression of heat shock protein in HeLa cells. Anticancer Drugs. 1991; 2(1): 45-8.
EP
[44] Riordan M, Sreedharan R, Wang S, Thulin G, Mann A, Stankewich M, et al. HSP70 binding modulates detachment of Na-K-ATPase following energy deprivation
AC C
in renal epithelial cells. American Journal of Physiology Renal Physiology. 2005; 288(6): F1236-42.
[45] Kroese LJ, Scheffer PG. 8-hydroxy-2'-deoxyguanosine and cardiovascular disease: a systematic review. Current Atherosclerosis Reports. 2014; 16(11): 452. [46] Adams TK, Saydam N, Steiner F, Schaffner W, Freedman JH. Activation of gene expression by metal-responsive signal transduction pathways. Environ Health Perspect. 2002; 110 Suppl 5: 813-7. [47] Dorian C, Gattone VH, Klaassen CD. Accumulation and degradation of the
ACCEPTED MANUSCRIPT protein moiety of cadmium-metallothionein (CdMT) in the mouse kidney. Toxicology and applied pharmacology. 1992; 117(2): 242-8. [48] Dalton T, Fu K, Enders GC, Palmiter RD, Andrews GK. Analysis of the effects of overexpression of metallothionein-I in transgenic mice on the reproductive
RI PT
toxicology of cadmium. Environmental Health Perspectives. 1996; 104(1): 68-76.
[49] Snyder RD. A review and investigation into the mechanistic basis of the genotoxicity of antihistamines. Mutation Research. 1998; 411(3): 235-48.
SC
[50] Abdel-Mageed AB, Agrawal KC. Activation of nuclear factor kappaB: potential role in metallothionein-mediated mitogenic response. Cancer Research. 1998; 58(11):
M AN U
2335-8.
[51] Lahiri A, Abraham C. Activation of pattern recognition receptors up-regulates metallothioneins, thereby increasing intracellular accumulation of zinc, autophagy, and bacterial clearance by macrophages. Gastroenterology. 147(4): 835-46.
TE D
[52] Galdes A, Hill HA, Kägi JH, Vasak M, Bremner I, Young BW. The 1H N.M.R. spectra of metallothioneins. Experientia Suppl. 1979; 34:241-8. [53] Tsujikawa K, Suzuki N, Shimaoka T, Iida K, Kohama Y, Otaki N, et
EP
al.Abnormal accumulation of copper-metallothionein in the hepatopancreas and kidney of Long-Evans rats with a cinnamon-like coat color (LEC rats).Biological &
AC C
pharmaceutical bulletin. 1994; 17(5): 591-5. [54] Zhao XL, Yu H, Kong LF, Li Q. Transcriptomic responses to salinity stress in the Pacific oyster Crassostrea gigas. PLoS One. 2012; 7(9): e46244. [55] Liu J, LiuY, MichalskaAE, ChooKH, KlaasenCD. Metallothionein plays less of a protective role in cadmium-metallothionein-induced nephrotoxicity than in cadmium chloride-induced hepatotoxicity. The Journal of Pharmacology. 1996; 276(3): 121623.
ACCEPTED MANUSCRIPT [56] Gazdag Z, Fujs S, Koszegi B, Kálmán N, Papp G, Emri T, et al. The abc1-/coq8respiratory-deficient mutant of Schizosaccharomyces pombe suffers from glutathione underproduction and hyperaccumulates Cd2+. Folia Microbiologica. 2011; 56(4): 353-
RI PT
9. [57] Dourado DF, Fernandes PA, Mannervik B, Ramos MJ. Glutathione transferase: new model for glutathione activation. Chemistry. 2008; 14(31): 9591-8.
SC
[58] Tan KH, Meyer DJ, Gillies N, Ketterer B. Detoxification of DNA hydroperoxide by glutathione transferases and the purification and characterization of glutathione
841-5.
M AN U
transferases of the rat hepatopancreas nucleus. Biochemical Journal. 1988; 254(3):
[59] Mazzullo M, Colacci A, Grilli S, Prodi G, Arfellini G. 1,1,2-Trichloroethane: evidence of genotoxicity from short-term tests. Japanese journal of cancer research : Gann. 1986; 77(6): 532-9.
TE D
[60] Murphy BJ, Andrews GK, Bittel D, Discher DJ, McCue J, Green CJ, et al. Activation of metallothionein gene expression by hypoxia involves metal response elements and metal transcription factor-1. Cancer Research. 1999; 59(6): 1315-22.
EP
[61] Eckelman BP, Salvesen GS, Scott FL. Human inhibitor of apoptosis proteins:
AC C
why XIAP is the black sheep of the family.EMBO Reports. 2006; 7(10): 988-94. [62] Richter BWM, Duckett CS. The IAP proteins: caspases inhibitors and beyond. Science STKE. 2000; 44: 1-4. [63] Ko DC, Gamazon ER, Shukla KP, Pfuetzner RA, Whittington D, Holden TD, et al. Functional genetic screen of human diversity reveals that a methionine salvage enzyme regulates inflammatory cell death. Proceedings of the National Academy of Sciences of the United States of America. 2012; 109(35): E2343-52.
ACCEPTED MANUSCRIPT
Figure legends
RI PT
Fig 1. Summary statistics of the counts of the different expression genes in gill (A) and hepatopancreas (B).
Fig 2. Classification of the annotated amino-acid sequences. Amino-acid sequences were grouped into different functional sub-categories: cellular component, molecular
SC
function and biological process.
Fig 3. COG function classification consensus sequence.
M AN U
Fig 4. KOG function classification consensus sequence.
Fig 5. Summary statistics of the counts of genes involving in different processes according to GO analysis. A: counts of genes down-regulated in gill 12 h vs control; B: counts of genes up-regulated in gill 12 h vs control; C: counts of genes downregulated in gill 48 h vs control; D: counts of genes up-regulated in gill 48 h vs control; E: counts of genes down-regulated in hepatopancreas 12 h vs control; F:
TE D
counts of genes up-regulated in hepatopancreas 12 h vs control; G: counts of genes down-regulated in hepatopancreas 48 h vs control; H: counts of genes up-regulated in hepatopancreas 48 h vs control.
Fig 6. Validation of relative expression level of GSH, MT and HSP32 by qRT-PCR.
EP
A: Relative expression level of GSH in gill; B: Relative expression level of MT in gill; C: Relative expression level of HSP32 in gill; D: Relative expression level of GSH in
AC C
hepatopancreas; E: Relative expression levelof MT in hepatopancreas; F: Relative expression level of HSP32 in hepatopancreas.
ACCEPTED MANUSCRIPT
Table 1. Summary statistics of the transcriptome sequencing for S. constricta from
Gill controlled
(12h)
(48h)
Gill treated (48h)
Hepatopancreas
Hepatopan-
Hepatopancreas
controlled
creas treated
Hepatopancreas controlled
(12h)
(12h)
(48h)
(48h)
SC
(12h)
Gill treated
M AN U
Gill controlled
RI PT
control group (PC) and treated group (PT)
treated
7,205,607
9,861,573 10,301,319
9,578,956
7,311,534
22,115,620
8,406,059
13,666,784
clean reads
3,291,760
4,189,228
4,346,584
4,375,850
3,618,758
11,665,491
3,476,846
6,720,966
clean reads (%)
45.68%
42.48%
42.19%
45.68%
49.49%
52.57%
41.36%
49.18%
mapped reads
4,616,632
6,685,160
mapped reads (%)
64.07%
67.79%
TE D
Raw sequencing reads
6,091,258
4,921,393
15,107,180
5,786,731
9,613,215
65.29%
63.95%
67.31%
68.31%
68.84%
70.34%
AC C
EP
6,725,731
ACCEPTED MANUSCRIPT Table 2. Summary statistics of the number of unigenes mapping all annotated databases
Unigene
≥300nt
≥1000nt
COG_Annotation
4,798
2,436
2,362
GO_Annotation
6,173
3,539
KEGG_Annotation
9,314
5,325
KOG_Annotation
11,244
6,408
Pfam_Annotation
12,971
7,107
Swissprot_Annotation
8,964
nr_Annotation
17,770
All_Annotated
18,330
SC
M AN U
TE D EP AC C
RI PT
Annotated databases
2,634 3,989 4,836
5,864
4,848
4,116
10,585
7,185
10,984
7,346
ACCEPTED MANUSCRIPT Table 3. Real-time quantitative PCR primer information in this study
Names
Sequences(5′-3′)
Application
ScHSP32
GGCATCTGGATTGGAACTATGATAG
Real-time PCR
CTTACCCCAACTGCATAAATGTTTC GTGTATGAGTCCGTGGTGGTGTGTA GAATGCAGGTCCTCTCGTATCTCAT ScMT
ATGAGGCATGTATAGACAGACAGTG
ScActin
CGGGATCCATGGCTGAGACAATG
EP
TE D
M AN U
CCCAAGCTTCTAGCTAAGGAGTTTCTGGTC
AC C
Real-time PCR
SC
GTCTCATTAAGGTCTCCATGGTAAC
Real-time PCR
RI PT
ScGSH
Real-time PCR
ACCEPTED MANUSCRIPT Table 4. Detailed annotations of differentially expressed genes associated with cadmium detoxification
Gene name
Highest absolute value log2(fold_change)
P-value
UN000905
HSP32 △
10.14
0.0001036
UN014045
MT △
5.6
0.0002495
Expression tissue/time
Gene Ontology
Up/Down
gill/48 h
F: catalytic activity (GO:0003824); P: metabolic process (GO:0008152); F: metal ion binding (GO:0046872)
up
P: response to metal ion (GO:0010038)
up
F: catalytic activity (GO:0003824)
up
gill/48 h
F: ATP binding (GO:0005524); P: response to stress (GO:0006950)
up
gill/24 h
P: response to stress (GO:0006950); F: unfolded protein binding (GO:0051082)
down
P: intracellular protein transport (GO:0006886); F: protein transporter activity (GO:0008565);P:vesic le-mediated transport (GO:0016192)
up
F: metal ion binding (GO:0046872)
down
F:sodium-potassiumexchani-ng ATPase activity (GO:0005391); C: integral component of membrane (GO:0016021); F: metal ion binding (GO:0046872)
down
SC
Hepatopancreas/48 h
UN106383
HSP70
2.76
HSP90
-1.78
0.000252 6
0.0003581
0.0004574
AP-1
2.34
AC C
UN068669
UN017282
Zn/Cu superoxide dismutase
Na+-K+ATPase
Hepatopa-
ncreas/48 h
Hepatopa-
0.0003513 ncreas/24 h
EP
UN066133
3.4
M AN U
UN081663
GSH △
TE D
UN005612
Hepatopa-2.50
0.0005741 ncreas/24 h
-1.91
RI PT
Transcript ID
0.0004097
gill/48 h
ACCEPTED MANUSCRIPT
MREs
3.2
0.0001533 ncreas/48 h
HepatopaUN017157
MTF-1
2.8
0.0003507
M AN U TE D EP AC C
ion
F: metal ion binding (GO:0046872); F: organic cyclic compound binding (GO:0097159); F: heterocyclic compound binding (GO:1901363)
SC
ncreas/48 h
F: manganese binding (GO:0030145)
RI PT
HepatopaUN009810
up
up
ACCEPTED MANUSCRIPT
Figure 1
SC
RI PT
A
AC C
EP
TE D
M AN U
B
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 2
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 3
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 4
ACCEPTED MANUSCRIPT
Figure 5
SC
RI PT
A
M AN U
B
AC C
C
EP
TE D
E
ACCEPTED MANUSCRIPT
SC
RI PT
D
EP AC C
F
TE D
M AN U
E
ACCEPTED MANUSCRIPT
SC
RI PT
G
AC C
EP
TE D
M AN U
H
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Figure 6
ACCEPTED MANUSCRIPT
Highlights: Eight cDNA libraries were constructed and sequenced in Cd2+ exposed Sinonovacula constricta.
RI PT
Total of 18,330 unigenes with an average length of 500 bp were annotated. Different expression unigenes were identified in gill and liver at two examined time points.
SC
ROS production-related genes were synchronously enriched in all experimental
AC C
EP
TE D
M AN U
samples