Zebrafish Whole-Adult-Organism Chemogenomics for Large-Scale Predictive and Discovery Chemical Biology

doi:10.1371/journal.pgen.1000121

PLoS Genet. 2008 July; 4(7): e1000121.

Published online 2008 July 11. doi: 10.1371/journal.pgen.1000121

PMCID: PMC2442223

Copyright Lam et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Zebrafish Whole-Adult-Organism Chemogenomics for Large-Scale Predictive and Discovery Chemical Biology

Siew Hong Lam,¹ Sinnakarupan Mathavan,² Yan Tong,¹ Haixia Li,² R. Krishna Murthy Karuturi,² Yilian Wu,¹ Vinsensius B. Vega,² Edison T. Liu,² and Zhiyuan Gong¹^*

¹Department of Biological Sciences, National University of Singapore, Singapore

²Genome Institute of Singapore, Agency for Science Technology and Research, Singapore

Stuart K. Kim, Editor

Stanford University Medical Center, United States of America

* E-mail: dbsgzy/at/nus.edu.sg

Conceived and designed the experiments: SHL RK VBV ETL ZG. Performed the experiments: SHL YT HL YW. Analyzed the data: SHL SM YT HL RK YW VBV ZG. Contributed reagents/materials/analysis tools: SHL SM YT HL RK YW VBV ETL ZG. Wrote the paper: SHL SM RK ETL ZG.

Received February 3, 2008; Accepted June 9, 2008.

This article has been cited by other articles in PMC.

Abstract

The ability to perform large-scale, expression-based chemogenomics on whole adult organisms, as in invertebrate models (worm and fly), is highly desirable for a vertebrate model but its feasibility and potential has not been demonstrated. We performed expression-based chemogenomics on the whole adult organism of a vertebrate model, the zebrafish, and demonstrated its potential for large-scale predictive and discovery chemical biology. Focusing on two classes of compounds with wide implications to human health, polycyclic (halogenated) aromatic hydrocarbons [P(H)AHs] and estrogenic compounds (ECs), we generated robust prediction models that can discriminate compounds of the same class from those of different classes in two large independent experiments. The robust expression signatures led to the identification of biomarkers for potent aryl hydrocarbon receptor (AHR) and estrogen receptor (ER) agonists, respectively, and were validated in multiple targeted tissues. Knowledge-based data mining of human homologs of zebrafish genes revealed highly conserved chemical-induced biological responses/effects, health risks, and novel biological insights associated with AHR and ER that could be inferred to humans. Thus, our study presents an effective, high-throughput strategy of capturing molecular snapshots of chemical-induced biological states of a whole adult vertebrate that provides information on biomarkers of effects, deregulated signaling pathways, and possible affected biological functions, perturbed physiological systems, and increased health risks. These findings place zebrafish in a strategic position to bridge the wide gap between cell-based and rodent models in chemogenomics research and applications, especially in preclinical drug discovery and toxicology.

Author Summary

To understand chemical-induced biological responses/effects, it is important to have large-scale and rapid capacity to investigate gene expression changes caused by chemical compounds at genome-wide scale in an adult vertebrate model; this capability is essential for drug development and toxicology. Small aquarium fish with vast genomic resources, such as zebrafish, will probably be the only vertebrate models that allow for cost-effective, large-scale, genome-wide determination of gene expression net changes in the entire adult organism in response to a chemical compound. Presently, such a whole adult organism approach is only feasible in invertebrate models such as the worm and fly, and not in rodent models, hence the usefulness of such an approach has not been demonstrated in a vertebrate. By using two classes of chemicals with wide implications to human health, we showed that capturing net changes of gene expression at a genome-wide scale in an entire adult zebrafish is useful for predicting toxicity and chemical classes, for discovering biomarkers and major signaling pathways, as well as for inferring human health risk and new biological insights. Our study provides a new approach for genome-wide investigation of chemical-induced biological responses/effects in a whole adult vertebrate that can benefit the drug discovery process and chemical toxicity testing for environmental health risk inference.

Introduction

Chemogenomics, application of genomic tools in pharmacology and toxicology, offers a promising approach that will enhance drug discovery (target identification/validation, lead identification, efficacy evaluation) and toxicity assessment [1],[2]. Presently, invertebrates such as the worm Caenorhabditis elegans and fly Drosophila melanogaster, are the only animal models that have benefited from whole-adult-organism expression chemogenomics [3]–[6]. The benefits of whole-adult-organism chemogenomics would usually translate into large-scale, high-throughput, high-content and cost-effective applications for chemical biology. It is highly desirable that the benefits of whole-adult-organism chemogenomics can be realized in a vertebrate model because of the many biological processes, health risks and diseases that are restricted to a mature vertebrate system including humans. The existing cell-, fly- and worm-based models, while suited for high-throughput chemogenomics, lacked the relevant physiological whole-organism setting of an adult vertebrate. This is especially important in the context of pharmacology and toxicology when many of the potentially targeted organ-systems such as the endocrine, digestive (liver in particular), immune, muscular-skeletal, vasculature, kidney are absent from the existing high-throughput models. In contrast, the rodent models, though providing in vivo adult vertebrate data, are not suited for high-throughput applications and are not cost-effective [7], thus creating a bottleneck situation when in vivo biological data, especially toxicology, is required for the high number of ‘hits’ generated from in vitro screenings or for the many newly emerging industrial compounds and waste that are coming into contact with the public and environment. We propose that whole-adult chemogenomics performed on a small vertebrate such as the zebrafish would be a strategy that is sufficiently high-throughput, cost-effective and would generate high content in vivo vertebrate data potentially useful for large-scale screening and toxicity testing purposes.

Conceptually, whole-adult-organism expression chemogenomics would capture the sum-total of the transcriptomic changes in an entire adult organism as a single biological entity responding to exogenous chemical cues. This, however, would have its inherent limitations such as loss of weak signals or signals from smaller tissues and loss of specific location of response, and they may be compounded further by the greater biological complexity in vertebrates compared to invertebrates. Thus, while whole-adult-organism chemogenomics had been shown to be useful in invertebrate models with regard to compound screening [3],[4] and identifying biological processes affected by specific compounds [5],[6], it is not known if chemogenomics data generated from a whole adult vertebrate will be useful. We hypothesized that since strong and well-represented expression signals are likely to be detected in whole-adult-organism chemogenomics, the expression signals that are captured would be robust for predictive chemical biology and for uncovering biology that is widely associated with the chemical-induced responses/effects in the adult vertebrate.

However, the idea of performing high-throughput whole-adult-organism chemogenomics on a vertebrate model was practically not feasible, if not unimaginable, until microarray technology was made available to small aquarium fish such as the zebrafish. The availability of the zebrafish in large numbers, its small size, low husbandry cost, vast genomic resources and its present use in disease modeling [8] and drug screening [9],[10], make the zebrafish ideal for high-throughput whole-adult-organism chemogenomics. Moreover, owing to their close physiological relationship with the environment, fish are highly sensitive to environmental changes particularly exogenous chemical cues; therefore the impact of chemical effects on fish system is more easily defined and readily studied than on terrestrial species [11]. Previously, we and others have shown that zebrafish responded biologically to chemicals, such as small molecules, drugs and environmental toxicants, in a similar manner as mammals [12]–[16]. In this study, we chose P(H)AHs [represented by Benzo[a]pyrene (BAP), 3-Methylcholanthrene (MC), 2,3,7,8-Tetrachlorodibenzodioxin (TCDD)] and ECs [represented by 17-beta estradiol (Es), Diethylstilbestrol (DES), Bisphenol A (Bis)] as model compounds because they represent two classes of chemicals with wide implications to human health. Both P(H)AHs and ECs are potent AHR and ER agonists, respectively, and these receptors are known to cross-talk and they are regulators of important cellular functions that are involved in various biological processes and have been associated with several patho-physiological conditions [17]–[19]. Some of these compounds have been used as drugs or investigated for therapeutic potential [20]–[22]. Moreover, both classes of compounds are also environmental carcinogens and endocrine disruptors that have generated considerable public health concern [23],[24]. By focusing on P(H)AHs and ECs in this study, we performed chemogenomics on whole adult zebrafish and demonstrated that it is good for large-scale predictive chemical biology, for discovering biomarkers and major signaling pathways, as well as useful for human health risk and biological insight inference. Our study placed zebrafish in a strategic position to bridge the gap between in vitro cell-based model and in vivo rodent model in chemogenomics.

Results/Discussion

Robust Predictive Power of Zebrafish Whole-Adult-Organism Chemogenomics

We generated 159 samples/arrays involving 28 treatment groups (4–7 replicates in each treatment group) in two experiments (‘A’ and ‘B’) and grouped them into four Datasets: I and II from experiment ‘A’ while III and IV from experiment ‘B’ (See Methods, Table S1 and Figure S1). Experiments ‘A’ and ‘B’ were performed one year apart using different batches of fish, reference RNA, reagents, array prints and experimental designs, which will test the robustness of the prediction models and help in identifying robust biomarkers. We trained six prediction models using Dataset I (Figure 1A and 1B) and validated them independently on Datasets II, III and IV (Figure 1C–1E). The validations were ‘independent’ in the sense that Datasets II, III and IV were totally left out and not used in the training of the prediction models that were tested on. Next, we trained another six prediction models using Dataset III (Figure 2A and 2B) and performed similar independent validation on the ‘unseen’ Datasets I, II and IV (Figure 2C–2E). The prediction models were trained using two supervised learning classifiers, k-nearest neighbours (kNN) and support vector machine (SVM), on selected discriminatory gene sets from Dataset I or III. During the training phase, the supervised learning classifiers used the individual discriminatory gene sets together with their expression data in Dataset I or III to produce respective sets of rules or reference weights that will serve as standards/models for the prediction of ‘unseen’ or unknown samples. The discriminatory gene sets for P(H)AH and EC classes were selected using threshold criteria of Q-value, FDR-value or P-value coupled with fold-difference between treated [combined representative groups of P(H)AHs or ECs] versus control samples in Dataset I or Dataset III (Figure 1A and and2A;2A; details of criteria in Table S2). The different statistical treatments and learning classifiers were used to examine if any idiosyncrasy associated with data processing affected the performances of the prediction models.

Figure 1

Performance of P(H)AH and EC prediction models trained using discriminatory gene sets from Dataset I.

Figure 2

Performance of P(H)AH and EC prediction models trained using discriminatory gene sets from Dataset III.

First, we evaluated the prediction models using ‘leave-one-out’ cross-validation approach to avoid the statistical problem of over-estimating prediction accuracy that occurs when a model is trained and evaluated with the same samples. In this procedure, one of the 30 or 25 samples from Datasets I or III, respectively, was withheld and the remaining 29 or 24 of the respective samples were used to train a prediction model based on a selected discriminatory gene set to predict the class of the withheld sample. The process was repeated until all 30 or 25 samples were predicted in turn, and all prediction models were tested similarly. We found that the ‘Predicted Results’ for all the samples matched the ‘Expected Results’, hence 100% of the samples were correctly classified by all the prediction models (Figure 1B and and2B).2B). The excellent performances of the prediction models were non-random (Fisher's Exact Test P-value

8.35×10⁻⁹−4.89×10⁻⁷; Table S3), indicating their predictive powers.

To test the robustness of their predictive powers, we sought to validate them independently. In this procedure, we used the 30 samples from Dataset I to train prediction models based on the selected discriminatory gene sets and predict each of the 25 ‘unseen’ samples from Dataset III. Remarkably, the prediction models for P(H)AHs and ECs trained by Dataset I classified 100% of the samples in Dataset III correctly (Figure 1C). We then reversed the order by training the prediction models with the 25 samples from Dataset III and testing them on each of the 30 samples from Dataset I. The prediction models for P(H)AHs trained by Dataset III classified all the samples in Dataset I correctly while the prediction models for ECs performed comparably well though with one or two false negatives (Figure 2C). The findings indicate that the prediction models were remarkably robust because both Datasets I and III were obtained through two separate experiments that contained substantial biological and technical variations.

Next, we tested the specificity (the ability to identify negative cases) of the above prediction models on 9 types of compounds from classes with acute mode-of-action and toxicity effects differing from that of P(H)AHs and ECs. Thus, we tested each of the prediction models generated using the 30 and 25 samples from the respective Datasets I and III, on each of the 46 ‘expected-negative’ samples from Dataset II, independently. As anticipated, with only one exception, all samples were predicted as ‘negative’ matching the ‘Expected Results’ (Figure 1D and and2D).2D). The robust performance of the prediction models generated from Dataset III is noteworthy since both Datasets II and III were derived from different experiments.

We also tested the performance of the prediction models on Dataset IV consisting of 58 samples obtained from fish exposed to multiple chemical mixture of BAP, DES and arsenic (As) at different concentrations and combinations. Arsenic was introduced to increase the complexity of the mixtures and to test if the prediction models could still perform well on samples exposed to mixtures consisting of a chemical not used in the training of the models. The performance of the prediction models for ECs was outstanding as they could classify 100% of the samples correctly notwithstanding the increase complexity that could arise from the mixture of compounds and that Datasets I and IV are from separate experiments (Figure 1E and and2E).2E). The prediction models for P(H)AHs displayed comparable performances in terms of specificity although the sensitivity (ability to identify positive cases) varied from 62.5%–97.5% (Table S3). Notably, prediction models of P(H)AHs trained from Dataset I performed poorer on Dataset IV compared to those trained from Dataset III, suggesting that when mixtures were involved, the performance of the P(H)AH models were affected by the inter-experimental variations. In this case, different statistical approaches can affect prediction performance, but with appropriate statistical tests, it is possible to generate prediction models that are sufficiently robust. This was observed in the case of SVM-trained P(H)AH models using discriminatory gene sets from Dataset I selected based on FDR-val (sensitivity

92.5%) and P-val (sensitivity

90.0%) which performed better compared to Q-val (sensitivity

62.5%) (Table S3).

Taken together, with the exception of Q-val_P(H)AH performance on Dataset IV, all the prediction models performed comparably well, in particular for most of ECs which scored 100% for specificity and sensitivity. The robust ability of the prediction models to discriminate compounds of the same class from those of different classes, even in a mixture of compounds, suggests that zebrafish whole-adult-organism chemogenomics is capturing genes associated with biological functions that are strongly affected by a class of compound. This demonstrates its potential use for compound screening, predictive chemical biology and biomarker discovery. The ability of the prediction models to tolerate a reasonable level of biological, technical and data-processing variations, indicate high amenability to real-life compound screening and predictive applications as such variations will inevitably occur over time, in different laboratories and experimental settings.

Identification of Potential Biomarker Genes and Signaling Pathways

Having demonstrated their predictive powers, we used the discriminatory gene sets to identify potential biomarker genes for P(H)AHs and ECs. To do so, we first consolidated the discriminatory gene sets into their corresponding two major groups P(H)AHs and ECs by combining the gene sets within their respective classes (including only genes with unique GenBank Identity and similar mean expression directionality). Then, we examined the consistency of their expression profiles throughout all the 28 treatment groups in this study and validated some of the responsive genes using quantitative real-time PCR. A two-way hierachical clustering showed that the consolidated gene sets were able to cluster tightly the respective treatment groups including those in the mixture groups from the non-P(H)AH (Figure 3A) or non-EC treatment groups (Figure 3B). However, the gene expression profiles formed in the EC gene set (Figure 3B) were more distinct compared to the P(H)AH gene set (Figure 3A). This is due to the large number of genes that are specifically affected by ECs compared to the presence of certain xenobiotic metabolism or stress response-associated genes shared between P(H)AHs and other compounds, as well as the presence of some genes whose expressions are affected by compound mixtures. A closer examination of the genes reveals that ahr2 and its known regulated/responsive genes [25] such as cyp1A1, NQO1 homolog, nfe2l2, TIPARP homolog, gstp1, cyp1C1 were among the genes found in a tight cluster that shows similar expression pattern across the P(H)AH treatment groups (Figure 3A). Likewise, esr1 and its known regulated/responsive genes [26],[27] such as vg1 and vg3, nots, XBP1 homolog, NUPRI/P8 homolog were among the genes found in a tight cluster that shows consistent expression pattern across the ECs (Figure 3B). The findings indicate that zebrafish whole-adult-organism chemogenomics is able to capture important genes associated with major signaling pathways such as AHR and ER that are deregulated by the compounds. More importantly, the presence of many unknown genes clustering together with the known P(H)AH- or estrogen-responsive genes, suggests that these are potential novel biomarkers for P(H)AHs and ECs, respectively. The consistent expression patterns across several treatment groups containing P(H)AHs or ECs from two different experiments highlights the robustness of these biomarker genes.

Figure 3

Identification of potential biomarker genes.

In addition, we have independently validated 27 and 28 genes that were significantly (P-value <0.05) deregulated in the BAP and DES groups, respectively, using quantitative real-time PCR. We later observed that 16 and 21 of these genes were in the consolidated discriminatory gene sets of P(H)AHs and ECs, respectively (Figure 3; Table S4). As expected, there are a large number of significantly deregulated genes that are excluded in the discriminatory gene sets due to non-fulfillment of the stringent selection criteria but they are nevertheless responsive to BAP and/or DES. The good concordance of the microarray and PCR results for these genes gave us the confidence to pin down these biomarkers at the tissue level.

Biomarker Discovery in Targeted Tissues

To confirm these biomarkers and obtain insights in targeted tissues, we performed a third independent experiment ‘C’ (Table S1) by treating zebrafish with BAP or DES followed by quantitative real-time PCR assay for the previously validated genes (Table S4) in seven tissue types (brain, gills, liver, gut, skin, testis and eyes). Out of the 27 and 28 validated genes for BAP and DES respectively, a total of 81 and 99 positive hits [significant gene deregulation (T-Test P-value <0.05)] were detected among the seven tissues (Table 1 and and2;2; Figure S2). About 37.0% (10/27) and 53.6% (15/28) of the validated genes in BAP and DES, respectively, were deregulated in 4–7 (>50%) of the selected tissues making them excellent biomarker genes for the respective class of compounds (Table 1 and and2).2). Among these 10 and 15 biomarkers for BAP and DES, respectively, 3 are well-known P(H)AH-responsive genes (ahr2, cyp1A1, and TIPARP homolog) and another 3 are well-known EC-responsive genes (esr1, vg1, and vg3), while the remainder 19 are potentially novel as their responsiveness to these compounds are relatively unknown/unreported especially within multiple tissue types. The P(H)AH biomarker genes are mainly associated with xenobiotic metabolism while the EC biomarker genes are mainly associated with molecular transport, metabolism and blood factors. Notably, about 90% (22/25) of these biomarker genes were found in the consolidated discriminatory gene sets. The remaining genes that were deregulated in 1-3 tissues may be useful biomarkers for tissue-specific analysis. Less than 10% of the validated genes in the tissues were found to be inconsistent with the whole fish data and this could be due to having not selected the appropriate responsive tissues or biological variations and cumulative effects of several tissues in the whole fish.

Table 1

Real-Time PCR validated genes for identification of biomarkers in selected tissue samples from fish treated with BAP (250ug/L) for 72 hours.

Table 2

Real-Time PCR of validated genes for identification of biomarkers in selected tissue samples from fish treated with DES (5ug/L) for 72 hours.

Among the selected tissues, eye and skin had the most BAP-responsive genes (63.0% and 51.9% of the validated genes, respectively) followed by gill, liver and testis (40.7% each tissue); these tissues yielded 79.1% of the total positive hits (Table 1; Figure S2). The findings are consistent with mammalian data as these organs (eye, skin, lung, liver and testis) are also known P(H)AH-targeted tissue in mammals [24],[28],[29]. As for tissues with the most DES-responsive genes, liver (92.9%) followed by gut (67.9%) and skin (60.7%) contributed to 62.6% of the total positive hits (Figure S2). The identification of these non-classical estrogen-targeted tissues is consistent with our previous data [30] as well as mammalian data [19],[31]. Activation of AHR or ER signaling pathways, as suggested by up-regulation of known responsive genes, was observed in many of these tissues. Interestingly, cyp1A which was up-regulated by BAP in all 7 tissues, was down-regulated by DES in 5 tissues, suggesting occurrence of similar inhibitory cross-talk between AHR and ER reported in mammalian cells [17],[18],[32]. Taken together, the findings show that the zebrafish shares similar biological responses in terms of molecules, signaling pathways and targeted tissues with mammalian system and is therefore a useful model for inference of chemical biology and health-risk inference in humans.

Biological Function and Human Health-Risk Inferences

To evaluate the potential for extracting biological insights and health-risk inferences in humans, we mapped the consolidated discriminatory gene sets for P(H)AHs and ECs to available corresponding human homologs as previously described [14] and used them for knowledge-based data mining via Ingenuity Pathway Analysis (IPA) software (Figure 4A and and5A).5A). Remarkably, the analysis of the human homologs from the consolidated P(H)AH and EC gene sets listed many affected molecular and cellular functions, perturbed physiological systems and human diseases/disorders that are known to be associated with these compounds [19], [23], [24], [28]–[31] (Figures 4 and and5;5; Tables S5 and S6). These also include canonical signaling pathways such as xenobiotic mechanism signaling (AHR, CYP1A, CYP1B1, CYP2C19, GSTP1, HSP90A, NFE2L2, NOS2A, NQO1, SULT2B1), ERK-MAPK signaling (EFL3, RAC1, STAT1, RPS6KA1) and PPARa/RXRa activation (CYP2C19, HSP90A, LPL, FABP1, NOS2A, ALAS1)for P(H)AH gene set and lipid metabolism (ACSL4, CYP1A1, CYP2C19, CYP51A1), protein ubiquitination pathway (PSMB6, PSMC2, PSMC5/SUG1, PSMC6, PSMD13) and coagulation cascade (FGB, SERPINA1, SERPINC1) for EC gene set. A correlation can be observed between these biological associations suggesting that prolonged or substantial perturbation of these biological functions (Figures 4B and and5B)5B) and physiological systems (Figures 4C and and5C)5C) by a compound would increase the susceptibility/risk of certain diseases/disorders (Figures 4D and and5D).5D). Significantly, cancer was listed among the top most (Fisher Exact Test P-value

1.98×10⁻⁵−4.76×10⁻²) associated disease as most P(H)AHs and some ECs are potent carcinogens. In addition, reproductive system disease, inflammatory disease, hematological disease and neurological disease were significantly associated with both P(H)AHs and ECs as the two classes of compounds are known to affect molecules and functions involved with reproductive system, inflammation, blood and nervous system (Figures 4C–D and 5C–D). Interestingly, the association of psychological disorders with ECs were also significant (P-value

3.34×10⁻⁴−3.72×10⁻²; Table S6). While it is well-known that estrogen affects mental health [33], the grouping of ESR1 with GPX4, PSMC6, DIABLO, FBXO9, XBP1 homologs which have been associated with bipolar disorder [34],[35] suggests that these molecules may also play a role in estrogen-related psychological disorders (Table S6). Several of these patho-physiological systems indicated to be affected by the compounds (Table S5 and S6), such as reproductive, respiratory, dermatological/connective tissue, digestive/metabolic, nervous/neurological and visual/ophthalmic, corroborated with our multiple targeted-tissue analysis which showed that many of the biomarkers were deregulated in these tissues (Tables 1 and and2;2; Figure S2), suggesting that they are good ‘biomarkers of effect’.

Figure 4

Human homolog mapping, biological function and health-risk inferences for P(H)AH consolidated gene set.

Figure 5

Human homolog mapping, biological function and health-risk inferences for EC consolidated gene set.

Thus, a whole-adult-organism representation fits well for human health-risk inferences as the biological information are obtained, not from single tissue but multiple tissues interacting in a complex biological system where diseases/disorders usually develop or off-target effects occur. Incidentally, while zebrafish has always been used to model after human diseases [8], here we demonstrate the potential of zebrafish for predicting disease susceptibility or health risk associated with the exposure to a compound, and this in turn can further help to develop chemical-induced zebrafish models of human diseases.

Novel Biological Insights Inferences

A closer examination of the top connected networks generated by IPA using P(H)AH and EC datasets revealed interesting biological insights (Figure 6). In the P(H)AH network (Figure 6A; P-value

1.00×10⁻³⁹) the well-established xenobiotic-responsive molecules (AHR, CYP1A1, CYP1B1, CYP2C19, HSP90AB1, NFE2L2, NQO1, TIPARP) were linked with key signaling molecules such as NR3C1/GR, STAT1 and IGF1 providing new insights into alternative mechanisms of how P(H)AHs could exert adverse effects resulting in various pathological conditions. Notably in Table S5, 6 categories of P(H)AH-related diseases/disorders were found to be associated with AHR, NR3C1/GR, STAT1 and IGF1, while all 24 categories of the related diseases/disorders were associated with at least one of the four molecules. Interestingly, it was only recently that stronger evidence of cross-talk between AHR and NR3C1/GR are emerging [36],[37]. Our analysis suggests that HSP90 [38] may be one of the mediators, as also observed in the relatively more studied AHR-ER cross-talk [39]. As both NR3C1/GR and STAT1 are known regulators of inflammatory and immune response [40], respectively, the AHR-NR3C1/GR cross-talk may be another pathway that could contribute to the known immuno-toxicity effects of P(H)AHs [41]. The opposing direction of expression for CYP1A1 and CYP2C19, as displayed in both P(H)AHs and ECs networks (Figure 6), are evidence of inhibitory cross-talk between AHR and ER. While CYP1A1 is a known targeted molecule of this inhibitory AHR-ER cross-talk [32], this is the first time a member of the CYP2 family that includes important drug metabolizing enzymes is implicated and its regulation appeared opposite to CYP1A1. Apart from the xenobiotic-responsive molecules, the ECs network (Figure 6B; P-value

1.00×10⁻⁵⁷) linked clusters of molecules associated with proteasome-mediated degradation (PSMC2, PSMC5/SUG1, PSMC6, PSMD13), endoplasmic-reticulum stress response (FKBP11, HM13, RRBP1, PDIA4, SRPRB, SSR1, SSR2, XBP1), and cell cycle/death (CISH, FKBP4, P8, PA2G4, PTCH1), providing insights into cellular homeostasis and pathology associated with ECs and ER [19],[42]. The linking of the 5 known estrogen-responsive transcription factors PSMC5/SUG1, XBP1, P8, PA2G4 and ESR1 suggests that, under the influence of ECs, these transcription factors play important roles in mediating endoplasmic-reticulum stress response, proteasome-mediated degradation and cell cycle/death, thus offering new insights into the regulation of ‘unfolded protein response’ which has been intensely studied due to its association with diseases, drug resistance and its potential as therapeutic targets [43]–[45]. While further investigation is warranted, the analysis demonstrate the discovery potential of zebrafish whole-adult-chemogenomics and serve the purpose of alerting the researcher of the potential molecular interactions and effects induced by the compounds at the early drug discovery stage.

Figure 6

Gene network analysis for biological insights inference.

In summary, we have demonstrated that zebrafish whole-adult-organism chemogenomics is practical and effective for large-scale predictive and discovery chemical biology. Specifically, we have generated robust prediction models, identified and validated biomarker genes in multiple targeted tissues, identified important signaling pathways and biological functions as well as inferred human health risks and biological insights for both P(H)AHs and ECs. As strong and well-represented expression signals are likely to be captured, this approach is valuable for acquiring a molecular snapshot of the chemical-induced biological state of an adult vertebrate, which includes biomarkers of effects, deregulated signaling pathways, as well as possible affected biological functions, perturbed physiological systems and increased health risks. Moreover, this approach allows for rapid sampling in large-scale experiments, abundant sample materials for assays (does not require pooling or amplification of samples) and easy scaling-up of experiments, hence affords greater statistical power for data analysis. These are essentials for successful high-throughput genomics applications and for building up large database for predictive chemical biology [2]. The zebrafish is more cost-effective than the rodent model for in vivo toxicology [12],[13] and there had been proposals and successful attempts of using chemical screens and toxicity testing in whole-adult zebrafish as there are biological processes and diseases that are mainly associated with adults [9],[46]. Therefore, with zebrafish, cost-effective in vivo adult vertebrate chemogenomics can be performed earlier in the drug discovery process and for industrial/environmental toxicology, where high resolution tissue-specific data is yet required but robust and informative in vivo toxicological data is deemed valuable. This will enable researchers to better understand the potential liabilities of new compounds before advancing them to clinical test and may help shift attrition upstream [2],[47],[48] or before allowing contact with the public or releasing them to the environment. This present study has provided a new strategy for genome-wide investigation of chemical-induced biological responses/effects in a whole-adult vertebrate model; where previously such whole-adult-organism chemogenomics approach were thought only feasible in invertebrate models. The realization of its potential can benefit the drug discovery process and toxicology, in particular chemical toxicity testing for environmental health-risk inference.

Materials and Methods

Chemical Exposure and Fish Samples

Three independent experiments (‘A’, ‘B’ and ‘C’) were performed in this study where different batches of adult male zebrafish were exposed to various chemical compounds at different concentrations (Table S1). The selected chemicals represent compounds with toxicological interests and/or environmental-health importance, and the concentrations used were based on available published data or our preliminary acute toxicity exposure experiments conducted for the compounds. Experimental procedures were performed within the guidelines of National University of Singapore's Institutional Animal Care and Use Committee (NUS-IACUC). The fish were immersed in the chemical solutions for 96 hours [experiment ‘A’] or 72 hours [experiment ‘B’ and ‘C’] at a density of 1 fish/200 ml at 27±2°C in a static condition. Control fish were kept in vehicle solution or water under similar condition. Chemical solutions and water were changed daily. At the end of experiment, individual whole fish were snap-frozen and pounded to powder in liquid nitrogen for subsequent total RNA extraction using Trizol reagent (Invitrogen, USA) protocol. For experiment ‘C’, specific tissues were snap-frozen in liquid nitrogen for total RNA extraction using RNeasy Mini Kit (QIAGEN, Germany) followed by DNAseI treatment and heat inactivation. Reference RNA was obtained by pooling total RNA extracted from male and female zebrafish. The integrity of RNA samples was verified by gel electrophoresis, and the concentrations were determined by UV spectrophotometer.

Zebrafish Oligonucleotide Microarray and Hybridization

The oligonucleotide probes for this array were designed by Compugen (USA) and synthesized by Sigma Genesis (USA). The arrays contained 16,416 oligonucleotide probes. The probes were resuspended in 3× SSC at 20 µM concentration and spotted onto in-house poly-L-lysine-coated microscope slides using a custom-built DNA microarrayer in the Genome Institute of Singapore (GIS). The arrays were spotted and quality controlled essentially as described by Eisen and Brown [49].

For fluorescence labeling of cDNAs, 20 µg of total RNA from the reference and sample RNAs were reverse transcribed in the presence of dNTPs mixed with Aminoallyl-dUTP (Sigma, USA) followed by coupling with mono-functional NHS-ester Cy3 and Cy5 dyes (Amersham, USA), respectively. A common reference design is used where equal amount of RNA samples from control and chemical-treated group were labeled with Cy5 and the same amount of common pooled reference RNA is labeled with Cy3. For each array, the Cy5-labeled samples from either the control or chemical-treated group was co-hybridized with the Cy3-labeled common reference. Thus, the respective paired Cy5- and Cy3-labeled cDNAs were pooled, concentrated, and resuspended in DIG EasyHyb (Roche Applied Science) buffer for hybridization at 42°C for 16 h in a hybridization chamber (Gene Machines). After hybridization, the slides were washed in a series of washing solutions (2× SSC with 0.1% SDS, 1× SSC with 0.1% SDS, 0.2× SSC and 0.05× SSC; 30 sec each), dried using low-speed centrifugation, and scanned for fluorescence detection.

Data Acquisition

The arrays were scanned using the GenePix 4000B microarray scanner (Axon Instruments, USA) and the generated images with their fluorescence signal intensities were analyzed using GenePix Pro 4.0 image analysis software (Axon Instruments, USA). All the arrays gave a mean signal to background ratio more than 5 and had >90% of the gene features that gave a measurable signal. Only gene features that were not flagged were extracted and subjected to Lowess normalization for further analyses. The microarray raw data have been formatted to be compliant with MIAME standard.

Statistical Procedures for Microarray Data

Statistical comparison of the relative mean expression level for each gene between test groups [combined representative groups of P(H)AHs or ECs] and their respective control groups from Dataset I or III were performed using Student's T-test and Significance Analysis of Microarray [50] (SAM) yielding respective P- and Q-values for each gene. The resulting P-values were further adjusted for Benjamini and Hochberg False Discovery Rate (FDR). The discriminatory gene sets were selected based on statistical threshold indicated by Q-value, FDR-value and P-value coupled with 1.5 fold-difference between treated versus control samples (Figure S1 and Table S2). The discriminatory gene sets together with their expression data in Dataset I or III were used to train two supervised learning classifiers, kNN and SVM, which generated prediction models for P(H)AHs class and ECs class. The procedure includes a training phase and a testing phase. In the training phase, the discriminatory gene sets together with their expression data from Dataset I or III were used as inputs to produce a set of rules or reference weights as standards/models for the testing phase. A ‘leave-one-out’ cross-validation is usually incorporated to validate the goodness of the model to avoid ‘over-fitting’ it (i.e. only good for predicting the training dataset but not sufficiently generalize to work well on other new and unknown datasets). The testing phase uses the standards created during ‘training’ to assign a discriminator score to each unseen (not used in the training) or unseen sample. Based on this score each sample is placed ‘into’ or ‘out of’ the class and their performances in terms of prediction specificity and sensitivity were determined (Figure 1 and and2;2; Table S3). Fisher's Exact Test was further used to determine that the performances of the prediction models were non-random.

Quantitative Real-Time PCR

Equal amounts of total RNA samples from reference, control and test groups were reverse transcribed to cDNA. The cDNA samples were used for quantitative real-time PCR analysis, performed using the Lightcycler system (Roche Applied Science) with Lightcycler-FastStart DNA Master SYBR Green 1 (Roche Applied Science) according to the manufacturer's instructions. Statistical comparison of the relative mean expression level for each gene between test and control groups was performed using Student's T-test and P-value<0.05 is considered significant.

Human Health-Risk and Biological Insights Inference via Knowledge-Based Data Mining

To evaluate the potential for human health-risk inference, the zebrafish genes were mapped to their corresponding human homologs using the GIS Zebrafish Microarray Annotation Database (http://giscompute.gis.a-star.edu.sg/govind/unigene_db/) http://giscompute.gis.a-star.edu.sg/govind/zebrafish/version2/as previously described [14]. The human homologs of the zebrafish genes from the consolidated gene sets P(H)AHs and ECs were used to mine the human database via Ingenuity Pathways Knowledge Base software (www.ingenuity.com). The ‘Biological and Disease Function Analysis’ was performed to identify biological functions and systems as well as diseases/disorders that were significantly associated with the gene sets. Fisher's Exact test was used to calculate P-values in determining the probability that each biological function, system and disease assigned to that data set is due to chance alone. P-value<0.05 is considered significant by the algorithm. Networks are generated from the available human homologs mapped from the discriminatory gene sets, by maximizing the specific connectivity of the human homologs, which is their interconnectedness with each other relative to all molecules they are connected to in Ingenuity's Knowledge Database. Networks are limited to 35 molecules each to keep them to a functional size and a network score is generated based on the hypergeometric distribution and is calculated with the right-tailed Fisher's Exact Test. In this study, only the top network for the consolidated gene sets P(H)AHs (P-value

1.00×10⁻³⁹) and ECs (P-value

1.00×10⁻⁵⁷) were used for further analysis.

Supporting Information

Figure S1

An overview of the workflow for large-scale predictive and discovery chemical biology using whole-adult zebrafish chemogenomics.

(1.12 MB PDF)

Click here for additional data file.^{(1.0M, pdf)}

Figure S2

Summary of real-time PCR validated genes for identification of biomarkers in 7 selected targeted tissues (brain, gill, liver, gut, skin, testis and eyes).

(0.38 MB PDF)

Click here for additional data file.^{(369K, pdf)}

Table S1

Array samples and concentration of chemical treatments used in the study.

(0.05 MB PDF)

Click here for additional data file.^{(51K, pdf)}

Table S2

Selection criteria and number of genes selected for discriminatory gene sets used for training of prediction models.

(0.01 MB PDF)

Click here for additional data file.^{(12K, pdf)}

Table S3

Summary of the performance of prediction models trained using Dataset I (Figure 1) and Dataset III (Figure 2).

(0.05 MB PDF)

Click here for additional data file.^{(50K, pdf)}

Table S4

Selected BAP- or DES- responsive genes validated to be significant (P<0.05) in whole fish using real-time PCR.

(0.03 MB PDF)

Click here for additional data file.^{(30K, pdf)}

Table S5

Biological associations of the P(H)AH-deregulated human homologs (68) generated by Ingenuity Pathway Analysis (IPA) software.

(0.05 MB PDF)

Click here for additional data file.^{(49K, pdf)}

Table S6

Biological associations of the EC-deregulated human homologs (70) generated by Ingenuity Pathway Analysis (IPA) software.

(0.05 MB PDF)

Click here for additional data file.^{(48K, pdf)}

Acknowledgments

We thank Jing Hu, Wen San Lim and Xiaoming Xu for helping out in some of the chemical treatments and microarray hybridization.

Footnotes

The authors have declared that no competing interests exist.

This work is supported by Biomedical Research Council (BMRC) and Environment and Water Industry (EWI) Development Council, Singapore.

References

1. Bredel M, Jacoby E. Chemogenomics: an emerging strategy for rapid target and drug discovery. Nat Rev Genet. 2004;5:262–275. [PubMed]

2. Ganter B, Tugendreich S, Pearson CI, Ayanoglu E, Baumhueter S, et al. Development of a large-scale chemogenomics database to improve drug candidate selection and to understand mechanisms of chemical toxicity and action. J Biotech. 2005;119:219–244.

3. Custodia N, Won SJ, Novillo A, Wieland M, Li C, et al. Caenorhabditis elegans as an environmental monitor using DNA microarray analysis. Ann NY Acad Sci. 2001;948:32–42. [PubMed]

4. Reichert K, Menzel R. Expression profiling of five different xenobiotics using a Caenorhabditis elegans whole genome microarray. Chemosphere. 2005;61:229–237. [PubMed]

5. Girardot F, Monnier V, Tricoire H. Genome wide analysis of common and specific stress responses in adult Drosophila melanogaster. BMC Genomics. 2004;30:74. [PMC free article] [PubMed]

6. Pedra JH, McIntyre LM, Scharf ME, Pittendrigh BR. Genome-wide transcription profile of field- and laboratory-selected dichlorodiphenyltrichloroethane (DDT)-resistant Drosophila. Proc Natl Acad Sci U S A. 2004;101:7034–7039. [PubMed]

7. Boess F, Kamber M, Romer S, Gasser R, Muller D, et al. Gene expression in two hepatic cell lines, cultured primary hepatocytes, and liver slices compared to the in vivo liver gene expression in rats: possible implications for toxicogenomics use of in vitro systems. Toxicol Sci. 2003;73:386–402. [PubMed]

8. Dooley K, Zon LI. Zebrafish: a model system for the study of human disease. Curr Opin Genet Dev. 2000;10:252–256. [PubMed]

9. Zon LI, Peterson RT. In vivo drug discovery in the zebrafish. Nat Rev Drug Discov. 2005;4:35–44. [PubMed]

10. Pichler FB, Laurenson S, Williams LC, Dodd A, Copp BR, et al. Chemical discovery and global gene expression analysis in zebrafish. Nat Biotechnol. 2003;21:879–883. [PubMed]

11. Cossins AR, Crawford DL. Fish as models for environmental genomics. Nat Rev Genet. 2005;6:324–333. [PubMed]

12. Parng C, Seng WL, Semino C, McGrath P. Zebrafish: a preclinical model for drug screening. Assay Drug Dev Technol. 2002;1:41–48. [PubMed]

13. Hill AJ, Teraoka H, Heideman W, Peterson RE. Zebrafish as a model vertebrate for investigating chemical toxicity. Toxicol Sci. 2005;86:6–19. [PubMed]

14. Lam SH, Wu YL, Vega VB, Miller LD, Spitsbergen J, et al. Conservation of gene expression signatures between zebrafish and human liver tumors and tumor progression. Nat Biotechnol. 2006;24:73–75. [PubMed]

15. Lam SH, Winata CL, Tong Y, Korzh S, Lim WS, et al. Transcriptome kinetics of arsenic-induced adaptive response in zebrafish liver. Physiol Genomics. 2006;27:351–361. [PubMed]

16. Yang L, Yang L, Kemadjou JR, Zinsmeister C, Bauer M, Legradi J, et al. Transcriptional profiling reveals barcode-like toxicogenomic responses in the zebrafish embryo. Genome Biol. 2007;8:R227. [PMC free article] [PubMed]

17. Safe S, Wormke M. Inhibitory aryl hydrocarbon receptor-estrogen receptor alpha cross-talk and mechanisms of action. Chem. Res. Toxicol. 2003;16:807–816. [PubMed]

18. Matthews J, Gustafsson JA. Estrogen receptor and aryl hydrocarbon receptor signaling pathways. Nucl Recept Signal. 2006;4:e016. [PMC free article] [PubMed]

19. Nilsson S, Gustafsson JA. Biological role of estrogen and estrogen receptors. Crit Rev Biochem Mol Biol. 2002;37:1–28. [PubMed]

20. Fried KW, Schneider CM, Schramm KW, Datta A, Chahbane N, et al. From Dioxin to Drug Lead-The Development of 2,3,7,8-Tetrachlorophenothiazine. ChemMedChem. 2007;2:890–897. [PubMed]

21. Safe S, McDougal A. Mechanism of action and development of selective aryl hydrocarbon receptor modulators for treatment of hormone-dependent cancers. Int J Oncol. 2002;20:1123–1128. [PubMed]

22. Dobs A, Darkes MJ. Incidence and management of gynecomastia in men treated for prostate cancer. J Urol. 2005;174:1737–1742. [PubMed]

23. Yang M, Park MS, Lee HS. Endocrine disrupting chemicals: human exposure and health risks. J Environ Sci Health C Environ Carcinog Ecotoxicol Rev. 2006;24:183–224. [PubMed]

24. Luch A. Nature and nurture-lessons from chemical carcinogenesis. Nat. Rev. Cancer. 2005;5:113–125. [PubMed]

25. Xu C, Li CY, Kong AN. Induction of phase I, II and III drug metabolism/transport by xenobiotics. Arch Pharm Res. 2005;28:249–268. [PubMed]

26. Rotchell JM, Ostrander GK. Molecular markers of endocrine disruption in aquatic organisms. J Toxicol Environ Health B Crit Rev. 2003;6:453–496. [PubMed]

27. Wang DY, McKague B, Liss SN, Edwards EA. Gene expression profiles for detecting and distinguishing potential endocrine-disrupting compounds in environmental samples. Environ Sci Technol. 2004;38:6396–6406. [PubMed]

28. Lewis RJ. New York: Van Nostrand Reinhold; 1996. Sax's Dangerous Properties of Industrial Materials. 9th ed. p. 357.

29. Agency for Toxic Substances and Disease Registry (ATSDR) Atlanta: U S Department of Health and Human Services, Public Health Service; 1995. Toxicological profile for Polycyclic Aromatic Hydrocarbons (PAHs). p. 458.

30. Wang H, Tan JT, Emelyanov A, Korzh V, Gong Z. Hepatic and extrahepatic expression of vitellogenin genes in the zebrafish, Danio rerio. Gene. 2005;356:91–100. [PubMed]

31. Zhang W, Andersson S, Cheng G, Simpson ER, Warner M, et al. Update on Estrogen Signaling. FEBS Letters. 2003;546:17–24. [PubMed]

32. Beischlag TV, Perdew GH. ER α-AHR-ARNT protein-protein interactions mediate estradiol-dependent transrepression of dioxin-inducible gene transcription. J Biol Chem. 2005;280:21607–21611. [PubMed]

33. Ostlund H, Keller E, Hurd YL. Estrogen receptor gene expression in relation to neuropsychiatric disorders. Ann NY Acad Sci. 2003;1007:54–63. [PubMed]

34. Kakiuchi C, Iwamoto K, Ishiwata M, Bundo M, Kasahara T, et al. Impaired feedback regulation of XBP1 as a genetic risk factor for bipolar disorder. Nature Genet. 2003;35:171–175. [PubMed]

35. Konradi C, Eaton M, MacDonald ML, Walsh J, Benes FM, et al. Molecular evidence for mitochondrial dysfunction in bipolar disorder. Arch Gen Psychiatry. 2004;61:300–308. [PubMed]

36. Dvorak Z, Vrzal R, Pavek P, Ulrichova J. An evidence for regulatory cross-talk between aryl hydrocarbon receptor and glucocorticoid receptor in HepG2 cells. Physiol Res. In press 2007

37. Sonneveld E, Jonas A, Meijer OC, Brouwer A, van der Burg B. Glucocorticoid-Enhanced expression of dioxin target genes through regulation of the rat aryl hydrocarbon receptor. Toxicol Sci. 2007;99:455–469. [PubMed]

38. Fang L, Ricketson D, Getubig L, Darimont B. Unliganded and hormone-bound glucocorticoid receptors interact with distinct hydrophobic sites in the Hsp90 C-terminal domain. Proc Natl Acad Sci U S A. 2006;103:18487–18492. [PubMed]

39. Caruso JA, Laird DW, Batist G. Role of HSP90 in mediating cross-talk between the estrogen receptor and the Ah receptor signal transduction pathways. Biochem Pharmacol. 1999;58:1395–1403. [PubMed]

40. Aittomäki S, Pesu M, Groner B, Jänne OA, Palvimo JJ, et al. Cooperation among Stat1, glucocorticoid receptor, and PU.1 in transcriptional activation of the high-affinity Fc gamma receptor I in monocytes. J Immunol. 2000;164:5689–1597. [PubMed]

41. Nebert DW, Puga A, Vasiliou V. Role of the Ah receptor and the dioxin-inducible [Ah] gene battery in toxicity, cancer, and signal transduction. Ann NY Acad Sci. 1993;685:624–640. [PubMed]

42. Gomez BP, Riggins RB, Shajahan AN, Klimach U, Wang A, et al. Human X-Box binding protein-1 confers both estrogen independence and antiestrogen resistance in breast cancer cell lines. FASEB J. 2007;21:4013–4027. [PubMed]

43. Ron D, Walter P. Signal integration in the endoplasmic reticulum unfolded protein response. Nat Rev Mol Cell Biol. 2007;8:519–529. [PubMed]

44. Kostova Z, Wolf DH. For whom the bell tolls: protein quality control of the endoplasmic reticulum and the ubiquitin-proteasome connection. EMBO J. 2003;22:2309–2317. [PubMed]

45. Ma Y, Hendershot LM. The role of the unfolded protein response in tumour development: friend or foe? Nat Rev Cancer. 2004;4:966–977. [PubMed]

46. Bayliss PE, Bellavance KL, Whitehead GG, Abrams JM, Aegerter S, et al. Chemical modulation of receptor signaling inhibits regenerative angiogenesis in adult zebrafish. Nat Chem Biol. 2006;2:265–273. [PMC free article] [PubMed]

47. Pritchard JF. Making better drugs: Decision gates in non-clinical drug development. Nature Rev Drug Discov. 2003;2:542–553. [PubMed]

48. Kola I, Landis J. Can the pharmaceutical industry reduce attrition rates? Nat Rev Drug Discov. 2004;3:711–715. [PubMed]

49. Eisen MB, Brown PO. DNA arrays for analysis of gene expression. Methods Enzymol. 1999;303:179–205. [PubMed]

50. Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A. 2001;98:5116–5121. [PubMed]

Articles from PLoS Genetics are provided here courtesy of
Public Library of Science