Identification of potential endocrine disrupting chemicals using gene expression biomarkers

Accepted Manuscript Identification of potential endocrine disrupting chemicals using gene expression biomarkers

J. Christopher Corton, Nicole C. Kleinstreuer, Richard S. Judson PII: DOI: Article Number: Reference:

S0041-008X(19)30291-1 https://doi.org/10.1016/j.taap.2019.114683 114683 YTAAP 114683

To appear in:

Toxicology and Applied Pharmacology

Received date: Revised date: Accepted date:

27 February 2019 5 July 2019 15 July 2019

Please cite this article as: J.C. Corton, N.C. Kleinstreuer and R.S. Judson, Identification of potential endocrine disrupting chemicals using gene expression biomarkers, Toxicology and Applied Pharmacology, https://doi.org/10.1016/j.taap.2019.114683

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Identification of potential endocrine disrupting chemicals using gene expression biomarkers

T

J. Christopher Corton1,* Nicole C. Kleinstreuer2 , and Richard S. Judson3

1

AN

US

CR

IP

Integrated Systems Toxicology Division, National Health and Environmental Effects Research Laboratory, US Environmental Protection Agency, Research Triangle Park, NC 27711 2 National Toxicology Program Interagency Center for the Evaluation of Alternative Toxicological Methods, National Institute of Environmental Health Sciences, Research Triangle Park, NC 3 National Center for Computational Toxicology, US Environmental Protection Agency, Research Triangle Park, NC * Corresponding author

PT

ED

M

Chris Corton Integrated Systems Toxicology Division National Health and Environmental Effects Research Lab US Environmental Protection Agency 109 T.W. Alexander Dr. MD-B105-03 Research Triangle Park, NC 27711

AC

CE

[email protected] 919-541-0092 (office) 919-541-0694 (fax)

Key words: toxicogenomics; biomarkers; endocrine disrupting chemicals; nuclear receptors; estrogen receptor; androgen receptor Email addresses of authors: Chris Corton: [email protected] Nicole Kleinstreuer: [email protected] Richard Judson: [email protected]

Page 1

ACCEPTED MANUSCRIPT Abstract Recent technological advances have moved the field of toxicogenomics from reliance on microarray platforms to high-throughput transcriptomic (HTTr) technologies that measure global gene expression. Gene expression biomarkers are emerging as useful tools for interpreting gene

IP

T

expression profiles to identify perturbations of targets of xenobiotic chemicals including those that act as endocrine disrupting chemicals (EDCs). Gene expression biomarkers are lists of

CR

similarly-regulated genes identified in global gene expression comparisons of cells or tissues 1)

US

exposed to known agonists or antagonists of the transcription factor (TF) and 2) after expression of the TF itself is knocked down/knocked out or overexpressed. Estrogen receptor α (ERα) and

AN

androgen receptor (AR) biomarkers have been shown to be very accurate at identifying both

M

agonists (94-97%) and antagonists (93-98%) in microarray data derived from human breast or prostate cancer cell lines. Importantly, the biomarkers have been shown to accurately replicate

ED

the results of computational models that predict ERα or AR modulation using multiple ToxCast

PT

HT screening assays. An integrated screening strategy using sets of biomarkers that simultaneously predict various EDC targets in relevant cell lines should simplify chemical

CE

screening without sacrificing accuracy. The biomarker predictions can be put into the context of

AC

the adverse outcome pathway framework to help prioritize chemicals with the greatest risk of potential adverse outcomes in the endocrine systems of animals and people.

Introduction Exposure to endocrine disrupting chemicals (EDCs) is a risk factor for disruption of the development of organ systems in humans and wildlife as well as for cancer and metabolic

Page 2

ACCEPTED MANUSCRIPT syndromes [1, 2]. In the 1990s the recognition that man-made chemicals can interfere with endocrine functions led to legislation in the United States which resulted in a mandate that the United States Environmental Protection Agency (EPA) develop a screening program for potential EDCs. In the resulting Endocrine Disruptor Screening Program (EDSP), approximately

T

10,000 existing chemicals would be evaluated for their potential to disrupt the estrogen,

IP

androgen, and thyroid signaling systems. Under the program guidelines, a battery of in vitro and

CR

short-term in vivo screening assays were developed for chemical hazard screening, to be followed by longer term, more definitive in vivo Tier 2 tests for endocrine disrupting activity

US

(http://www.epa.gov/endocrine-disruption/endocrine-disruptor-screening-and-testing-advisory-

AN

committee-edstac- final).

M

The battery of assays included those that assess the activity of a number of nuclear receptor family members that when inappropriately activated or inhibited can interfere with

ED

normal endocrine signaling. These receptors including two estrogen receptors (ERα, ERβ), the

PT

androgen receptor (AR), and two thyroid hormone receptors (TRα, TRβ) act as ligand-binding transcription factors (TFs) that can be activated or suppressed by xenobiotic chemicals, resulting

CE

in altered gene expression in susceptible tissues. The EPA’s vision for the EDSP in the 21st

AC

Century (EDSP21) to prioritize chemicals includes eventually replacing or providing alternatives to some or all of the current EDSP Tier 1 screening assays, including via the use of in vitro highthroughput screening (HTS) assays coupled with computational modeling (https://www.federalregister.gov/articles/2015/06/19/2015-15182/use-of-high-throughputassays-and-computational-tools-endocrine-disruptor-screening-program-notice; last accessed August 10, 2018). The EPA ToxCast screening program (http://www.epa.gov/chemicalresearch/toxicity-forecasting) and the cross-agency Tox21 program Page 3

ACCEPTED MANUSCRIPT (http://www.ncats.nih.gov/tox21) currently carry out HTS assays in chemical safety evaluation programs. These programs have screened more than 1800 chemicals in as many as 700 HTS assays [3] and the assays and resultant data streams have been invaluable in prioritizing chemicals for further testing and for identification of potential EDCs.

IP

T

One drawback of the HTS program is that the full suite of commercially-available, off-

CR

the-shelf ToxCast assays examine only a small fraction of the potential chemical targets in the genome (~350 gene products out of ~20,000) [3], potentially missing effects on important ED-

US

related pathways. Due in part to this data gap, there is increased interest in utilizing highthroughput transcriptomic (HTTr) approaches that can assess the expression of the entire genome

AN

in HT formats as part of chemical screening programs. In this vision, HTTr approaches would be

M

used as “Tier 0” screens followed by selected Tier 1 assays to confirm predictions derived from the gene expression analysis. Opportunities for performing HTTr screening have come from

ED

recent advances in HTS approaches to measure gene expression in chemically-treated cells. Past

PT

chemical assessment studies almost exclusively utilized commercially-available microarrays (e.g., Agilent and Affymetrix microarrays) which assess the expression of the entire genome but

CE

have low-throughput. Investigators have recently turned to a number of sequencing-based

AC

approaches. Whole transcriptome RNA-Seq allows a deeper interrogation of complex changes in gene expression patterns including assessment of non-mRNA species (microRNAs and long noncoding RNAs) as well as alternative splicing, but most of the standard protocols are labor intensive and cost prohibitive for high-throughput use [4]. Higher-throughput alternatives to RNA-Seq have been recently described including PLATE-Seq which allows samples in 96-well plates to be simultaneously profiled at a lower cost per sample [5] and Digital RNA with pertUrbation of Genes (DRUG-Seq) which is a massively parallelized, automated, method which Page 4

ACCEPTED MANUSCRIPT drives down the price per sample to as low as $2-4 [6]. Targeted RNA-Seq is an alternative strategy to assess gene expression. The TempO-Seq platform (BioSpyder, Inc.) [7] has been adapted to HTS to allow measurement of gene expression from lysates of treated cells and is currently used by the EPA [8]. These HTTr techniques promise to generate ever increasing

T

amounts of data and predictions of molecular effects that have the potential to replace more

IP

conventional HTS assays and transform toxicity testing.

CR

This review describes an approach for using HTTr data to identify chemical- induced

US

changes in the activity of important targets of EDCs that is based on our published studies. The approach uses gene expression biomarkers built from genes regulated by the factor to be

AN

measured. Due to their accuracy and relative ease of use and discussed in detail below, these

M

approaches have the potential of complementing or replacing current batteries of HTS tests and will allow an integrated assessment of most of the major factors that represent targets for EDCs.

ED

While we focus much of the review on using the previously published gene expression

PT

biomarkers for ER and AR [9, 10], we discuss how the strategy can be extended to include measurement of essentially all important targets of environmental chemicals. We do not cover

CE

systems biology approaches to identifying groups of genes that define similarity in chemical-

AC

induced effects (e.g.,[11, 12] ) as these approaches are not TF focused. Challenges in moving toward more routine HTTr testing will be discussed.

Building gene expression biomarkers for transcription factors linked to ED A large number of computational methods have been used in past studies to identify groups of differentially-expressed genes (sometimes called biomarkers, signatures, or

Page 5

ACCEPTED MANUSCRIPT classifiers) that allow the prediction of a predefined set of changes in molecular, cellular, and tissue endpoints often linked to health or disease. The following is a definition of a biomarker from the Biomarkers Definitions Working Group relevant to the use of biomarkers in toxicity testing; “a characteristic that is objectively measured and evaluated as an indicator of normal

T

biological processes, pathogenic processes, or pharmacologic responses to a therapeutic

IP

intervention” [13]. Biomarkers including those constructed from gene expression profiles should

CR

ideally have three properties. First, the biomarker must be highly accurate in predicting the endpoint it was designed for in the context of its use, which could be either a rodent tissue or a

US

human cell line used in screening. Second, the biomarker must make predictions using a minimal

AN

number of samples from the experimental system tested and a relatively small number of genes (i.e., tens rather than thousands), ensuring feasibility and potentially reducing cost. Third, the

M

biomarker must be insensitive to the way in which investigators assess gene expression,

ED

including differences in sample processing or assay technologies to ensure its generalizability across labs. As discussed below, the previously described biomarkers for ER and AR exhibit

PT

these properties.

CE

Our lab has recently employed a relatively straightforward weight of evidence (WOE)

AC

approach to building gene expression biomarkers that relies on the known biology of the factor to be predicted. Past studies have shown that the use of biological information adds a layer of validation and prioritization that can be exploited for biomarker discovery [14]. Two biomarker panels have been built and tested in the lab that measure human TFs (ERα and AR) important in ED. These biomarkers and the only other gene expression biomarkers used to predict changes in transcription factor activity are described in Table 1. The general approach can also be used to build biomarkers that predict pathways and phenotypic effects relevant to

Page 6

ACCEPTED MANUSCRIPT cancer [15]. The biomarkers are built starting with statistically-filtered lists of genes derived from microarray comparisons of chemical or genetic conditions known to activate or inhibit the TF of interest (Figure 1). The published biomarkers have been built using a series of filtering steps. In the first filter, genes are identified that exhibit consistent changes in expression (either

T

increased or decreased) after exposure to TF agonists. This step often yields hundreds of genes

IP

and utilizes comparisons from related cell lines in which the TF is active (e.g., human ERα-

CR

positive breast cancer cell lines). The second filter identifies genes that exhibit consistent but opposite expression after exposure to antagonists [9, 10]. Identification of these genes increases

US

the potential that the biomarker can be used as a tool to find both agonists and antagonists in the

AN

same assay. The third filter, implemented for most of the biomarkers, employs a genetic filter to select genes that are TF-dependent. For example, the AR biomarker genes were filtered for those

M

with consistent expression in prostate cancer cells expressing constitutively active mutants of AR

ED

or after knockdown of the AR gene using gene-specific siRNAs. Biomarkers characterized in the lab that are useful for understanding chemical effects in the rodent liver have employed a similar

PT

genetic strategy using comparisons between chemical-treated wild-type and TF-null mouse livers

CE

[16-22]. In the last step, genes are filtered for those with relatively robust expression changes (usually an absolute average fold-change across the agonist exposures of ≥ 1.5-fold). Not

AC

surprisingly, the resulting biomarker genes exhibit a pattern of expression remarkably similar to the agonists used in the first filtering step. The constitutively active mutants of AR resulted in a profile similar to that of AR agonists, while knockdown of AR expression resulted in a profile similar to that of AR antagonists [10]. For each biomarker, there is a consistent hierarchy of expression changes across the agonist treatments from highly induced/inhibited to genes with more modest changes. The resulting biomarker genes are considered TF-dependent but chemical-

Page 7

ACCEPTED MANUSCRIPT independent, i.e. they are regulated regardless of how the structurally-diverse chemicals interact with the TF. This is an important point; the biomarker is most useful if the genes are not regulated by only a subset of chemicals that work through a narrow set of chemical-TF interactions (described in greater detail below).

IP

T

Biomarker gene lists are evaluated a number of ways to confirm that identified genes are regulated by the TF. Typically, the genes are examined for enrichment with respect to predicted

CR

pathways. For example, the top pathways identified by Ingenuity Pathways Analysis enriched

US

with ERα biomarker genes included the pathways “Ovarian Cancer Signaling,” “Endometrial Cancer Signaling,” and “Estrogen-Dependent Breast Cancer Signaling” [9]. An analysis of

AN

upstream activators predicted to regulate the genes (by Ingenuity Pathway Analysis) has also

M

been used to confirm that the gene set is regulated by the TF. A reexamination of archived chromatin immunoprecipitation coupled with DNA sequencing (ChIP-Seq) datasets is used to

ED

determine if the biomarker genes are physically associated with the TF and thus likely to be

PT

directly regulated. The ChIP-Seq approach has consistently shown that the TF directly regulates most of the identified biomarker genes [9, 10, 21]. Consistent with the biomarker requirement of

CE

relatively small sets of genes to be used for prediction discussed above, the resultant gene lists

AC

are usually in the range of 50-150 genes. Once the biomarker has been built and characterized, the genes can be used in the analysis of global gene expression data by comparing the biomarker profile to lists of statistically- filtered genes after chemical treatment. These tests ask whether the chemicalinduced profile exhibits similarities in direction and magnitude of gene expression changes to the biomarker. To perform these comparisons computationally, we have employed the rank-based nonparametric Running Fisher algorithm [23] which in a pair-wise manner determines the Page 8

ACCEPTED MANUSCRIPT significance of the correlation of the overlapping genes and determines if the correlation is positive or negative. Thus, if the correlation is positive, the TF is activated, and if the correlation is negative, the TF is suppressed. For chemicals, the implication is that activation is caused by an agonist and suppression is caused by an antagonist; however, further testing is required to

T

distinguish direct binding from other mechanisms (discussed below). In the Running Fisher test,

IP

the genes are first ranked by absolute fold-change which removes dependence on the range of

CR

absolute values of fold-change values between filtered gene lists derived from different platforms, minimizes some of the effects of the normalization methods used, and facilitates

US

comparisons between different studies, platforms, and analysis methods [23]. The Running

AN

Fisher test is conceptually similar to the Gene Set Enrichment Analysis (GSEA) method [24, 25], but differs from GSEA in the assessment of the statistical significance, as p-values are computed

ED

M

by a Fisher exact test rather than by permutations [23].

PT

Determining predictive accuracy of the biomarkers

CE

To determine predictive accuracy of the biomarkers as well as to perform HTS, we have capitalized on the availability of thousands of microarray comparisons that are found in a large

AC

commercial genomic database (Illumina’s BaseSpace Correlation Engine). The database is built of lists of statistically- filtered genes derived from archived raw microarray files mostly from the genomic databases Gene Expression Omnibus or ArrayExpress [23]. The database consists of ~137,000 gene lists derived from ~22,400 studies and comparisons from treatment of ~4680 chemicals (as of November, 2018). The database also includes results derived from genetic perturbations including wild-type vs. nullizygous mouse comparisons, human cells treated with scrambled control siRNAs vs. siRNAs against specific gene targets, and comparisons between Page 9

ACCEPTED MANUSCRIPT the wild-type gene and the corresponding gene with mutations that have diverse effects on the encoded protein. The gene lists are derived from microarray comparisons using standard accepted methods for microarray analysis after invoking a number of quality control filters which remove approximately one-third of all studies [23]. Species-specific databases for mouse,

T

rat, and human have been compiled for biomarker investigations, which includes relevant

IP

information about each microarray comparison, facilitating downstream interpretation of

CR

biomarker findings. The comparisons are annotated with the broad category (e.g., chemical, gene, etc.) and the specific factor being examined (e.g, estradiol, AR). Additional information

US

includes the GEO or ArrayExpress study number, tissue/cell line used in the experiment, and

AN

information about the chemical treatment (concentration/dose, time of exposure).

M

To determine predictive accuracy and screening, compendia of tissue- and speciesspecific microarray comparisons have been assembled that are used in conjunction with specific

ED

biomarkers (Figure 2, top). For example, the ERα biomarker was tested and used for screening

PT

in a human breast cancer cell line compendium consisting mostly of comparisons generated in the ERα-positive MCF-7 cell line [9]. This compendium contains over 1400 biosets from

CE

chemically-treated MCF-7 cells, most of which came from the CMAP 2.0 drug study [25]. There

AC

are also hundreds of comparisons of hormone effects mostly of 17β-estradiol (E2) treatments as well as experiments which examined the effects of overexpression or knocking down expression of approximately 200 genes including the ESR1 gene itself. The AR biomarker was characterized and used in the context of a prostate cancer cell line compendium derived from human prostate cancer cell lines, including over 1600 biosets from chemically-treated cells [10], most of which also came from the CMAP 2.0 study using the AR-negative PC3 cell line [25]. Together these compendia of microarray profiles have been invaluable to identify comparisons that can be used

Page 10

ACCEPTED MANUSCRIPT to build biomarkers, to perform predictive accuracy tests, and to identify novel chemicals and other experimental conditions that modulate the TF of interest (Figure 2, bottom).

T

Use of gene expression biomarkers to identify potential EDCs in vitro: accurate replication

IP

of ToxCast model predictions

CR

While the initial focus of the lab was to use biomarkers to make predictions of in vivo effects [16-22, 26, 27], we have increasingly focused our efforts on predicting effects in

US

microarray profiles derived from in vitro cultures of human cells. This effort was due in large

AN

part to the increased emphasis in the Office of Research and Development at the EPA on in vitro screening and prioritizing chemicals and on using dose-response data with exposure modeling to

M

perform high-throughput risk assessments [28]. Using the approach discussed above, we are

ED

presently building predictive biomarkers of high value targets of xenobiotic chemicals that include those involved in ED as well as xenobiotic metabolism, cancer signaling, and stress

PT

response [15]. The goal is to build a large battery of predictive biomarkers that will allow

CE

accurate assessment of key events associated with a wide range of different toxicities. However, to make that a reality, it was necessary to carry out proof-of-principle experiments to

AC

demonstrate that the approach can be applied within the context of HTTr screening data.

Using the ER biomarker to replicate the results of the ER computational model Our first effort to create human biomarkers focused on prediction of ERα activity [9]. The 46 ERα biomarker genes were evaluated as a predictive tool by comparison to a number of annotated gene expression datasets from experiments using MCF-7 cells, including those Page 11

ACCEPTED MANUSCRIPT evaluating the transcriptional effects of hormones and chemicals. The biomarker was evaluated for the ability to predict ERα modulation by chemicals with known activity. The ERα reference chemicals were taken from the OECD TG457 BG1 guidance document ([29]; https://doi.org/10.1787/9789264185395-en). Out of the 45 reference chemicals listed, there were

T

21 chemicals represented by 98 microarray comparisons in the MCF-7 compendium, including

IP

13 agonists, 3 antagonists, and 5 inactives. The biomarker correctly classified 18 out of 21 (86%)

CR

of these reference chemicals including those that are classified as “weak” and "very weak" agonists. Using a larger set of 141 comparisons from chemical- and hormone-treated cells, the

US

biomarker gave balanced accuracies for prediction of ERα activation or suppression of 94% and

AN

93%, respectively.

The biomarker predictions were compared to predictions from an ERα computational

M

model based on 18 in vitro HTS assays measuring different key events in the ERα pathway (e.g.

ED

binding, transcription, proliferation) carried out as part of the EDSP HTS program [30]. The

PT

biomarker was able to correctly classify 105 of the 114 chemicals in common. Accuracy tests for ERα activation or suppression gave balanced accuracies of 95% or 98%, respectively. The

CE

chemical chrysin was the only false negative for activation or suppression of the chemicals tested

AC

and was found to be a very weak modulator in the Judson et al. study with a score very near the cutoff used [30]. Of the 4 “false positive” chemicals for ERα activation, there were two chemicals with evidence of estrogenicity; 4-nonylphenol was shown to activate ERα in MCF-7 cells [31], and theobromine was identified as an estrogenic chemical in a study of drug repurposing [32]. Of the four false positive chemicals for ERα suppression, digoxin has been linked to estrogenicity in other cellular contexts [33], indicating that digoxin can act as a selective ERα modulator. These results from our previous work demonstrated that the ERα

Page 12

ACCEPTED MANUSCRIPT biomarker can accurately identify ERα modulators in large collections of microarray data derived from MCF-7 cells.

T

Using the AR biomarker to replicate the results of the AR computational model

IP

A similar analysis of the predictive capabilities of the AR biomarker was carried out in

CR

our previous publication [10]. The 51-gene biomarker was evaluated as a predictive tool using annotated sets of chemicals and hormone treatments. Because very few chemicals were

US

evaluated in the compendium and environmental chemicals were found to more likely be AR

AN

antagonists than AR agonists [34], we generated a microarray dataset in the prostate cancer cell line LAPC-4 of 28 AR reference antagonists at ~10-20 uM that were also examined in ToxCast

M

assays. LAPC-4 cells encode an AR gene that is not mutated [35] as in LNCaP cells [36] nor

ED

amplified as in VCaP cells [37]. AR biomarker predictions of the chemicals were compared to OECD AR reference chemical classifications [34, 38, 39]. The biomarker correctly classified 16

PT

out of the 17 AR reference antagonists including those that are “weak” and “very weak”.

CE

Overall, 18 out of 21 OECD classifications (86%) were correctly classified by the biomarker with the 2 false negatives marginally inactive. In an examination of a larger set of 158 chemical

AC

and hormone treatments including 95 chemicals, the biomarker gave balanced accuracies for prediction of AR activation or AR suppression of 97% or 98%, respectively [10]. Predictions based on microarray profiles from the LAPC-4 cells treated with the 28 antagonists described above were compared with the predictions from an AR computational model based on 11 in vitro HTS assays which examined different key events in AR pathway activity [34]. Accuracy tests for AR suppression gave a balanced accuracy of 93%. There was 1

Page 13

ACCEPTED MANUSCRIPT false positive and 1 false negative. The false negative cyproterone acetate approached significance using the biomarker approach, indicating that if a full dose response similar to the ToxCast protocol was used (up to 100 uM), cyproterone acetate would have reached significance. The false positive was the 5α-reductase inhibitor finasteride, an OECD reference

T

antagonist for the rodent Hershberger assay, and identified as an antagonist in the two Tox21 AR

IP

antagonist assays [40]. Although assessment of the predictive accuracy of a larger set of

CR

chemicals is needed, the results indicate that under some screening conditions, blocking 5αreductase activity may indirectly inhibit AR by decreasing the levels of the more potent AR

US

activator dihydrotestosterone (DHT) compared with testosterone. If finasteride was considered a

AN

true positive for antagonism, then the balanced accuracy increases to 95% [10].

M

The high concordance between the biomarker classifications and the HTS methods was remarkable considering the deficiencies inherent in some of the microarray comparisons used.

ED

For the MCF-7 comparisons used in the ERα analysis, the biosets from the CMAP 2.0 dataset

PT

included chemical comparisons in which statistically-significant gene lists were derived using a t-test comparison between multiple control samples vs. 1 treated sample examined at only one

CE

time point (6 hours) and concentration level (approximately 10 µM for most chemicals). As the

AC

concentration of a chemical and time of exposure are critical factors determining effects, evaluation of a range of concentrations and time points is necessary to reduce the risk of false negatives in toxicity testing. It is possible that if each chemical evaluated by microarrays was examined under multiple conditions, the accuracy of the biomarker method would approach 100%. Despite these inherent weaknesses in the dataset used, our approach did not misclassify many true positives. For present and future HTTr screens, the ability to identify chemicals of concern should be greatly improved over the present analysis when multiple biological

Page 14

ACCEPTED MANUSCRIPT replicates, concentrations, and times of exposure are examined. It should be noted that using our methods, the biomarkers are accurate independent of the platform used to carry out the gene expression analysis which included a variety of microarrays (mostly Affymetrix, Agilent, and Illumina) and RNA-Seq. We recently reported that the accuracy of one biomarker (the TGx-DDI

T

biomarker; Table 1) is similar across platforms including microarrays, RT-PCR, and direct

IP

digital nCounter data [15], satisfying one of the requirements for a biomarker. Based on our

CR

previous studies, our approach using biomarkers and the Running Fisher test accurately identified chemicals in gene expression data that modulate AR and ERα, two of the major targets

US

of EDCs. These studies satisfy the accuracy requirements of a biomarker and strongly support

AN

the use of our method as an alternative to the multiple HTS assays run in ToxCast to identify AR

ED

M

and ERα active chemicals.

Biomarkers can identify modulation of AR and ERα by multiple mechanisms

PT

The well annotated microarray compendia we used in our studies allowed us to determine

CE

the ability of the biomarkers to detect different chemical- induced mechanisms of AR and ERα modulation. Like other nuclear receptor family members, AR and ERα regulate target gene

AC

expression through multiple mechanisms, some of which are not completely characterized. The best known classical ligand-dependent pathway involves ligand binding to the receptor by agonists, receptor interaction with co-activator proteins, and binding of the transcriptional complexes to estrogen or androgen response elements (ERE or ARE) in the DNA [41, 42]. For either AR or ERα the transcriptional complexes regulate an overlapping set of genes determined by ligand- induced receptor conformations, the population of bound co-regulator proteins, and position of the receptors on the DNA. Antagonists inhibit receptor activity by promoting Page 15

ACCEPTED MANUSCRIPT interactions with co-repressors that block productive interactions with the transcriptional machinery. Chemicals that regulate unique sets of genes in tissue-specific contexts are called selective estrogen receptor/androgen receptor modulators (SERMs/SARMs) [43, 44]. Ligands can also activate AR- and ERα-dependent transactivation through non-genomic “tethered”

T

signaling interactions in which the ligand bound receptor interacts with DNA-bound

IP

transcription factors [41, 42]. Ligand-independent signaling also exists. AR and ERα are

CR

activated through kinase-dependent cascades after exposure to cytokines (interleukin-6 for AR; interleukin-1β and tumor necrosis factor α for ERα) and growth factors (epidermal growth factor

US

(EGF) and insulin-like growth factor 1 for AR; EGF for ERα) that result in phosphorylation of

AN

the receptors [45, 46]. Lastly, AR and ERα can associate with the inner side of the cell membrane (membrane-bound receptors) and after ligand exposure activate other transcription

M

factors through signaling cascades [42]. The signaling pathway involving membrane-bound ERα

ED

is distinct from those regulated through the GPR30/GPER protein which belongs to the G-

PT

protein coupled class of receptors [47].

With these different mechanisms of regulation in mind, it is useful to consider what

CE

mechanisms of activation the biomarkers can and cannot detect. There is overwhelming evidence

AC

presented above that the biomarkers are useful to identify chemicals that work through classical receptor ligand-dependent mechanisms. This was by design, as the way in which the genes in the biomarker were identified would favor selection of a subset of genes that are commonly regulated despite differences in chemical-receptor-coregulator interactions. As an example of this, the ERα biomarker has been used to validate the activation of ERα by glyphosate-based herbicide constituents [48] and bisphenol A alternatives [49].

Page 16

ACCEPTED MANUSCRIPT The biomarkers can also identify conditions in which there are decreases in endogenous activators. Comparing cells grown in media depleted of steroids (charcoal-stripped media) versus cells grown in unstripped media gave gene expression profiles similar to exposure to ERα antagonists [9]. A number of related conditions were identified that led to suppression of AR

T

likely due to differences in the levels of activating steroids [10]. These included prostate cancer

IP

cell line xenografts grown in castrated mice versus intact mice, cells grown under androgen

CR

deprivation for up to 12 months versus no deprivation, and human prostate tissue from patients undergoing androgen ablation therapy (consisting of AR antagonists and androgen synthesis

US

inhibitors) versus before ablation therapy. These results indicate that the biomarkers can detect

AN

perturbation of AR and ERα by depletion of steroids and imply that chemicals that modulate endogenous steroid levels would secondarily affect receptor activity. An example of this was

M

discussed above (5α-reductase inhibition by finasteride). A variety of chemicals may also induce

ED

hepatic xenobiotic metabolizing enzymes, resulting in increased clearance of endogenous steroids, or in the case of animal models such as the rodent Hershberger assay, stimulating

PT

androgen metabolism, yielding “anti-androgenic” effects that are indirect in nature [38, 39]. It

CE

remains to be seen if such chemicals would result in changes in receptor activity that could be

AC

identified via the biomarker.

An additional mechanism that emerged from our biomarker studies is the effect of chemical- induced changes in receptor gene expression. The first indication of this came from comparisons between the ERα biomarker and profiles from cells knocked down by ESR1 siRNAs resulting in profiles that were similar to ERα antagonists [9]. Similarly, there were profiles in which AR was knocked down in prostate cancer cells that were similar to AR antagonists [10]. These findings indicate that the biomarkers can detect changes in the levels of

Page 17

ACCEPTED MANUSCRIPT ESR1 and AR expression and imply that chemicals that alter receptor expression could be identified in HTTr screens. In support of this, our AR screen [10] identified chemicals that act as AR antagonists in the absence of any evidence of direct binding. These chemicals may exert their effects at the level of AR expression. Treatment with retinoic acid receptor–related orphan

T

receptor (ROR) γ-specific antagonists, XY011 and SR2211 inhibit the expression of AR and

IP

regulated genes in prostate cancer cell lines and tumors [50]. Trichostatin A inhibited AR gene

CR

expression in prostate cancer cell lines [51, 52]. Urolithin A, derived from walnut pedunculagin, decreased expression of AR mRNA and protein [53]. Thus, our method identifies potential EDCs

US

that cause decreases in expression of the factor being studied. Examples of chemicals that

AN

increase the expression of the receptors and thus would act like agonists have yet to be identified.

M

The ability of our methods to reliably identify chemicals that modulate AR or ERα activity through activation or inhibition of key steps in steroid metabolism remains to be

ED

determined. The MCF-7 cell line may not be appropriate to identify chemicals that alter ERα

PT

activity through effects on aromatase (CYP19) which catalyzes the aromatization of testosterone to estradiol. Two aromatase inhibitors were included in the human breast cancer cell line

CE

compendium (letrozole, aminoglutethimide) but both had no effects on ERα [9] consistent with

AC

studies that show that endogenous aromatase is expressed at levels that do not allow inhibition effects to be observed [54]. To circumvent this problem, an MCF-7 cell line that constitutively expresses human aromatase, the MCF-7aro cell line has been used as a model to screen for aromatase inhibitors [55]. It would be useful to carry out a systematic study with a set of steroidogenesis inhibitors in cell lines used for HTTr to determine the number of potential targets that could contribute to AR or ERα modulation. In summary, our published screening methods identified 4 conditions in which AR and ERα were modulated: agonism, antagonism, depletion

Page 18

ACCEPTED MANUSCRIPT of activating steroids, and inhibition of receptor expression all of which work through the classical pathway. All of these conditions would also be likely observed as alterations in activity using trans-activation assays in cells that express the intact receptor. Identification of environmental chemicals that pose risks as EDCs have focused on those

IP

T

that interact directly with the receptors. The impact of environmental chemical exposure on activation of nonclassical pathways is an understudied area. Currently, there are no assays within

CR

the ToxCast battery (or for that matter used by other institutions/labs) that systematically

US

evaluate the effects of chemicals on these nonclassical pathways. Our methods for gene selection allows opportunities to create biomarkers specific for nonclassical pathways to fill these data

AN

gaps. For example, genes regulated by EGF through ERα could be identified by exposure of wild-type and ESR1-null MCF-7 cells to EGF at different time points or concentrations. After

M

removal of genes that overlap with those regulated through the classical pathway, the “EGF-

ED

ERα” biomarker could be used for screening HTTr data. A similar approach could be taken to

PT

create and use for screening a GPER-specific biomarker. This is increasingly pertinent as GPER is bound and activated by several EDCs [56] and when activated leads to a number of adverse

CE

phenotypes including E2-induced meiotic arrest of oocytes in fish [57] and increases in breast

AC

tumorigenesis and metastasis [47]. This set of ligand-dependent and -independent biomarkers would allow a systematic analysis of ED pathways mediated by AR, ERα, and GPER in current and future HTTr screens.

An integrated strategy for identification of potential ED chemicals in HTTr profiles

Page 19

ACCEPTED MANUSCRIPT Given the rapid advancements in the field, we can envision that in the near future HTTr profiling will be used as an alternative to some HTS assays, including those used in ToxCast [58], and eventually as an alternative to low-throughput in vitro methods or even animal tests. In the near-term scenario, HTTr profiling of chemicals in dose-response format will be carried out

T

as “Tier 0” assays prior to any ToxCast Tier 1 screening (Figure 3). Predictions using a battery

IP

of biomarkers will allow prioritization of further testing in targeted Tier 1 assays. For example,

CR

with AR and ERα positives, additional assays might be carried out to determine the mechanism of modulation (direct binding, alteration of metabolism, etc.) allowing a better assessment of

US

toxicological significance. Transcriptional points of departure for the genes in the biomarkers

AN

could be determined to estimate minimal concentrations of chemical that lead to activation or suppression of the factor similar to methods carried out to determine transcriptional points of

M

departure for pathways [59]. These values could then be used for in vitro to in vivo extrapolation

ED

to tissue concentrations of the chemical in humans. Exposure models could be used to determine the margin of exposure between estimated human exposures and those that induce responses,

PT

allowing a preliminary assessment of “safe” levels of exposure [28]. These in vitro to in vivo

CE

extrapolation (IVIVE) methods are already well established at the EPA and used to interpret ToxCast HTS data [60-64]. Quantitative IVIVE comparisons have demonstrated that in vitro

AC

approaches can be used to accurately predict effect levels from guideline animal studies [65], and could be similarly applied to activity concentrations derived from biomarkers. The large number of predictions coming from batteries of biomarkers including those that measure important targets of EDCs, could be integrated into the molecular initiating events (MIEs) and key events (KEs) in networks of AOPs allowing prioritization of chemicals preliminarily based on number of MIEs/KEs modulated and relevance of the adverse outcome

Page 20

ACCEPTED MANUSCRIPT pathways (AOPs) for impacting human health. Subcellular, cellular, and tissue adverse outcomes predicted from the integrated responses could then be confirmed by medium throughp ut

T

organotypic tissue assays, model organisms such as Zebrafish, or short-term rodent tests.

IP

Challenges for implementing a HTTr screening strategy

CR

There are a number of challenges that need to be overcome before HTTr screening can be routinely used to predict toxicological effects. One challenge for HTTr is choosing the best set of

US

cell lines for screening [66, 67]. The cell lines will determine the spectrum of potential chemical

AN

targets that can be measured dependent on the expression of not only the targets of interest but other associated proteins that determine target-dependent chemical responses. At least for EDC

M

screening, a set of cell lines that can measure both activation and suppression of the greatest

ED

number of the most important targets (AR, ERα, ERβ, thyroid hormone receptors, progesterone receptor, glucocorticoid receptor, etc.) would be most ideal. (Very few environmental chemicals

PT

appear to directly modulate thyroid hormone receptors but rather affect key steps in the

CE

synthesis, transport, and metabolism of thyroid hormone.) Ideally, cell lines would be selected for screening using a number of criteria including adequate expression of functional nonmutated,

AC

nonamplified genes encoding the most important xenobiotic targets. For example, in place of screening for ERα modulators in MCF-7 cells [9] and AR modulators in LAPC-4 cells [10], it would be more efficient to screen for ERα and AR modulators simultaneously in the same cell type. In fact, MCF-7 cells express AR and are responsive to androgens (although potentially not anti-androgens) [68]. An MCF-7 AR biomarker could be developed using a strategy similar to that discussed above involving comparisons between wild-type and AR-null cell lines. This strategy could be used to develop other biomarkers for HTTr screening in MCF-7 cells or other Page 21

ACCEPTED MANUSCRIPT cell lines. It is important to recognize that ideally, biomarkers can be designed to screen for factor modulation in any cell line. However, due to tissue-specific TF responses, there may not be enough regulated genes to create these cell-independent biomarkers. For example, the ERα biomarker works well in ERα-positive human breast cancer cell lines but not in ERα-positive cell

IP

T

lines from other human tissues (e.g., Ishikawa cells) [9]. An important criteria for choosing a set of cell lines would be the ability to screen for

CR

both agonists and antagonists in the same assay. While chemical agonists of ERα and AR were

US

easily identified, the MCF-7 and LAPC-4 cell lines diverged in their ability to simultaneously identify antagonists [9, 10]. Identification of both ERα agonists and antagonists in the same

AN

assay was most likely due to the low level constitutive activation of ERα and ERα-regulated

M

genes that could be reversed in the presence of an ERα antagonist [9] and likely due to the presence of estrogens in the media. As mentioned above, the depletion of steroids in the media

ED

gave profiles similar to ERα antagonists [9]. In contrast to MCF-7 cells, LAPC-4 cells had little,

PT

if any background activation of the AR-regulated genes. Thus for antagonist screening, the cells required cotreatment with the potential antagonist in the presence of a low level of an AR agonist

CE

[10]. Other prostate cancer cell lines may be more useful for simultaneous screening of agonists

AC

and antagonists. The prostate cancer cell line LNCaP is typically cultured in 10% fetal calf serum supplemented media which contains 55−98 pM of total testosterone [69]. LNCaP cells concentrate and metabolize testosterone to produce a physiological (i.e., approximately 10 nM) level of intracellular DHT which stimulates the growth of these cells in vitro [69] and likely contributes to the background level of AR activation and the ability to identify antagonists without adding a reference agonist. Thus, LNCaP cells might be useful to screen for both agonists and antagonists simultaneously in the absence of exposure to a reference agonist. The

Page 22

ACCEPTED MANUSCRIPT caveat is that the LNCaP cells possess an AR gene that contains a mutation in the ligand binding domain leading to promiscuous activation by some chemicals [36]. It might be possible to create a cell culture cocktail containing low levels of agonists of TFs that lack background activation ensuring that agonists and antagonists can be identified in the same HTTr screening assay.

IP

T

There are a number of technical challenges that need to be met before HTTr screening can be routinely used. Despite a number of recent technological advances in HTTr including

CR

assessment of gene expression in cells on 384-well plates, use of cell lysates in place of purified

US

RNA, and decreases in sequencing costs, performing HTTr on a large number of samples is still beyond the reach of a typical investigator. As mentioned above, DRUG-Seq is the latest

AN

technological innovation that may help to reduce the costs of HTTr [6]. The technique relies on

M

reducing the depth of genome coverage to ~2 million (M) sequencing reads/sample, far below that of a typical RNA-Seq experiment of >20M reads/sample. The difference in number of genes

ED

detected between a standard RNA-Seq of ~42M reads (17,000 genes) and ~2M reads (11,000

PT

genes) was not dramatic [6]. DRUG-Seq does not use methods such as attenuators that prevent highly expressed genes from occupying a disproportionate amount of available read space and

CE

decreasing the ability to quantify low abundance transcripts [7]. The question that remains

AC

unanswered is whether the level of genome coverage will allow assessment of all important pathways or will the generated profiles only be useful to examine pathways containing genes that are expressed at a sufficiently high level. As we found that decreasing the number of genes in the ERα biomarker decreases sensitivity of the predictions of both ERα activation and suppression [9], an insufficient sequencing depth may result in reduced predictive accuracy of biomarkers. A major challenge is how to comprehensively predict the effects of chemicals on all potential targets and pathways important in toxicity. This problem can be addressed with two Page 23

ACCEPTED MANUSCRIPT approaches viewed as either bottom up or top down. The bottom up approach involves a systematic effort to build, test and use for screening high value human biomarkers. As mentioned earlier, there are already efforts to create a large set of biomarkers from existing microarray data [15]. Given the accuracy of biomarkers developed using microarray profiles generated from cells

T

after chemical and genetic modulation (Table 1), it is worthwhile to consider strategies for

IP

building a large number of biomarkers that incorporate Crispr-Cas9 knockout of specific genes

CR

in human cell lines. These efforts would be similar to the wild-type vs. knockout mouse studies that resulted in a set of very accurate biomarkers used to assess chemical mode of action in livers

US

of treated mice [16-19, 21]. The complementary top down approach, uses systems biology

AN

approaches to build networks of interacting or co-expressed genes grouped into modules that often carry out similar functions and are coordinately regulated by chemical exposures that lead

M

to the same phenotypes. For example, the DrugMatrix dataset was used to build a network model

ED

of gene co-expression in the rat liver that could be linked to chemical-dependent and independent pathologies in not only rats but humans [70]. Other systems biology approaches

PT

could also be used that classify chemicals into different mechanism of action classes, such as was

CE

used to identify vascular disrupting compounds acting as developmental toxicants [71]. These systems approaches would help to generate hypotheses of mechanism of action through

AC

similarities between reference chemical behavior and those being characterized. There are significant time and resource requirements for building such models including the need for a large number of microarray profiles from the cell line or cell lines of interest. Furthermore, it will be important to determine whether the models can be used across different cell lines or whether individual models need to be built for each cell line. Regardless of the top-down or bottom-up approach, prediction of cell phenotypes (e.g., cell proliferation, apoptosis, mitochondrial

Page 24

ACCEPTED MANUSCRIPT toxicity) from altered gene expression would allow moving beyond predicting only MIEs to also include the downstream KEs mechanistically closer to the adverse outcomes. Currently, there are no examples of biomarkers in the literature that would be useful in HTTr screening to predict phenotypic effects in human cell cultures. Characterization of these biomarkers would require

T

parallel measurements of gene expression and alterations in cell phenotype, for example by high

IP

content imaging.

CR

While the Running Fisher test used with the biomarkers is consistently accurate, the

US

method is embedded in a commercial database that is not publicly available. Ideally, when an investigator generates microarray data and wishes to identify chemical targets, there is a

AN

publicly-available web interface to go to that would allow batch uploads of hundreds or

M

thousands of microarray profiles. These profiles could then be compared to batteries of biomarkers using a number of methods including an open source Running Fisher test. As an

ED

example of how this interface might look, a web tool has been recently described that allows

PT

comparison of uploaded microarray data with the TGx-DNA Damage Induction (DDI) biomarker used as an indicator of DNA damage [72]. Each uploaded gene list is compared to the

CE

expression of the 64 TGx-DDI biomarker genes in the 28 reference treatments used to create the

AC

biomarker. The analysis uses the Nearest Shrunken Centroids model, principal component analysis, and hierarchical clustering for comparison and gives a summary classification probability of whether the treatment causes DNA damage. In principle, this web tool could be used to simultaneously compare any gene list to a battery of characterized biomarkers. Another challenge is how to organize and integrate the predictions of the HTTr profiles into frameworks that allow linkage to adverse outcomes on a cellular, tissue, and organismal level. This problem is not unique to HTTr screening but is also relevant to ToxCast assays. Page 25

ACCEPTED MANUSCRIPT Currently there is little published work in this area. The AOPWiki (https://aopwiki.org/; accessed November 15, 2018) contains descriptions of 224 AOPs that cover a broad range of adverse outcomes in different tissues. Information about the key events (e.g., gene/protein target) can be extracted and in principle be used to create network graphs of AOPs linked by shared key

T

events. Recently, the AOP-DB, an exploratory database resource has been described that

IP

aggregates association relationships between genes and their related chemicals, diseases,

CR

pathways, species orthology information, ontologies, and gene interactions [73]. The AOP-DB

US

may be a useful interface between the AOP framework and existing and future HTS and HTTr screening data. In principle, key events could be populated by predictions from HTTr data

AN

allowing prioritization of chemicals that for example affect the greatest number of AOPs or those

M

deemed most toxicologically significant.

Lastly, the field could benefit from having large numbers of gene expression biomarkers

ED

and standardized computational methods accepted by regulatory agencies as tools for prediction

PT

in gene expression analysis. While there are gene expression biomarkers used to predict cancers including those that consist of multiple genes [74, 75], there are no gene expression biomarkers

CE

accepted by any regulatory agency specifically for toxicity testing. As genotoxicity testing is a

AC

critical part of chemical risk assessment, the Health and Environmental Sciences Institute's (HESI) Technical Committee for the application of genomics to mechanism‐ based risk assessment developed the TGx-DDI biomarker that can be used to classify chemicals as either DNA damage inducing (DDI) or non‐ DDI [76]. The biomarker can accurately classify 45 test agents across a broad set of chemical classes using the nCounter high‐ throughput cell‐ based testing platform and probability analysis, principal components analysis, and two‐ dimensional clustering [77]. In particular, the approach classified most compounds with irrelevant positive in Page 26

ACCEPTED MANUSCRIPT vitro chromosome damage as negative. The TGx-DDI biomarker has been recently submitted to the Food and Drug Administration to identify chemicals that cause DNA damage (https://www.fda.gov/downloads/Drugs/DevelopmentApprovalProcess/DrugDevelopmentTool sQualificationProgram/BiomarkerQualificationProgram/UCM605364.pdf ; accessed November

T

26, 2018). The HESI Genomics Committee is currently considering additional gene expression

IP

biomarkers that could be used for potential chemical characterization and classification in the

US

CR

context of short-term animal tests.

AN

Summary

Using our previously published work as a basis, we describe a gene expression

M

biomarker-based approach to screen for chemicals that modulate the activity of important targets

ED

of EDCs in HTTr screens. The approach can be extended to essentially measure most of the important targets in toxicology. To identify the biomarker genes, we used a WOE approach

PT

incorporating knowledge of TF-gene interactions including chemical and genetic perturbation

CE

experiments crucial to identify predictive genes. There are a number of advantages to using our WOE approach as compared to other methods such as machine learning [78]. The method

AC

is simple and does not require any specialized software or training to identify the genes. The approach is transparent as the genes are selected based on a predefined set of criteria. This contrasts to machine learning methods that are not transparent to the filtering criteria used to identify genes and the fact that filtering criteria cannot be readily altered by the user. Furthermore, the approach can be used to build predictive biomarkers from relatively small data sets as opposed to machine learning models that generally require larger numbers of

Page 27

ACCEPTED MANUSCRIPT samples to train the model. The resulting lists consist of a high percentage of genes that are known to be directly regulated by the factor. The biomarker genes can then be compared to essentially any statistically- filtered gene list that contains fold-change values. Overall, our studies consistently show that we can use a WOE approach to identify genes linked to the TF of

IP

T

interest. Our approach was shown to be accurate in identifying conditions of TF modulation. The

CR

ERα biomarker accurately replicated the results of a model based on 18 HTS ToxCast assays,

US

and likewise, the AR biomarker accurately replicated the results of a model based on 11 AR ToxCast assays. HTTr profiling for ERα and AR modulators, ideally in the same cell line, could

AN

complement the current screening paradigm by serving as Tier 0 screens which would be followed by more targeted Tier 1 assays to uncover the underlying mechanism of action. The

M

procedures, at least for ERα, have the advantage of simultaneously assessing agonist-like and

ED

antagonist- like activity in a single assay system. As detailed in a U.S. Federal Register Notice

PT

(https://www.federalregister.gov/articles/2015/06/19/2015-15182/use-of- high-throughput-assaysand-computational-tools-endocrine-disruptor-screening-program- notice; accessed 23 June 2015),

CE

three assays in the EDSP Tier I battery could be replaced by in vitro ERα assays based on the

AC

ability of the assays used in the Judson et al. model [30] to accurately predict uterotrophic results in mice and rats [79]. Our preliminary studies indicate that the ERα biomarker approach can also be used to accurately predict uterotrophic results (manuscript in preparation). Thus, the methods described here could not only be used as a more streamlined alternative to the ToxCast assays but also provide a general strategy for the identification of ERα and AR modulators that would meet the needs of EDSP stakeholders.

Page 28

ACCEPTED MANUSCRIPT Looking beyond the AR and ERα biomarkers and applying what we have learned so far, it is conceivable that future HTTr screening will utilize batteries of gene expression biomarkers that predict important targets in not only ED but other toxicities. This bottom up approach will complement more top down systems biology approaches to identify potential targets. A number

T

of challenges remain for HTTr, including how to systematically construct and characterize

IP

biomarkers that capitalize on more methodical approaches for gene knockout, chemical

CR

exposures, and gene expression; how to identify appropriate cell models for HTTr screening that can cover most, if not all important targets/pathways linked to adverse outcomes in a range of

US

developmental stages, tissues and species; how to predict the greatest number of targets at the

AN

lowest cost; how to move towards the use of computational systems of analysis that are transparent, publicly available, and can be used to easily compare generated gene lists to

M

biomarkers using a number of accepted methods. Finally, a major challenge for HTTr screening

ED

as well as Tox21/ToxCast assays is how to integrate the biomarker predictions into the growing network of MIEs and KEs in AOPs to allow linkage to adverse outcomes. Given the rapid

PT

advances in the field, it would not be surprising to see that many of these problems will be

AC

CE

solved in the next few years.

Acknowledgments

This study was carried out as part of the EPA High-throughput Testing project within the chemical safety for sustainability (CSS) program. The information in this document has been funded in part by the U.S. Environmental Protection Agency. The views expressed in this paper are those of the authors and do not necessarily reflect the statements, opinions, views, conclusions, or policies of the United States EPA. The authors declare they have no actual or Page 29

ACCEPTED MANUSCRIPT potential competing financial interests. We thank Mr. Chuck Gaul for assistance in creating the figures, and Drs. Imran Shah and Warren Casey for critical review of the manuscript.

T

Conflict of Interest

CR

IP

The authors declare no conflicts of interest.

References

5. 6. 7. 8.

9.

10. 11. 12. 13. 14.

US

AN

M

ED

4.

PT

3.

CE

2.

Diamanti-Kandarakis, E., et al., Endocrine-disrupting chemicals: an Endocrine Society scientific statement. Endocr Rev, 2009. 30(4): p. 293-342. Gore, A.C., et al., EDC-2: The Endocrine Society's Second Scientific Statement on EndocrineDisrupting Chemicals. Endocr Rev, 2015. 36(6): p. E1-e150. Judson, R., et al., In vitro and modelling approaches to risk assessment from the U.S. Environmental Protection Agency ToxCast programme. Basic Clin Pharmacol Toxicol, 2014. 115(1): p. 69-76. Merrick, B.A., Next generation sequencing data for use in risk assessment. Current Opinion in Toxicology, 2019(This issue). Bush, E.C., et al., PLATE-Seq for genome-wide regulatory network analysis of high-throughput screens. Nat Commun, 2017. 8(1): p. 105. Ye, C., et al., DRUG-seq for miniaturized high-throughput transcriptome profiling in drug discovery. Nat Commun, 2018. 9(1): p. 4307. Yeakley, J.M., et al., A trichostatin A expression signature identified by TempO-Seq targeted whole transcriptome profiling. PLoS One, 2017. 12(5): p. e0178302. Shah, I., Harrill J, Setzer RW, Haggard D, Karmaus A, Martin MT, Thomas RS. , Predicting Chemical Mechanisms of Action using High-Throughput Transcriptomic Data. . The Toxicologist, 2018. 162(1): p. 1689. Ryan, N., et al., Moving Toward Integrating Gene Expression Profiling Into High-Throughput Testing: A Gene Expression Biomarker Accurately Predicts Estrogen Receptor alpha Modulation in a Microarray Compendium. Toxicol Sci, 2016. 151(1): p. 88-103. Rooney, J.P., et al., Identification of Androgen Receptor Modulators in a Prostate Cancer Cell Line Microarray Compendium. Toxicol Sci, 2018. Iorio, F., et al., Discovery of drug mode of action and drug repositioning from transcriptional responses. Proc Natl Acad Sci U S A, 2010. 107(33): p. 14621-6. Woo, J.H., et al., Elucidating Compound Mechanism of Action by Network Perturbation Analysis. Cell, 2015. 162(2): p. 441-451. Biomarkers and surrogate endpoints: preferred definitions and conceptual framework. Clin Pharmacol Ther, 2001. 69(3): p. 89-95. Moreau, Y. and L.C. Tranchevent, Computational tools for prioritizing candidate genes: boosting disease gene discovery. Nat Rev Genet, 2012. 13(8): p. 523-36.

AC

1.

Page 30

ACCEPTED MANUSCRIPT

23. 24. 25. 26. 27. 28. 29. 30.

31. 32. 33. 34. 35. 36.

T

IP

CR

22.

US

21.

AN

20.

M

19.

ED

18.

PT

17.

CE

16.

Corton, J.C., Identification of Potential Chemical Carcinogens in Compendia of Gene Expression Profiles. The Toxicologist, 2018: p. 1690. Oshida, K., et al., Identification of chemical modulators of the constitutive activated receptor (CAR) in a gene expression compendium. Nucl Recept Signal, 2015. 13: p. e002. Oshida, K., et al., Screening a mouse liver gene expression compendium identifies modulators of the aryl hydrocarbon receptor (AhR). Toxicology, 2015. 336: p. 99-112. Oshida, K., et al., Identification of modulators of the nuclear receptor peroxisome proliferatoractivated receptor alpha (PPARalpha) in a mouse liver gene expression compendium. PLoS One, 2015. 10(2): p. e0112655. Oshida, K., et al., Disruption of STAT5b-Regulated Sexual Dimorphism of the Liver Transcriptome by Diverse Factors Is a Common Event. PLoS One, 2016. 11(3): p. e0148308. Oshida, K., D.J. Waxman, and J.C. Corton, Chemical and Hormonal Effects on STAT5b-Dependent Sexual Dimorphism of the Liver Transcriptome. PLoS One, 2016. 11(3): p. e0150284. Rooney, J., et al., Activation of Nrf2 in the liver is associated with stress resistance mediated by suppression of the growth hormone-regulated STAT5b transcription factor. PLoS One, 2018. 13(8): p. e0200004. Rooney, J.P., et al., Chemical Activation of the Constitutive Androstane Receptor (CAR) Leads to Activation of Oxidant-Induced Nrf2. Toxicol Sci, 2018. Kupershmidt, I., et al., Ontology-based meta-analysis of global collections of high-throughput public data. PLoS One, 2010. 5(9). Lamb, J., The Connectivity Map: a new tool for biomedical research. Nat Rev Cancer, 2007. 7(1): p. 54-60. Lamb, J., et al., The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science, 2006. 313(5795): p. 1929-35. Rooney, J.P., et al., From the Cover: Genomic Effects of Androstenedione and Sex-Specific Liver Cancer Susceptibility in Mice. Toxicol Sci, 2017. 160(1): p. 15-29. Rooney, J., et al., Adverse outcome pathway-driven identification of rat liver tumorigens in shortterm assays. Toxicol Appl Pharmacol, 2018. 356: p. 99-113. Sipes, N.S., et al., An Intuitive Approach for Predicting Potential Human Health Risk with the Tox21 10k Library. Environ Sci Technol, 2017. 51(18): p. 10786-10796. OECD, OECD Test No. 457: BG1Luc Estrogen Receptor Transactivation Test Method for Identifying Estrogen Receptor Agonists and Antagonists. OECD Publishing, Paris, 2012. Section 4. Judson, R.S., et al., Integrated Model of Chemical Perturbations of a Biological Pathway Using 18 In Vitro High-Throughput Screening Assays for the Estrogen Receptor. Toxicol Sci, 2015. 148(1): p. 137-54. Vivacqua, A., et al., The food contaminants bisphenol A and 4-nonylphenol act as agonists for estrogen receptor alpha in MCF7 breast cancer cells. Endocrine, 2003. 22(3): p. 275-84. Iskar, M., et al., Characterization of drug-induced transcriptional modules: towards drug repositioning and functional understanding. Mol Syst Biol, 2013. 9: p. 662. Biggar, R.J., Molecular pathways: digoxin use and estrogen-sensitive cancers--risks and possible therapeutic implications. Clin Cancer Res, 2012. 18(8): p. 2133-7. Kleinstreuer, N.C., et al., Development and Validation of a Computational Model for Androgen Receptor Activity. Chem Res Toxicol, 2017. 30(4): p. 946-964. Watson, P.A., V.K. Arora, and C.L. Sawyers, Emerging mechanisms of resistance to androgen receptor inhibitors in prostate cancer. Nat Rev Cancer, 2015. 15(12): p. 701-11. Veldscholte, J., et al., The androgen receptor in LNCaP cells contains a mutation in the ligand binding domain which affects steroid binding characteristics and response to antiandrogens. J Steroid Biochem Mol Biol, 1992. 41(3-8): p. 665-9.

AC

15.

Page 31

ACCEPTED MANUSCRIPT

45. 46. 47. 48. 49.

50. 51. 52. 53. 54.

55. 56.

57.

58.

T

IP

CR

44.

US

43.

AN

42.

M

41.

ED

40.

PT

39.

CE

38.

Makkonen, H., et al., Androgen receptor amplification is reflected in the transcriptional responses of Vertebral-Cancer of the Prostate cells. Mol Cell Endocrinol, 2011. 331(1): p. 57-65. Browne, P., et al., Development of a curated Hershberger database. Reprod Toxicol, 2018. 81: p. 259-271. Kleinstreuer, N.C., et al., Evaluation of androgen assay results using a curated Hershberger database. Reprod Toxicol, 2018. 81: p. 272-280. Lynch, C., et al., Identifying environmental chemicals as agonists of the androgen receptor by using a quantitative high-throughput screening platform. Toxicology, 2017. 385: p. 48-58. Lamont, K.R. and D.J. Tindall, Minireview: Alternative activation pathways for the androgen receptor in prostate cancer. Mol Endocrinol, 2011. 25(6): p. 897-907. Barone, I., L. Brusco, and S.A. Fuqua, Estrogen receptor mutations and changes in downstream gene expression and signaling. Clin Cancer Res, 2010. 16(10): p. 2702-8. Hah, N. and W.L. Kraus, Hormone-regulated transcriptomes: lessons learned from estrogen signaling pathways in breast cancer cells. Mol Cell Endocrinol, 2014. 382(1): p. 652-664. Narayanan, R., C.C. Coss, and J.T. Dalton, Development of selective androgen receptor modulators (SARMs). Mol Cell Endocrinol, 2018. 465: p. 134-142. Stender, J.D., et al., Structural and Molecular Mechanisms of Cytokine-Mediated Endocrine Resistance in Human Breast Cancer Cells. Mol Cell, 2017. 65(6): p. 1122-1135.e5. Lupien, M., et al., Growth factor stimulation induces a distinct ER(alpha) cistrome underlying breast cancer endocrine resistance. Genes Dev, 2010. 24(19): p. 2219-27. Prossnitz, E.R., GPER modulators: Opportunity Nox on the heels of a class Akt. J Steroid Biochem Mol Biol, 2018. 176: p. 73-81. Mesnage, R., et al., Evaluation of estrogen receptor alpha activation by glyphosate-based herbicide constituents. Food Chem Toxicol, 2017. 108(Pt A): p. 30-42. Mesnage, R., et al., Editor's Highlight: Transcriptome Profiling Reveals Bisphenol A Alternatives Activate Estrogen Receptor Alpha in Human Breast Cancer Cells. Toxicol Sci, 2017. 158(2): p. 431443. Wang, J., et al., ROR-gamma drives androgen receptor expression and represents a therapeutic target in castration-resistant prostate cancer. Nat Med, 2016. 22(5): p. 488-96. Rokhlin, O.W., et al., Mechanisms of cell death induced by histone deacetylase inhibitors in androgen receptor-positive prostate cancer cells. Mol Cancer Res, 2006. 4(2): p. 113-23. Welsbie, D.S., et al., Histone deacetylases are required for androgen receptor function in hormone-sensitive and castrate-resistant prostate cancer. Cancer Res, 2009. 69(3): p. 958-66. Sanchez-Gonzalez, C., et al., Urolithin A causes p21 up-regulation in prostate cancer cells. Eur J Nutr, 2016. 55(3): p. 1099-112. Zhou, D.J., D. Pompon, and S.A. Chen, Stable expression of human aromatase complementary DNA in mammalian cells: a useful system for aromatase inhibitor screening. Cancer Res, 1990. 50(21): p. 6949-54. Chen, S., et al., Cell-Based High-Throughput Screening for Aromatase Inhibitors in the Tox21 10K Library. Toxicol Sci, 2015. 147(2): p. 446-57. Thomas, P. and J. Dong, Binding and activation of the seven-transmembrane estrogen receptor GPR30 by environmental estrogens: a potential novel mechanism of endocrine disruption . J Steroid Biochem Mol Biol, 2006. 102(1-5): p. 175-9. Pang, Y. and P. Thomas, Involvement of estradiol-17beta and its membrane receptor, G protein coupled receptor 30 (GPR30) in regulation of oocyte maturation in zebrafish, Danio rario. Gen Comp Endocrinol, 2009. 161(1): p. 58-61. Thomas, R.S., et al., The Next Generation Blueprint of Computational Toxicology at the U.S. Environmental Protection Agency. Toxicol Sci, 2019. 169(2): p. 317-332.

AC

37.

Page 32

ACCEPTED MANUSCRIPT

67. 68. 69. 70. 71. 72.

73. 74.

75.

76.

77.

78.

T

IP

CR

US

66.

AN

65.

M

64.

ED

62. 63.

PT

61.

CE

60.

Phillips, J.R., et al., BMDExpress 2: Enhanced transcriptomic dose-response analysis workflow. Bioinformatics, 2018. Pearce, R., et al., httk: R Package for High-Throughput Toxicokinetics. Journal of Statistical Software, 2017. 79(4): p. 1-25. Rotroff, D.M., et al., Incorporating human dosimetry and exposure into high-throughput in vitro toxicity screening. Toxicol Sci, 2010. 117(2): p. 348-58. Wambaugh, J.F., et al., Toxicokinetic Triage for Environmental Chemicals. Toxicol Sci, 2015. Wetmore, B.A., et al., Incorporating High-Throughput Exposure Predictions With DosimetryAdjusted In Vitro Bioactivity to Inform Chemical Toxicity Testing. Toxicol Sci, 2015. 148(1): p. 121-36. Wetmore, B.A., et al., Relative impact of incorporating pharmacokinetics on predicting in vivo hazard and mode of action from high-throughput in vitro toxicity assays. Toxicol Sci, 2013. 132(2): p. 327-46. Casey, W.M., et al., Evaluation and Optimization of Pharmacokinetic Models for in Vitro to in Vivo Extrapolation of Estrogenic Activity for Environmental Chemicals. Environ Health Perspect, 2018. 126(9): p. 97001. Judson, R., Order from Chaos: Pattern Recognition in Challenging Human Health Datasets. Toxicological Sciences, 2019. xx: p. xx. Setzer, R., Thomas, RS, Strategies for Choosing Biological Diversity for High Throughput Toxicity Testing. Toxicological Sciences, 2016. xx: p. xx. Yeh, S., et al., Abnormal mammary gland development and growth retardation in female mice and MCF7 breast cancer cells lacking androgen receptor. J Exp Med, 2003. 198(12): p. 1899-908. Sedelaar, J.P. and J.T. Isaacs, Tissue culture media supplemented with 10% fetal calf serum contains a castrate level of testosterone. Prostate, 2009. 69(16): p. 1724-9. Sutherland, J.J., et al., Toxicogenomic module associations with pathogenesis: a network-based approach to understanding drug toxicity. Pharmacogenomics J, 2018. 18(3): p. 377-390. Kleinstreuer, N., et al., A computational model predicting disruption of blood vessel development. PLoS Comput Biol, 2013. 9(4): p. e1002996. Jackson, M.A., et al., The TGx-28.65 biomarker online application for analysis of transcriptomics data to identify DNA damage-inducing chemicals in human cell cultures. Environ Mol Mutagen, 2017. 58(7): p. 529-535. Pittman, M.E., et al., AOP-DB: A database resource for the exploration of Adverse Outcome Pathways through integrated association networks. Toxicol Appl Pharmacol, 2018. 343: p. 71-83. Kamel, H.F.M. and H. Al-Amodi, Exploitation of Gene Expression and Cancer Biomarkers in Paving the Path to Era of Personalized Medicine. Genomics Proteomics Bioinformatics, 2017. 15(4): p. 220-235. Martinez-Ledesma, E., R.G. Verhaak, and V. Trevino, Identification of a multi-cancer gene expression biomarker for cancer clinical outcomes using a network-based algorithm. Sci Rep, 2015. 5: p. 11966. Li, H.H., et al., Development of a toxicogenomics signature for genotoxicity using a doseoptimization and informatics strategy in human cells. Environ Mol Mutagen, 2015. 56(6): p. 50519. Li, H.H., et al., Development and validation of a high-throughput transcriptomic biomarker to address 21st century genetic toxicology needs. Proc Natl Acad Sci U S A, 2017. 114(51): p. E10881-e10889. Libbrecht, M.W. and W.S. Noble, Machine learning applications in genetics and genomics. Nat Rev Genet, 2015. 16(6): p. 321-32.

AC

59.

Page 33

ACCEPTED MANUSCRIPT

T IP CR US AN M ED PT

81.

CE

80.

Browne, P., et al., Screening Chemicals for Estrogen Receptor Bioactivity Using a Computational Model. Environ Sci Technol, 2015. 49(14): p. 8804-14. Corton, J.C., A. Williams, and C.L. Yauk, Using a gene expression biomarker to identify DNA damage-inducing agents in microarray profiles. Environ Mol Mutagen, 2018. Rooney JP, C.J., Mining a Human Transcriptome Database for Nrf2 Modulators. The Toxicologist, 2016. 2624.

AC

79.

Page 34

ACCEPTED MANUSCRIPT

Table 1. Characteristics of gene expression biomarkers for HTTr screening of human cells. Key event predicted

Domain of applicability

Number of genes

Chemical or other filter

Genetic filter

ERα

46

17α-estradiol, clomifene, diethylstilbestrol, 17β-estradiol, estrone, fulvestrant, genistein, hexestrol, mestranol, raloxifene

None

51

R1881, RT I-6413-018, androstenedione, dexamethasone, bicalutamide, triticonazole, vinclosalin, tebuconazole

64

28 treatments4

Wild-type vs. AR knockdown or AR constitutive activation mutants None

Nrf2

ERα positive human breast cancer cell lines AR positive human prostate cancer cell lines p53 positive human cancer cell lines Many cell lines

143

oltipraz, sulforaphane, CDDO-Im

NFkB

Many cell lines

108

T umor necrosis factor α

AR T Gx-DDI

Wild-type vs. Keap1 knockdown; Wild-type vs. Nrf2 knockdown Wild-type vs. IkB overproduction

A

U N

1

M

Predictive accuracy inhibition 98%

Reference

97%

98%; 93%2

[10]

87-97%3

Not applicable

[80]

92%

Not applicable Not applicable

[81]

T P

I R

C S

Predictive accuracy activation 95%1

92-96%

[9]

In preparation

For the ER biomarker, used 114 chemicals that allowed a comparison to the ER computational model. For the AR biomarker predictions, the first prediction of accuracy of the suppression was for all 158 chemicals tested in a number of prostate cancer cell lines and the second was for 28 chemicals tested in LAPC-4 cells. 3 For the TGx-DDI biomarker which predicts effects of chemicals on p53, accuracies are shown as a range dependent on the threshold for determining DDI positives. 4 The 28 treatments used to create the TGx-DDI biomarker are found in [77]. 2

D E

T P

E C

C A

Page 35

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

Figure 1. Construction of gene expression biomarkers that predict targets of potential EDCs.

D E

M

Biomarker genes are identified from the differentially-expressed genes in experiments in which cells are exposed to chemicals known

T P

to either activate or inhibit the TF important in ED. Additional microarray comparisons help to identify TF-dependent genes, including those in which the TF is either knocked down by siRNA, knocked out by Crispr-Cas9 technologies, or activated by mutations or

E C

overexpression. After generation of statistically- filtered gene lists, the biomarker genes are derived using a number of filtering steps

C A

(see text). Final biomarkers include lists of genes and associated fold-change values which represent the average expression across the agonist/activator comparisons. Post hoc analysis includes examination of the linkage of the genes to pathways and determination if they are primary targets of the TF by queries of chromatin immunoprecipitation coupled with DNA sequencing (ChIP-Seq) datasets.

Page 36

ACCEPTED MANUSCRIPT

T P

I R

C S

U N

A

D E

M

T P

E C

C A

Figure 2. Approaches to testing biomarker accuracy and use in in silico screening.

Page 37

ACCEPTED MANUSCRIPT (Top) Comparisons are made between the list of biomarker genes and each microarray comparison in the database found in BaseSpace Correlation Engine (Illumina) using a pair-wise rank-based algorithm (the Running Fisher test). The test results in direction of correlation and p-value for each comparison. Predictive accuracy of the biomarker is determined by comparison to profiles from cells

T P

exposed to chemicals with known effects on the factor to be analyzed. Chemicals are categorized as those that activate (+), suppress (-

I R

), or do not affect (o) the TF and used in the accuracy test which ranks the comparisons by -Log(p-value)s of the correlations. An

C S

example of the ranking and associated heatmap of biomarker genes is shown.

U N

(Bottom) To perform in silico screening to identify novel chemicals or other factors that modulate the TF, the biomarker genes are

A

compared to the profiles from cells exposed to uncharacterized chemicals or factors. Predictions can then be confirmed using other

M

techniques including comparisons in wild-type and TF-null cell lines.

D E

T P

E C

C A

Page 38

ACCEPTED MANUSCRIPT Figure 3. Putting biomarker predictions into networks of adverse outcome pathways. In ongoing and future HTTr profiling studies, cells are exposed to chemicals or environmental mixtures at different concentrations. Gene expression changes are compared to a large battery of validated biomarkers, each of which can accurately predict the activation

T P

or suppression of not only TFs and other adverse outcomes important in ED but effects on pathways (e.g., tumor necrosis factor α) and

I R

phenotypes (e.g., mitochondrial dysfunction). Predictions can then be used to populate network models of AOPs that cover a broad

C S

range of toxicities and mechanisms and used to rank and prioritize chemicals or mixtures for validation of TF targets in a series of

U N

wild-type and TF-null cell lines.

A

D E

M

T P

E C

C A

Page 39

ACCEPTED MANUSCRIPT Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

T P

I R

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

C S

A

U N

D E

M

T P

E C

C A

Page 40