Reduced Zebrafish Transcriptome Atlas toward Understanding

May 21, 2018 - Transcriptomic approaches monitoring gene responses at genome-scale are increasingly used in toxicological research and help to clarify...
0 downloads 0 Views 6MB Size
Download PDF
Article Cite This: Environ. Sci. Technol. 2018, 52, 7120−7130

pubs.acs.org/est

Reduced Zebrafish Transcriptome Atlas toward Understanding Environmental Neurotoxicants Kun Zhang†,‡ and Yanbin Zhao*,†,‡ †

School of Environmental Science and Engineering, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China Brigham and Women’s Hospital, Harvard Medical School, 60 Fenwood Road, Boston, Massachusetts 02115, United States

‡

Environ. Sci. Technol. 2018.52:7120-7130. Downloaded from pubs.acs.org by UNIV OF SOUTH DAKOTA on 08/15/18. For personal use only.

S Supporting Information *

ABSTRACT: Transcriptomic approaches monitoring gene responses at genome-scale are increasingly used in toxicological research and help to clarify the molecular mechanisms of adverse effects caused by environmental toxicants. However, their applications for chemical assessment are hampered due to high expenses required and more importantly the lack of indepth data mining and mechanistic perspectives. Here, we described a reduced transcriptome atlas (RTA) approach which integrates transcriptomic data sets and a comprehensive panel of genes generated to represent neurogenesis and the early neuronal development of zebrafish, to determine the potential neurodevelopmental toxicities of environmental chemicals. Transcriptomic data sets of 74 chemicals and 736 related gene expression profiles were integrated resulting in 135 exposure signatures. Chemical prioritization demonstrated four sets of hits to be neurotoxic: neuro-active chemicals (representatively, Valproic acid, VPA and Carbamazepine, CAR), xenoestrogens (Bisphenol A, BPA; Genistein, GEN; 17-α ethinylestradiol, EE2), microcystins (cyanopeptolin, CP1020; microcystin-LR, MCLR) and heavy metals (AgNO3, AgNPs). The enriched biological pathways and processes were distinct among the four sets, while the overlapping functional enrichments were observed within each set, for example, over 25% differentially expressed genes and four of top five KEGG pathways were shared between VPA and CAR. Furthermore, gene expression index (GEI) analysis demonstrated that a gene panel with 300 genes was sufficient to effectively characterize and cluster chemicals and therefore offer an efficient and cost-effective tool for the prioritization of neurotoxicants. Thus, the RTA approach provides novel insights into the understanding of the in-depth molecular mechanisms of environmental neurotoxicants and can be used as an indication for potential adverse outcomes.

■

INTRODUCTION Transcriptomics technologies can characterize and quantify the entirety of gene transcripts simultaneously and consequently are useful to identify changes in molecular patterns related to environmental stressors, such as toxic chemicals. Transcriptomics have increasingly been applied in toxicological research. They help to clarify the molecular mechanisms of adverse effects and predict the potential toxic effects caused by environmental chemicals.1−3 The integration of transcriptomics and high-throughput screening contribute to the adverse outcome pathway analysis for chemical risk assessments4,5 and assist with the classification and prioritization in large toxicity screening program, for instance, Tox21 and the human toxome project.6,7 High throughput transcriptomic approaches have received more and more attention recently, however, several issues still limit their extensive utilizations in general toxicology. One major concern that significantly hampered the progress is the lack of data mining tools and advanced knowledge discovery.8 Due to the data complexity, most of the transcriptomics-based toxicology studies still focus on the signaling pathways and © 2018 American Chemical Society

biomarker genes related to the already known adverse effects of the toxicants, which helps to explain the molecular mode of actions of the toxic outcomes observed at higher organizational levels (e.g., the histological and physiological toxicities),9,10 whereas other meaningful molecular responses were either overlooked or ignored as not directly relevant to the expected mode of actions. For example, bisphenol A (BPA), with weak estrogenic activities, has been classified as an endocrine disruptor and a “substance of very high concern”.11 Transcriptomics highlighted its dysregulations on endocrine systems such as on the receptor binding and estrogen stimulus, steroidogenesis and cell proliferation.12,13 However, the significant changes were also observed for other pathways, such as the DNA metabolic process, oxidative stress, apoptosis, neurogenesis, and nervous system development. These dysregulations indicate that BPA would have unanticipated Received: Revised: Accepted: Published: 7120

March 12, 2018 May 17, 2018 May 21, 2018 May 21, 2018 DOI: 10.1021/acs.est.8b01350 Environ. Sci. Technol. 2018, 52, 7120−7130

Article

Environmental Science & Technology

Figure 1. Construction of the reduced transcriptome atlas (RTA) approach in zebrafish. (A) Design and the workflow of the RTA approach development. (B) L2000 gene set represents the extrinsic and intrinsic factors involved in neurogenesis and early neuronal development of zebrafish.

regulators/targets rather than the whole transcriptome.19 Especially, the LINCS (NIH library of integrated networkbased cellular signatures) program as well as the related L1000 expression profiling assay has generated a large-scale public data set that indicates how human cells respond to various genetic and environmental stressors, which already included over 27 000 perturbagens and produced over 1.3 million L1000 profiles.20,21 As a powerful and effective tool, the LINCS program greatly accelerated the drug discovery. For toxicology, a human S1500 sentinel gene set were generated from 3339 gene expression series by Tox21 program (Toxicology Testing in the 21st Century) recently with full pathway coverages. It can be used to accurately predict pathway perturbations and biological relationships for samples under study.22 Similarly, a reduced human transcriptome (RHT) approach and a reduced zebrafish transcriptome (RZT) approach were sequentially developed via the integration of 1200 human Entrez genes and 1637 zebrafish Entrez genes, respectively, to comprehensively represent the biological pathways and toxicologically relevant processes.23,24 They were proven to be promising for profiling the molecular pathways that modulated by reference chemicals and the water samples ranging from wastewater to drinking water.23,24 The enriched pathways in RHT and RZT analysis could be used to evaluate and prioritize chemicals for future assessment. Specifically, the most sensitive biological pathways could be linked with adverse outcomes at higher organizational levels and could identify the early molecular responses for adverse outcome pathway (AOP) construction and assess the hazard potencies of environmental toxicants in early developmental stages.23,24 On this basis, the aim of present study was to develop a reduced transcriptome atlas (RTA) by integrating a large number of transcriptomic data sets with a specific gene panel (L2000) to comprehensively represent the neurogenesis and early neuronal development in zebrafish. This approach was aimed to uncover and prioritize the potential neurodevelopmental toxicants by in-depth transcriptomic data mining and interpretation. Meanwhile, a specific gene set was aimed to be generated to effectively characterize and prioritize chemicals as L2000 and therefore serves to establish a cost-effective reduced zebrafish transcriptome approach that specifically targets the nervous system.

toxicities on DNA repair process and zebrafish neurodevelopment that are beyond the adverse outcomes on the endocrine system. Meanwhile, due to the high biological complexity, many adverse outcomes in animals under chemical exposure cannot be easily detected. The lack of advanced knowledge and technologies to qualitatively and quantitatively characterize the traits of different phenotypes hampered the progress to characterize and prioritize toxic chemicals. Especially for the nervous system, there are limited strategies that were commonly applied to assess and prioritize the neuronal disruptive chemicals during neurogenesis and early neuronal development, for example, by immunohistochemistry, in situ hybridization and the emerging behavioral testing (e.g., locomotor activities in teleost).14,15 Transcriptomics approaches can help to predict these unanticipated toxic effects, however, an effective data interpretation is still hard to perform due to the limited knowledge of gene-function relationships, regulatory mechanisms and the immature data mining tools.8,14 For instance, functional enrichment analysis of transcriptome can clarify the dysregulated biological pathways and processes, such as the altered mRNA metabolic process, cellar protein catabolic process and chromosome organization caused by silver ions.9 Nevertheless, limited perspectives can be obtained based on these general biochemical alterations and it is still difficult to uncover the related toxic outcomes at higher organizational levels. Thus, the challenges on data analysis and interpretation call for novel methodological strategies, especially for the complex nervous system where the toxicity outcomes were hard to determine qualitatively and quantitatively. Another obstacle confronting the utilization of transcriptomics is the relatively high cost. In spite of falling prices, transcriptomic profiling remains expensive to perform, especially for a large scale toxicological research, such as the high throughput screening for chemical clustering and prioritization. In recent years, alternative strategies were developed with substantially reduced costs.16−19 Of them, a reduced transcriptome approach was increasingly used. It includes a small number of genes as a gene panel which can effectively represent the whole-genome expressions, such as the L1000 landmark genes18 and 917 human pathway reporter gene panel.19 In this way, the activity of entire signaling networks can be assessed based on these established key 7121

DOI: 10.1021/acs.est.8b01350 Environ. Sci. Technol. 2018, 52, 7120−7130

Article

Environmental Science & Technology

■

MATERIALS AND METHODS RTA Gene Panel Development. The RTA gene panel was aimed to comprehensively represent the biological pathways and processes involved in neurogenesis and early neuronal development in zebrafish (Danio rerio). Thus, we propose three criteria to get a maximum coverage in the present study: Representing the entire neural networks, maximal coverage of biological pathways and processes and integrating novel neural factors (Figure 1). Our current knowledge and large-scale signaling pathway databases provide us the opportunity to retrieve most of the currently known neural factors involved in neural networks and signaling pathways. With a comprehensive survey of the up-to-date literatures, novel annotated neural genes can be further integrated. Therefore, at first, we constructed a general framework for the neural signaling pathways and related genes based on the current knowledge,25−27 which included most of the classic extrinsic and intrinsic factors, such as neurotransmitters and growth factors. Then, the associated genes within each signaling pathway were curated and integrated from two databases: Kyoto Encyclopedia of Genes and Genomes (KEGG) database and Molecular Signatures Database (MSigDB). For the pathways generated for other species, like Homo sapiens and Mus musculus, the corresponding zebrafish homologous genes were retrieved from Ensembl zebrafish genome database (http://useast. ensembl.org/Danio_rerio/Info/Index). The Gene Ontology (GO) data set from the European Bioinformatics Institute (EMBL-EBI) and MSigDB helps to uncover more neural factors. Nervous system development (GO biological process, GO: 0007399) and the child GO terms, such as neurogenesis (GO: 0022008) and neural tube development (GO:0021915), provide a large set of genes involved in neurogenesis and early neuronal development. This data set was further reviewed and clustered. The neural genes beyond the classic pathways described above were further integrated manually. Meanwhile, the gene sets from Mammalian Neurogenesis Gene Ontology (MANGO), a database of genes mapped to adult neurogenesis via a comprehensive survey of the literatures,28 were also checked to include some novel annotated neural genes. At last, all genes were integrated and clustered manually and the related detailed annotations were extracted from Ensembl database and Genebank database using R (version 3.4.3). This information was stored in FASTA format. In total, the gene panel in the present study contained 2015 genes (termed L2000) (Supporting Information (SI) Table S1), which were divided into six groups: extrinsic morphogens, growth factors, neurotransmitters and synaptic vesicle cycle, intrinsic transcription factors, cytoskeletal proteins and epigenetic regulators. At last, the representation of RTA gene panel on neural development was further validated by the enrichment analysis of KEGG_Pathway and GO_Biological Process on the online bioinformatics resources database: The Database for Annotation, Visualization and Integrated Discovery (DAVID). Data Sets Collection. The basic criteria for data inclusion were as follows: (I) Exposure durations representing the neurogenesis processes/early neuronal development. In zebrafish embryos, the brain ventricles have formed by 48 h postfertilization (hpf). Neurogenesis peaks between 18 hpf and 36 hpf.29 Therefore, the exposures covered the 0−48 hpf neurogenesis or exposures after the first 48 hpf to represent the early neuronal development were selected for the present study.

(II) Studies have appropriate experimental setups. Treatment groups and vehicle controls should have sufficient replicates for the statistical analysis. (III) Complete raw/normalized data set with sufficient annotation. The raw data or normalized data set should be completed and can be retrieved from public repositories or related literature. Annotation repository should be supplied with sufficient information to cover all gene probes (Figure 1). The criteria for data inclusion basically followed the approaches proposed by Ramasamy et al.30 and the MIAMEGuidelines.31 A database query was then conducted in the publicly accessible NCBI Gene Expression Omnibus (GEO) and EMBL-EBI ArrayExpress online repositories. Data sets satisfying the above restrictions were manually selected and integrated. Besides, the PubMed literature database was also employed to search for relevant studies. Data set was included if the criteria met. In total, the searching resulted in 27 studies, which contain 74 environmental chemicals and 736 related gene expression profiles (SI Table S1). Data Normalization and Analysis. The raw data and platform information were downloaded from GEO or Array Express. For data normalization, if available, the processed signal intensity was employed. The raw data were log2 transformed and then the sample qualities were checked using R package arrayQualityMetrics. The differentially expressed gene (DEG) was identified with a moderated t test using R-package “limma”. Cutoff criterium for DEGs was an adjusted P-value of 3, that is, xenoestrogens (EE2, GEN), microcystins (MCLR), a neuro-drug (VPA) and heavy metals (AgNO3, AgNPs) (Figure 2b). Hierarchical clustering map depicts 135 exposure signatures for the L2000 gene set (Figure 2c). Of these signatures, only a small part displayed clear similar gene expression patterns. They can be roughly divided into three sets: (a) the highest prioritization occurred for pharmaceuticals, that is, VPA, CAR, retinoic acid (RA), and caffeine (CAFF). They are all neurodrugs or neuroactive compounds. (b) the second prioritized group consisted of extrinsic environmental hormones, including the classic xenoestrogens, like EE2, GEN, BPA, 17β-estradiol (E2), and others such as pentachlorophenol, triclosan, propiconazole, and diisobutyl_phthalate. (c) the third group consisted of microcystins and heavy metals, including CP1020, AgNO3, and the related silver nanoparticles (Ag50, Ag150, and PVP) (Figure 2d). Therefore, the chemical prioritization analysis provides a list of hit chemicals that displayed remarkable dysregulations on genes involved in neurogenesis and early neuronal development in zebrafish embryos. The most frequently occurred substances were: neuro-drugs and neuroactive compounds (representa7124

DOI: 10.1021/acs.est.8b01350 Environ. Sci. Technol. 2018, 52, 7120−7130

Article

Environmental Science & Technology

Figure 3. Functional enrichment analysis of the 12 hit chemicals. (A) KEGG pathway enrichment. Top five KEGG pathways are displayed for each of the 12 chemicals (valproic acid (VPA), carbamazepine (CAR), fluoxetine (FLO), bisphenol A (BPA), 17-α ethinylestradiol (EE2), genistein (GEN), cyanopeptolin (CP1020), microcystin-LR (MCLR), AgNO3, Ag150, cobalt and 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD)). (B) Heatmap depicts the transcriptional responses of key extrinsic signaling (morphogens, growth factors, neurotransmitters and synaptic vesicle cycle) under four types of chemicals exposure (VPA, BPA, MCLR, and Ag150). (C) Polar area diagram shows the enriched neural signaling pathways in zebrafish embryos under VPA exposure. The top 10 enriched pathways are described on the right. The percentage value shows the proportion of different expressed genes compared to the related gene panel for each signaling pathway recruited for the present study.

enriched and identical Wnt signaling and TGF-beta signaling were observed for CP1020 and MCLR and the enriched and identical Wnt signaling, VEGF signaling and Adrenergic signaling were observed for AgNO3 and AgNP-150. In comparison, the GO_ Biological Processes were varied among the hit chemicals (SI Figure S2). It should be noted that Wnt signaling was the most enriched pathway and biological process among the hit chemicals, which was not surprising because it is indispensable in virtually every aspect of embryonic development.47 However, the in-depth analysis of L2000 displayed quite distinct dysregulations on Wnt signaling components as well as other pathways among different types of chemicals. The transcriptional changes of the representative extrinsic signaling (i.e., morphogens, growth factors, neurotransmitters and synaptic vesicle cycle) response to four types of chemicals (VPA, BPA, MCLR, and Ag150) is depicted in Figure 3b. For VPA (Figure 3b), almost all signals were remarkably dysregulated; the ratio ranges from 54.2% (13/24, BMP signaling) to 86.4% (19/22, FGF signaling). For example, 71.1% (59/83) genes involved in Wnt signaling were significantly changed, with 31 genes up-regulated (representatively, sfrp2, wif1, and nlk2) and 28 genes down-regulated (representatively, apc2, wnt4b, gsk3b, and ctnnd2b). This result indicated that VPA had comprehensive adverse effects on the signals of neuronal development in zebrafish embryos. In contrast, BPA had a specific enrichment on the extrinsic signals. Compared to the enriched Wnt signal (15.6%, 13/83), the enrichment on neurotransmitters (Glutamate signaling, 48.9%, 22/45) and synapse vesicle cycle process (52.3%, 46/88) was much stronger. For MCLR, it displayed less effect on most of the extrinsic signals such as growth factors and neurotransmitters, the weak enrichments were only found for Wnt signaling (26.5%, 22/83) and synapse vesicle cycle (22.7%, 22/

actions and the expression patterns during the hierarchical clustering. On the contrary, the exposures with relatively high concentrations (such as pharmacological or lethal concentration) generally did well. This phenomenon was quite consistent with the observations in a similar report previously on zebrafish transcriptome meta-analysis.32 Thus, to be able to derive biologically meaningful information on different substances via transcriptome analysis, only the exposures with reported visible effects, sublethal concentrations or pharmacological concentrations were included in the following functional enrichment analysis, as described above. Functional Enrichment Analysis. To systemically characterize the hit chemicals, the DEGs were subjected to DAVID for KEGG pathway and Gene Ontology enrichment analysis. The top five KEGG pathways and GO Biological Processes enriched for each of the 12 hit chemicals (VPA, CAR, Fluoxetine (FLO), BPA, EE2, GEN, CP1020, MC-LR, AgNO3, Ag150, Cobalt, and TCDD) were depicted (Figure 3a, SI Figure S2). VPA and CAR, the commonly administered antiepileptic drugs, showed very similar dysregulations on the KEGG pathways, four of the top five pathways were identical between them. They were Wnt signaling, GnRH signaling, MAPK signaling and Melanogenesis. In contrast, the antidepressant drug of the selective serotonin reuptake inhibitor (SSRI) class, FLO, displayed a distinct dysregulation. Only the weak dysregulated Wnt signaling and GnRH signaling were identical with others. BPA and EE2 also showed similar patterns, with three of the top five pathways identical, that is, Wnt signaling, GnRH signaling, and Notch signaling. Interestingly, the enriched pathways of BPA and EE2 also resembled the neuro-drugs (VPA and CAR). On the contrary, the microcystins and heavy metal forms showed quite different patterns compared to the neuro-drugs and xenoestrogens. The 7125

DOI: 10.1021/acs.est.8b01350 Environ. Sci. Technol. 2018, 52, 7120−7130

Article

Environmental Science & Technology

Figure 4. Overlapping functional enrichments of the neuro-active chemicals. (A) Venn diagram shows the overlap of different expressed RTA genes relative to control among the four neuro-active chemicals (VPA, CAR, RA, and CAFF). (B) The top four shared KEGG pathways are shown at below. (C) Overlapped gene ontology terms (BP: Biological process; CC: cellular component; MF: molecular function) among the four neuroactive chemicals shown in a bubble plot using the Goplot package. Top five GO terms with adjusted p-value