Significance Cells change their mRNA expression in response to biologically active substances in a dose-dependent manner. Because different genes in a cell show distinct sensitivities to the same substance, changes in the genome-wide mRNA expression profile induced by low and high doses of a substance are essentially different, but this notion has been commonly overlooked in previously published studies. Using a human cell culture model and microarray, we performed genome-wide determinations of gene sensitivities to hormonally active substances with statistically rigorous approaches. Our study provides a conceptual and methodological framework for the systematic examination of gene sensitivities and demonstrates effective detection of nonmonotonic dose-dependent responses, introducing the importance of gene sensitivity analysis to pharmacogenomic and toxicogenomic research.

Abstract Although biological effects of endocrine disrupting chemicals (EDCs) are often observed at unexpectedly low doses with occasional nonmonotonic dose–response characteristics, transcriptome-wide profiles of sensitivities or dose-dependent behaviors of the EDC responsive genes have remained unexplored. Here, we describe expressome analysis for the comprehensive examination of dose-dependent gene responses and its applications to characterize estrogen responsive genes in MCF-7 cells. Transcriptomes of MCF-7 cells exposed to varying concentrations of representative natural and xenobiotic estrogens for 48 h were determined by microarray and used for computational calculation of interpolated approximations of estimated transcriptomes for 300 doses uniformly distributed in log space for each chemical. The entire collection of these estimated transcriptomes, designated as the expressome, has provided unique opportunities to profile chemical-specific distributions of ligand sensitivities for large numbers of estrogen responsive genes, revealing that at low concentrations estrogens generally tended to suppress rather than to activate transcription. Gene ontology analysis demonstrated distinct functional enrichment between high- and low-sensitivity estrogen responsive genes, supporting the notion that a single EDC chemical can cause qualitatively distinct biological responses at different doses. Expressomal heatmap visualization of dose-dependent induction of Bisphenol A inducible genes showed a weak gene activation peak at a very low concentration range (ca. 0.1 nM) in addition to the main, strong gene activation peak at and above 100 nM. Thus, expressome analysis is a powerful approach to understanding the EDC dose-dependent dynamic changes in gene expression at the transcriptomal level, providing important information on the overall profiles of ligand sensitivities and nonmonotonic responses.

The endocrine disrupting chemicals (EDCs) are environmental pollutants that interfere with the endocrine system to disturb biological processes such as development, reproduction, and metabolism (1⇓⇓-4). EDCs consist of a wide variety of manmade and natural compounds with highly diversified chemical structures (4). Some EDCs are agonistic ligands that directly bind to hormone receptors, whereas others inhibit receptor functions as antagonists. For example, xenoestrogens such as Bisphenol A (BPA) (1) and genistein (5) are EDCs that activate estrogen receptors by direct binding. Tributyltin is an obesogenic EDC that binds to peroxisome proliferation-activated receptor γ (PPAR-γ) to activate the retinoid X receptor (RXR)–PPAR heterodimer nuclear receptor complex, enhancing adipogenesis (6). Some other EDCs more indirectly disrupt the endocrine system by affecting hormone synthesis, metabolism, sensitivity, or the negative feedback system (4). For example, thyroid hormone EDCs such as polychlorinated biphenyl congeners reduce the circulating thyroid hormone levels, and some of them may also bind to the thyroid hormone receptors as functional ligands (7).

Evidence is accumulating that the EDCs may cause significant biological effects in humans or animals at doses far lower than the exposure limits set by regulatory agencies (8, 9). In addition to such low-dose effects, an increasing number of studies also support the concept of the nonmonotonic EDC effects, whose dose–response curves show U shapes or inverted-U shapes (8⇓-10). However, whereas nonmonotonic dose–responses are often observed in the endocrine system, this notion conflicts with one of the Bradford–Hill criteria for cause-and-effect relationships that requires stronger effects with greater degrees of exposure (11), leading to controversies when biological effects are observed most strongly with lower rather than higher doses of EDCs (12⇓-14). On the other hand, when an EDC does not show significant biological effects at a high dose, that chemical is generally assumed safe at lower doses, and possible nonmonotonic dose–response characteristics are seldom appreciated in risk assessment. Thus, as the standard toxicological approaches commonly adopted by the industry and regulatory apparatus do not presently assume the low-dose–specific toxicity or nonmonotonic dose–response relationship, it is urgently desired to establish solid scientific frameworks for proper handling of data possibly demonstrating the nonstandard dose-dependent effects. Because most published studies on the EDC effects involve limited numbers of doses for each chemical species, research on mechanisms of the low-dose and nonmonotonic actions is still in its infancy (15, 16). Obviously, the lack of widely applicable standard approaches to study the low-dose and nonmonotonic EDC effects on gene expression is a critical obstacle to understanding the mechanisms of the genetic and genomic actions of EDCs.

Our present study attempts to introduce frameworks for comprehensive analyses and data visualization of the EDC dose-dependent transcriptomal dynamism and its possible nonmonotonic characteristics. We describe the expressome as a library of interpolated approximations of transcriptomal profiles for hundreds of doses uniformly distributed in the log space within the range of doses for which the seed transcriptomes are experimentally determined (ACR, analyzed concentration range). Our expressome analysis of representative xenoestrogens has determined transcriptome-wide profiles of ligand sensitivities of estrogen responsive genes in MCF-7 cells and provided visualization of the nonmonotonic aspects of the transcriptomal effects of BPA inducible genes. We propose that expressome analysis is a powerful approach for the comprehensive and statistically rigorous examination of the dose-dependent dynamic aspects of transcriptomal responses to EDCs.

Discussion An increasing number of studies suggest that the traditional practice of toxicological risk assessments overlook important characteristics of EDCs such as low-dose effects and nonmonotonic dose–response (3, 4, 14). Because animal studies often demonstrate significant biological effects of EDCs at concentrations far lower than the safety limits set by regulatory agencies, appropriateness of the use of “no observable effects limit” for EDC risk assessment has been questioned (8, 25). The nonmonotonic dose–response may result in paradoxical observations that a lower EDC dose can be more toxic than a higher dose of the same chemical (8, 10). Accumulation of transcriptome-wide data on EDC sensitivities of individual EDC responsive genes will be critical for the mechanistic understanding of these extraordinary phenomena. However, there have been no conceptual or methodological frameworks for guiding systematic investigations on this subject. In an effort to develop such a framework, in the present study we propose the expressome analysis along with its two practical applications—namely, the SDC analysis of sensitivities of EDC responsive genes (Fig. 3, Fig. S4, and Table S2) and the expressomal heatmap visualization of nonmonotonic dose–response (Fig. 4). These unique tools will provide not only the research communities but also the governmental regulatory apparatus with opportunities to perform comprehensive, objective, and precise dose-dependent effects of toxic substances on mammalian genomes. As the costs and the accuracy of transcriptomal analyses are rapidly improving due to the technological revolutions such as deep sequencing, it is practical for toxicological risk assessments to take the anomalous dose–response relationships and changes in sensitivities to natural or other exogenous substances into consideration. As our present study demonstrates, the nonmonotonicity of the dose–effect relationship possibly characterizing certain toxic substances may not be clearly evident when each series of dose–response data are independently examined, due to the relatively stochastic nature of this phenomenon. To detect the nonmonotonicity, data obtained from multiple repeated experiments should be processed without subjective removal of apparently “outlier” data and summarized using an appropriate statistical approach. The expressomal analysis provides unique opportunities to perform these tasks with proper objectivity and sensitivity. To perform expressomal analyses, a sufficient number of dose-dependent transcriptomal changes (seed transcriptomes) have to be experimentally determined. Credibility and resolution power of expressomal analyses are primarily determined by the number and distribution of seed transcriptomes, which should be examined in more detail by future studies. In the present study, the MCF-7 exposure time to estrogens was fixed at 48 h to focus on the dose-dependent aspects of gene expression. Responses of the so-called “early responder genes” (e.g., MYC) were no longer detected after 48 h of exposure to estrogens (21, 22). Future studies should attempt to examine time-dependent changes in the expressomes, which may be visualized in the 3D space whose x, y, and z axes are dose, mRNA expression, and time, respectively. The expressomal SDC analysis introduced in this study provides unique opportunities to obtain the following information: (i) transcriptome-wide profiles of sensitivities of EDC responsive genes (Fig. 3, Fig. S4, and Table S2), (ii) statistical significance of differential sensitivities between two or more groups of genes (such as EDC inducible genes and suppressible genes; Table S2), (iii) identification of genes with distinct EDC sensitivities (Fig. 3, Fig. S4, and Dataset S1), and (iv) demonstration of functional heterogeneities among genes with differential EDC sensitivities (Dataset S1). Because genes responsive to a single EDC have significantly varying sensitivities, distinct subsets of EDC responsive genes are induced or suppressed by different concentrations of an EDC. Consequently, genomic responses to a low-dose EDC exposure are presumably distinct from responses to a higher dose, and this may cause low-dose–specific biological effects of EDCs. Future studies should examine presumable biological consequences of the differentially expressed EDC responsive genes at different EDC doses. Molecular mechanisms determining sensitivities of estrogen responsive genes can be very diverse, as depicted in Fig. S6. Transcriptional activation or suppression by estrogen receptors (ERs) involves differential recruitments of specific coregulators on different target genes (26). Specific DNA base sequences at and around the ER binding sites serve as allosteric modulators of ER conformation and thereby affect selection of coregulator recruitment and transcriptional activities (27). Multiple other transcription factors and chromatin-modifying mechanisms also affect the genomic actions of ERs, which may bind to chromatin hundreds of kilobases away from the transcription initiation sites (28). As a type I nuclear receptor, ERs do not interact with DNA until ligand binding. Therefore, it is possible that different ligands, including endogenous ligands such as 27-hydroxycholesterol (29), may alter the affinity of ERs to various DNA targets and coregulators, resulting in a highly complex, and possibly nonmonotonic, dose–response relationship. The expressomal SDC analysis may provide important clues to recognize characteristic sets of coregulators and other molecular events differentially involved in the high-sensitivity and low-sensitivity transcriptional responses to estrogens. On the other hand, Kortenkamp et al. reported that the cell culture effects of estrogenic chemical mixtures are well predictable from the strength of estrogenicity of each component by the concentration addition model (30), implying a single or relatively limited number of mechanisms involved in the estrogenic actions of xenoestrogens. However, they also observed the “effect modulation” phenomenon, in which inclusion of low estrogenicity chemicals results in overall reduction in the estrogenicity strength of the mixture (30). The expressomal SDC analysis may be able to determine detailed transcriptomal effects of EDC mixtures, providing clues to understand mechanisms of the concentration additive features and the effect modulation phenomenon. Mechanisms involved in nonmonotonic dose–response of EDC effects are largely unknown. In the whole animal, the negative feedback mechanisms of the endocrine system and other factors affecting metabolism or distributions of hormones may contribute to the nonmonotonicity. Our present study provides evidence that nonmonotonic EDC dose–response can be observed even in cell culture models (Fig. 4 C and G and Fig. S5H). Because MCF-7 cells express only minimal amounts of ERβ (31), the observed nonmonotonicity was unlikely due to differential transcriptional regulations by the two isoforms of estrogen receptors. Because the low-dose gene expression peaks (orange rectangles in Fig. 4 C and G) of BPA inducible genes were observed only for genes showing the highest BPA sensitivities in the major mRNA expression peaks, mechanisms of the nonmonotonicity may be related with the high BPA sensitivity. Future studies are needed to determine transcriptional coregulators, epigenetic marks, and other transcription factors involved in the high-sensitivity BPA activation of transcription and their possible roles in the nonmonotonic responses. In summary, the present study introduced the expressome, a collection of transcriptome-wide approximations of continuous dose-dependent mRNA expression curves generated based on limited numbers of experimentally determined transcriptomes. Analyses of the expressome provide unique information about genome-wide profiles of gene sensitivities to transcriptionally active substances and nonmonotonic aspects of dose-dependent gene expression. Expressomal analyses may contribute to the development of computational handling of genomics data relevant to gene sensitivities of EDCs for pharmaco- and toxico-genomics research and risk assessment.

Methods Cell Culture and Determination of Transcriptomes by Affymetrix Microarray. MCF-7 cell culture, 10-d cell proliferation assay for estrogenicity, determination of transcriptomal profiles using Affymetrix Human Genome U-133 plus 2.0 microarray, and hierarchical clustering and heatmap visualization were described in our preceding studies (21, 22). Outlier transcriptomal data were identified by box plot and Q–Q plot analyses and eliminated from the study. All experimentally determined transcriptomal data are available from the National Center for Biotechnology Information Gene Expression Omnibus database (accession no. GSE50705) and our website (http://mplwebserver.partners.org). TaqMan real-time qPCR determination of BIK and WISP2 mRNA expression was also performed as previously described (22). All estrogenic agents were highest grade chemicals purchased from Sigma-Aldrich. Computational Generation and Analyses of the Expressome Models. Computational tasks were performed using R/CRAN package mgcv 1.6–1 (23) and homemade scripts written in Python 2.7 or Ruby 1.9 programming languages. Details are described in SI Text.

Acknowledgments This study was supported by Susan G. Komen for Cure Grants FAS0703860 and KG090515 (to T.S.) and NIH Grants U41-HG004059 and R01-HL086601 (to V.J.C.).

Footnotes Author contributions: T.S. and K.J.I. designed research; T.S., N.F.R., and K.R.C. performed research; T.S., M.S., and V.J.C. contributed new reagents/analytic tools; T.S., M.P., M.M., and V.J.C. analyzed data; and T.S. and K.J.I. wrote the paper.

The authors declare no conflict of interest.

Data deposition: The data reported in this paper have been deposited in the Gene Expression Omnibus (GEO) database, www.ncbi.nlm.nih.gov/geo (accession no. GSE50705).

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1315929110/-/DCSupplemental.