RNA-Sequencing Analysis of TCDD-Induced Responses in Zebrafish Liver Reveals High Relatedness to In Vivo Mammalian Models and Conserved Biological Pathways

TCDD is one of the most persistent environmental toxicants in biological systems and its effect through aryl hydrocarbon receptor (AhR) has been well characterized. However, the information on TCDD-induced toxicity in other molecular pathways is rather limited. To fully understand molecular toxicity of TCDD in an in vivo animal model, adult zebrafish were exposed to TCDD at 10 nM for 96 h and the livers were sampled for RNA-sequencing based transcriptomic profiling. A total of 1,058 differently expressed genes were identified based on fold-change>2 and TPM (transcripts per million) >10. Among the top 20 up-regulated genes, 10 novel responsive genes were identified and verified by RT-qPCR analysis on independent samples. Transcriptomic analysis indicated several deregulated pathways associated with cell cycle, endocrine disruptors, signal transduction and immune systems. Comparative analyses of TCDD-induced transcriptomic changes between fish and mammalian models revealed that proteomic pathway is consistently up-regulated while calcium signaling pathway and several immune-related pathways are generally down-regulated. Finally, our study also suggested that zebrafish model showed greater similarity to in vivo mammalian models than in vitro models. Our study indicated that the zebrafish is a valuable in vivo model in toxicogenomic analyses for understanding molecular toxicity of environmental toxicants relevant to human health. The expression profiles associated with TCDD could be useful for monitoring environmental dioxin and dioxin-like contamination.


Introduction
Dioxin-like compounds are major environmental contaminants that could pose serious threats to public health and the ecosystem [1]. Since 1980s, it has been increasingly documented that dioxinlike compounds cause various biological effects in laboratory animals and human [2,3]. Among them, 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) is the most potent toxicant and it is produced from both natural and anthropogenic processes including incineration of chlorine-containing substances, bleaching of paper, manufacturing of specific organochlorine chemicals, volcanoes, and forest fires [4]. As an aromatic hydrocarbon, TCDD has a long biological half-life and is heavily accumulated in the food chain, which causes adverse effects on human health at environmental levels [4,5]. Till now, many animal models, including fishes, have been used in TCDD studies [6,7]. However, most of these studies have focused on the physiological-biochemical parameters and typical molecular markers, such as the gene cyp1a and the Aryl hydrocarbon receptor (AhR) signaling pathway [7] and detailed molecular toxicity of TCDD remains to be elucidated.
Genomic approaches have been increasingly used in toxicological research in the past decade. Previous toxicogenomic studies mostly used DNA microarray technology to capture the global gene expression data and to evaluate the effects of toxicant exposure [8,9]. However, there are several drawbacks in microarray analysis, e.g. limited sensitivity, probe cross-hybridization, incomplete genome coverage and a prerequisite for sequence information in order to include new probes [10]. Recently, the advent of next-generation sequencing (NGS) technologies has significantly accelerated genomic research and provided a better alternative for transcriptomic analysis [11]. By high-throughput RNA sequencing, it is feasible to measure transcript abundance and transcriptomic profiles with a broad dynamic range, therefore providing a powerful tool to determine the potential adverse effects of environmental contaminants on public health [8].
Comparative studies across different taxonomic groups are important not only in understanding of organism diversity but also for inferring important biological responses due to evolutionary conservation [4]. In toxicology, animal models are widely used to infer human responses to chemical exposure. Particularly, it is valuable to use multiple animal models for comparative analyses in order to determine conserved adverse responses or molecular events. As lower vertebrates in evolution, fishes are particularly important models for such comparative studies [12,13]. Now the zebrafish (Danio rerio) has become an increasingly popular model in human disease studies because of its many advantageous properties in laboratory experiments, such as its small size and easy availability in a large number, short generation time for genetic manipulation, and cost effectiveness for high-throughput studies. More importantly, many molecular and developmental studies have shown that zebrafish and human share many common genes in conserved developmental pathways in organogenesis and related physiological processes as well as in carcinogenesis [14,15].
As there is so far no NGS based transcriptomic analysis for molecular response to TCDD treatment, in the present study we carried out RNA sequencing analyses of TCDD-treated in order to carry out a genome-wide identification of novel TCDD responsive genes and pathways, which should provide a comprehensive understanding of the molecular mechanism of TCDD-induced toxicity in mammals and human. We further compared the zebrafish response with those from mammalian systems and thus identified common molecular pathways deregulated in both fish and mammals. We found that the zebrafish model is more similar to mammalian in vivo models than in vitro models, thus indicating a validity of the zebrafish as an emerging in vivo model in comparative toxicological research.

Zebrafish
Experimental procedures were carried out following the approved protocol by Institutional Animal Care and Use Committee of National University of Singapore (Protocol 079/07). Adult zebrafish (3-month old), which are Singapore wild type zebrafish, were purchased from a local aquarium farm (Mainland Tropical Fish Farm, Singapore) and acclimated for at least one week in our aquarium before chemical exposure experiment. Fish were maintained at ambient temperature of around 28uC with a 14-h light and 10-h dark cycle in a flow-through water system.

Chemical exposure
TCDD was purchased from Sigma-Aldrich (USA) and dissolved in dimethyl sulfoxide (DMSO). Male adult fish were used in the exposure experiments by immersing in the TCDD water (10 nM) for 96 h at ambient temperature (28uC) in a static condition. Control fish were kept in water with 0.01% DMSO (vehicle) under the same condition. Water was changed daily throughout the treatment. After 96 hours of chemical exposure, treated and control fish were sacrificed and liver samples were dissected from each fish.

RNA sample preparation and SAGE library sequencing
Total RNA was extracted from livers (excluding gall bladders) of individual fish using TRIzolH Reagent (Invitrogen) and treated with DNase I (Invitrogen) to remove genomic DNA contamination. For RNA sequencing, RNA was pooled equally from 9 fish for each group (TCDD and control). Poly A+ RNA was purified using DynabeadsH Oligo (dT) EcoP (Invitrogen) and subjected to cDNA synthesis. Synthesized cDNA was digested by NlaIII and sequencing adapters were added to the cDNA fragments. SAGE (serial analysis of gene expression) sequencing (tag length = 27 nucleotides) was performed by Mission Biotech (Taiwan) with ABI SOLiD TM System 2.0 (Applied Biosystems). The RNA-sep data reported in the present study was submitted to Gene Expression Omnibu with an access number GSE49915.

Gene annotation and selection of differentially expressed genes
All SAGE tags were mapped to the zebrafish Reference Sequence database (http://www.ncbi.nlm.nih.gov/RefSeq) with maximum 2 nucleotide mismatch. Uniquely mapped tag counts for each transcript were normalized to TPM (transcripts/tags per million). For biological implication analyses, genes with only marginal expression, as defined by TPM,10 in both control and TCDD groups, were excluded. As it has been previously reported that the actual measured quantity of differential expression (fold change or ratio) is more consistent and reproducible in identifying differentially expressed genes than the statistical significance (pvalue) [16,17], in the present study, we selected differentially expressed genes based on fold change.2 and TPM.10.

Gene Ontology enrichment analysis
Gene ontology enrichment analysis was performed using DAVID (The Database for Annotation, Visualization and Integrated Discovery) with the total zebrafish genome information as the background and p-values based on a modified Fisher's exact t-test. Gene Ontology Fat categories were used for this analysis and the cut-off p-value is 0.05.

Real-time PCR
Real-time PCR was performed using the LightCycler system (Roche Applied Science) with LightCycler FastStart DNA Master SYBR Green I (Roche Applied Science) according to the manufacturer's instruction. For comparison between real-time Analysis of the hepatic enriched gene list by GSEA and IPA analysis GSEA (Gene Set Enrichment Analysis) pre-ranked option was used to analyze the entire set of differentially expressed genes (651 up-and 407 down-regulated genes). Briefly, the gene symbols of human homologs of the enriched zebrafish Unigene clusters were ranked using logarithm transformed fold change (base 2). The number of permutation used was 1000. Pathways with false discovery rate (FDR) ,0.25 were considered statistically significant. For IPA (Ingenuity Pathway Analysis), the same set of differentially expressed genes was uploaded to online Ingenuity Pathways Knowledge Base for functional implication analyses.

Cross-species comparison
Six sets of transcriptomic data for in vivo and in vitro mammalian models treated by TCDD (GSE10082, GSE10083, GSE10769, GSE10770, GSE14555 and GSE34251) were retrieved from GEO (Gene Expression Omnibus). GSEA was used to establish the relatedness between zebrafish and mammalian models. The zebrafish hepatic transcriptome lists were converted into human and mouse homolog Unigene clusters. The statistical significance of the enrichment score was estimated by using an empirical phenotype-based permutation test. An FDR value was provided by introducing adjustment of multiple hypothesis testing.

Results and Discussion
General features of SAGE tags in TCDD-treated and control groups There were 11.7 million and 17.9 million SAGE sequence tags generated from the DMSO vehicle control and TCDD-treated groups, respectively. About 2.7 million tags in the control group and 3.2 million tags in the TCDD group could be mapped with the known transcripts. The mapped tags were normalized to transcripts per million (TPM) and the expression level of genes ranged from 0.3 to 110844.2 TPM in the two groups, indicating a dynamic range of more than six orders of magnitude in transcript abundance. As shown in Figure 1A, transcript abundance profiles in both groups were very similar. The transcriptomes consist of a small number of high abundant transcripts and a large number of low abundant transcripts, similar to those reported in many previous RNA-seq studies [18,19]. The percentages of the transcripts with relatively high expression level (TPM.256) were 3.8% in the TCDD group and 3.1% in the control group, but contributed to 80.4% and 83.0% of the total transcripts respectively. In contrast, the percentages of the transcripts with the low expression level (TPM,16) were 67.1% and 76.3% in the TCDD group and control group, contributing only 3.1% and 3.2% of the total transcripts.
Differentially expressed genes in response to TCDD and their gene ontology profile Figure 1B shows the relationship between the dynamic transcript fold-change of hepatic transcriptome and its expression level in zebrafish treated by TCDD. Differentially expressed genes after TCDD treatment were first determined by comparison of the two sets of mapped SAGE tags with fold-change.2 and TPM .10. In total, 1,058 genes were identified, including 651 upregulated genes and 407 down-regulated genes. The top 20 most up-regulated transcripts based on fold-change are listed in the Table 1. After compared with TCDD responsive genes in CTD (The Comparative Toxicogenomics Database, http://ctdbase. org/), ten gens (cyp1a, ahrra, actc1b, odf3b, mid1ip1, pklr, lrrc23, ccdc113, sdf2, and fgf13) have been reported to be related in fish and/or mammalian models after TCDD exposure among the top 20 genes. Not surprisingly, cyp1a, the best known molecular marker for TCDD exposure [20], was the most up-regulated gene (151.4 fold) in the list, indicating that our TCDD treatment in the experiment was effective. It is interesting to note that adrra (aryl hydrocarbon receptor repressor a) is the second highest upregulated gene, while the two aryl hydrocarbon receptor genes (ahr1a and ahr2) were not significantly up-regulated in our RNAseq data with a modest 50% increase for adr2 and slightly decrease for adr1a (data not shown), 10 novel TCDD-responsive genes, eight annotated (rab36, sycp3l, rnasel2, tnni2b.2, hormad1, tuba7l, irak1bp1 and zgc:152916,) and two unannotated (zgc:193690 and si:dkey-190g11.7), have been identified from the top 20 most upregulated genes. To confirm their inducibility of the 20 genes by TCDD, real-time PCR was carried out with individual fish liver samples from an independent set of TCDD treatment experiment. As shown in Figure 2, 17 out of 20 genes showed up-regulation in at least three individual samples; among them, 13 genes had upregulation in all of the five individual samples. Thus, our RT-qPCR data indicated a strong agreement with the RNAsequencing data and, more importantly, the validation was from five independent biological samples. Interestingly, all of the 10 novel TCDD responsive genes found in the present study were validated by RT-qPCR ( Figure 2). Thus, our RNA-seq data provided additional biomarker genes for TCDD exposure. Among these validated biomarker genes, several of them apparently encodes secreted protein (e.g. fgf13b and rnasel2) and their protein products are also likely present in the circulating blood and may offer convenient non-invasive assays for detection of TCDD exposure. A notable exception in the validation experimental data was lrrc23 that was constantly down-regulated in all five individuals, indicating that irrc23 is not a reproducibly upregulated gene by TCDD. In future, these novel molecular markers for TCDD exposure can be further tested for their timeand dose-responsiveness. Besides cyp1a, some other genes with high abundance were also up-regulated by TCDD exposure in zebrafish (Table S1). The upregulation of pklr (pyruvate kinase, liver and red blood cell) in our experiment was consistent with several previous reports that TCDD affects pyruvate utilization for energy and thus glycolysis and gluconeogenesis [21,22]. Gene expression of alas1, the first and rate-limiting enzyme involved in heme biosynthesis [23], was upregulated in our study, suggesting that the heme biosynthesis was deregulated in hepatic cells of fish exposed to TCDD. Hspa5 is a member of the heat shock protein 70 (Hsp70) family and it has been used as an endoplasmic reticulum (ER) stress sensor [24]. In our study, the expression of hspa5 was up-regulated with 10.2 fold, implying TCDD led to the ER stress and this event was a highly conserved adaptive response induced by environmental toxicants. Moreover, the gene pck1 (phosphoenolpyruvate carboxykinase 1), as a key factor in gluconeogenesis, was up-regulated 6.8 fold, indicating that TCDD exposure also affected the gluconeogenesis by altering the expression of genes encoding key gluconeogenic enzymes [25].
Through gene ontology analysis, the differentially expressed genes were classified into different categories based on GO database using DAVID software ( Table 2). Among the up-regulated GO categories, the most significant ones included Endoplasmic reticulum (p-value = 1.19E-13) and Proteasome complex (p-value = 2.67E-09). Since most P450s (e.g. CYP1A enzymes) are involved in xenobiotic metabolism and primarily located in endoplasmic reticulum, it is not surprising that the TCDD-induced physiological stress occurred mainly in these cellular components, which was also supported by the up-regulated 19 genes of CYP family in our study (Table S2). Moreover, other histological studies also found TCDD led to hypertrophy of hepatocytes, glycogen depletion and lipidosis of liver in zebrafish and other fish [26], which were related to the endoplasmic reticulum (ER) stress via some cellular receptor or/and protein [27,28]. In subsequent analysis, our results also showed that the proteasome related pathways (Proteasome pathway and HSA03050 Proteasome, see Table 3) were significantly up-regulated in zebrafish after TCDD treatment. Consistent with this, several GO categories involved in Proteasome deregulation, such as Protein folding and Cellular protein catabolic process in the Biological Process category, as well as Unfolded protein binding in the Molecular Function category, further indicating that the proteasome related bio-functions were deregulated.

Change of transcription factor networks by TCDD exposure
By IPA analysis, six significant transcription factors were found to be related to TCDD-induced hepatotoxicity (Table 4), including Xbp1, Nfe2l2, Nr5a2, Ptf1a, Tp53 and Mycn. In particular, the two up-regulated Xbp1 and Nfe2l2 networks were highly enriched with P-values of 1.47E-16 and 4.73E-11 respectively, In the Xbp1 network, 38 target genes were found from the differentially expressed gene list and almost all of them (36) were up-regulated ( Figure 3A), indicating that Xbp1 plays a central role in regulating a battery of genes responsible for protein trafficking and secretion [29]. In the Nfe2l2 network, 36 up-and 14 down-regulated genes were induced in the differentially expressed gene list ( Figure 3B). This is in consistence with a previous study [30] that most of the regulated genes (such as mgst1, usp14, herpud1, dnajc3, actg1, hmox1, dnajb11, etc.) by Nfe2l2 are related to oxidative stress and mediate transcriptional events that facilitate protective responses in animal models exposed to xenobiotic.   (Table S1)

Deregulated pathways in zebrafish livers treated by TCDD
To investigate the change of molecular pathways induced by TCDD treatment, GSEA was performed using the set of differentially expressed genes (651 up and 407 down) and 163 up-and 123 down-regulated pathways were identified (Table S3). Among all pathways, five up-and eight down-regulated pathways had significant FDR values less than 0.25 and are shown in Table 4. The first three pathways are associated with cell cycle progression, which have also been reported in mammalian models under TCDD stress [31,32]. However in the present study, cell cycle and related pathways were significantly enhanced while some mammalian studies showed a decrease of these activities, suggesting different regulatory mechanism in various models as well as by different treatment regimes.
The mechanisms responsible for TCDD toxicity are always associated with its ability to disrupt endocrine functions. In our data, ubiquitin-proteasome pathway was significantly up-regulated in fish liver after TCDD treatment. Consistent with another study in the mammalian system [33], the up-regulated proteasome and related pathways in our study further indicated that the TCDDinduced ubiquitin-proteasomal degradation of AhR influenced the nucleus transcription by controlling the level of ligand-activated AhR. As ubiquitin-proteasome pathway is significantly induced in zebrafish and all mammalian models analyzed in our study, it is an apparently conserved mechanism of TCDD-induced toxicity between fish and mammals. Prostaglandin synthesis regulation pathway was suppressed in TCDD-treated zebrafish, consistent with a previous report that prostanoid synthesis pathway could be regulated by COX2-TBXS-TP (cyclooxygenase2-thromboxane A synthase1-thromboxane receptor) in zebrafish after TCDD treatment [34].
In the present study, several pathways related to signal transduction were significantly inhibited in zebrafish after TCDD treatment, including Ca 2+ regulation pathways, Fas signaling pathway, T-cell signal transduction and T-cell receptor signaling pathway, ribosomal protein and its related pathways. TCDD has been reported to significantly increase intracellular free Ca 2+ in many cell culture systems [35] and to trigger Ca 2+ -mediated endonuclease activity leading to apoptosis [36]. Our study, consistent with several previous reports [37][38][39][40], indicated that Ca 2+ regulation pathways could be commonly involved in TCDD action. Furthermore, Fas signaling pathway, T-cell signal transduction and T-cell receptor signaling pathway were also suppressed in our study. Fas, a member of the tumor necrosis factor receptor (TNFR) family, contains a death domain which is essential for the delivery of the death signal [41]. Thus, Fas signaling regulation could be a mechanism of impaired T-cell related pathway induced by TCDD, which are consistent with another previous report [41]. Moreover, ribosomal proteins, in conjunction with rRNA, are involved in the cellular process of translation [42]. Interestingly, the ribosomal protein and its related pathway were down-regulated in zebrafish liver after TCDD exposure, but were induced in other mammalian models [43].
Furthermore, three immune-related pathways were significantly repressed in the zebrafish liver after TCDD exposure (Table 4), including natural killer cell mediated cytotoxicity, T-cell signal transduction and T-cell receptor signaling pathway. Moreover, several B-cell related pathways were also found to be inhibited with FDR.0.25 (Table S3). Collectively, our data indicate that hepatic immune-related functions were impaired in fish exposed to TCDD.

Toxicogenomic comparison between zebrafish and mammalian models
In order to gain insight into the common molecular toxicity of TCDD between fish and mammals as well as the validity of the zebrafish model to predict chemical toxicity for risk assessment for human health, our transcriptomic data were compared with available transcriptomic data from both in vivo and in vitro mammalian studies with TCDD. Among 26 related series in the GEO database using human cell lines, rat and mouse tissues and cells based on microarray studies, we found six series (GSE10082, GSE10083, GSE10769, GSE10770, GSE14555 and GSE34251) with effective comparability (Table 5) while other series were not included in the comparative analysis due to incompatibility in data uploading, platform, experimental strategy, etc. The list of 650 upregulated genes from the current zebrafish study was used to represent the zebrafish transcriptome in GSEA. As shown in Figure 4A, based on normalized enrichment scores (NES) and false discovery rate (FDR), the zebrafish hepatic genes expression showed more resemblance to the in vivo models (average NES is 1.54) than in vitro models (average NES is 0.93). The comparison with all of the four in vivo data, FDR was smaller than 0.25, while none of the in vitro data showed such significance. These observations indicate that the zebrafish data is more similar to the in vivo mammalian data than in vitro data, further enforcing the validity of the zebrafish model as a potentially high-throughput and economic in vivo experimental models for studies relevant to human health.
To further analyze the correlation of zebrafish and other mammalian models after TCDD treatment, comparison of pathways were made by GSEA pre-ranked function. One upregulated (Proteasome pathway) and three down-regulated pathways (HSA04650 Natural killer cell mediated cytotoxicity, HSA04660 T cell receptor signaling pathway, HSA04020  Calcium signaling pathway) in zebrafish showed the same regulation direction as those in all other mammalian models, indicating that the four pathways were well conserved in vertebrates. Meanwhile, HSA03010 Ribosome pathway displayed completely opposite regulated trends between zebrafish and all mammalian models we compared ( Figure 4B), indicating the species-specific patterns between fish and mammals after TCDD treatment, which need further analysis in future. Moreover, after comparing the conserved pathways (the same regulation direction pathways) among zebrafish model, in vivo and in vitro mammalian models, we noted that the degree of regulation is more consistent between the zebrafish and in vivo mammalian models than between the zebrafish and in vitro mammalian models, as indicated by the color codes in Figure 4B. In summary, by RNA-sequencing based transcriptomic analyses, we have carried out detailed analyses of TCDD-induced molecular changes in zebrafish. Other than the well characterized AhR pathway and cyp1a1 biomarker gene induced by TCDD, we also found some new biomarker genes that have been validated from independent experiments. Interestingly, some of the new biomarker genes encode secreted proteins such as Fgf13b, Sdf2 and RNasel2, which may by analyzed from serum samples by a non-invasive approach. Further and comparative transcriptomic analyses of our zebrafish data and available mammalian data from TCDD treatment experiments indicate several well conserved TCDD-responsive pathways, including up-regulated proteomic pathway and several down regulated pathways such as calcium signaling pathway, natural killer cell mediated cytotoxicity, T-cell receptor signaling pathway etc. Furthermore, GSEA analyses indicate that the TCDD-induced zebrafish transcriptomic data is more similar to in vivo mammalian data than in vitro data, thus indicating the validity of the zebrafish model as a valuable in vivo model to infer molecular toxicity relevant to human health.

Supporting Information
Table S1 Complete list of differentially expressed genes in zebrafish liver exposed to TCDD with a cut off of fold change.2 and TPM .10.

(XLSX)
Table S2 Up-regulated CYP family genes in zebrafish treated by TCDD.

(XLSX)
Table S3 Complete list of pathways in zebrafish exposed to TCDD as revealed by GSEA analysis. (XLSX)