In vitro transcriptomic analyses reveal pathway
perturbations, estrogenic activities, and potencies
of data-poor BPA alternative chemicals
Geronimo Matteo ,1,2 Karen Leingartner,1,2 Andrea Rowan-Carroll,1,2 Matthew Meier,1,2 Andrew Williams,1,2 Marc A. Beal,3,4

Matthew Gagn�e,4 Reza Farmahin,4 Shamika Wickramasuriya,4 Anthony J. F. Reardon,4 Tara Barton-Maclaren,4

J. Christopher Corton ,5 Carole L. Yauk ,2,* Ella Atlas 1,6,*
1Environmental Health Science and Research Bureau, Healthy Environments and Consumer Safety Branch (HECSB) Health Canada, Ottawa, Ontario K2K 0K9,
Canada
2Department of Biology, Faculty of Science, University of Ottawa, Ottawa, Ontario K1N 9A7, Canada
3Bureau of Chemical Safety, Health Canada
4Existing Substances Risk Assessment Bureau, Healthy Environments and Consumer Safety Branch, Health Canada, Ottawa, Ontario K2K 0K9, Canada
5Center for Computational Toxicology and Exposure, US Environmental Protection Agency, Research Triangle Park, North Carolina, USA
6Department of Biochemistry, Faculty of Medicine, University of Ottawa, Ottawa, Ontario K1H 8M5, Canada

*To whom correspondence should be addressed. E-mail: ella.atlas@hc-sc.gc.ca and carole.yauk@uottawa.ca.

Disclaimer: 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. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.

Abstract

Since initial regulatory action in 2010 in Canada, bisphenol A (BPA) has been progressively replaced by structurally related alternative
chemicals. Unfortunately, many of these chemicals are data-poor, limiting toxicological risk assessment. We used high-throughput
transcriptomics to evaluate potential hazards and compare potencies of BPA and 15 BPA alternative chemicals in cultured breast
cancer cells. MCF-7 cells were exposed to BPA and 15 alternative chemicals (0.0005–100 mM) for 48 h. TempO-Seq (BioSpyder Inc) was
used to examine global transcriptomic changes and estrogen receptor alpha (ERa)-associated transcriptional changes. Benchmark
concentration (BMC) analysis was conducted to identify 2 global transcriptomic points of departure: (1) the lowest pathway median
gene BMC and (2) the 25th lowest rank-ordered gene BMC. ERa activation was evaluated using a published transcriptomic biomarker
and an ERa-specific transcriptomic point of departure was derived. Genes fitting BMC models were subjected to upstream regulator
and canonical pathway analysis in Ingenuity Pathway Analysis. Biomarker analysis identified BPA and 8 alternative chemicals as ERa
active. Global and ERa transcriptomic points of departure produced highly similar potency rankings with bisphenol AF as the most
potent chemical tested, followed by BPA and bisphenol C. Further, BPA and transcriptionally active alternative chemicals enriched
similar gene sets associated with increased cell division and cancer-related processes. These data provide support for future read-
across applications of transcriptomic profiling for risk assessment of data-poor chemicals and suggest that several BPA alternative
chemicals may cause hazards at similar concentrations to BPA.

Keywords: TempO-Seq; BPA; biomarker; benchmark concentration; potency; toxicogenomics

Bisphenol A (BPA) is an endocrine disruptor that interacts with
nuclear hormone receptors, as well as altering non-hormonal
pathways, to exert its effects (Cimmino et al., 2020). Exposure to
BPA is associated with multi-organ toxicity and negative health
outcomes including metabolic and endocrine dysfunction, repro-
ductive and developmental disorders, and hormone-related can-
cer (Ma et al., 2019). Canada was the first country to take
regulatory action to limit exposure to BPA in 2008, followed by
the United States, the EU, and others (Rogers, 2021). As pressures
to phase out BPA increase, manufacturers are relying more on
BPA alternative chemicals.

BPA alternative chemicals are detected in the environment,
consumer products, and in humans (Chen et al., 2016). Many of
these chemicals are structurally similar to BPA but differ in their
adjoining atom and aryl substituents. A systematic review has

highlighted the ability of bisphenol S (BPS), a sulfone-based
bisphenol, and bisphenol F (4,40-BPF), a BPA analogue missing the
geminal methyl groups, to act as endocrine disruptors at similar
concentrations to BPA (Rochester and Bolden, 2015). Some substi-
tutes are less structurally related to BPA, like Pergafast201 and
4,400-bis-(p-tolylsulfonyl ureido)-diphenylmethane (BTUM); these
are non-phenolic sulfonyl urea-based alterative chemicals used
primarily in thermal paper (Björnsdotter et al., 2017). Several of
these alternative chemicals are also associated with endocrine
disruption, reproductive toxicity, and carcinogenicity (den
Braver-Sewradj et al., 2020). Thus, there is a clear need to identify
potential hazards and determine relative potencies of BPA alter-
native chemicals.

High-throughput transcriptomics (HTTr) is increasingly used
to identify hazards and establish bioactivity thresholds for

VC His Majesty the King in Right of Canada, as represented by the Minister of Health, 2022.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which
permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Toxicological Sciences, 2023, 191(2), 266–275

https://doi.org/10.1093/toxsci/kfac127
Advance Access Publication Date: 19 December 2022

Research article

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://orcid.org/0000-0003-0819-4471
https://orcid.org/0000-0002-6197-2036
https://orcid.org/0000-0002-6725-3454
https://orcid.org/0000-0003-4410-4402


regulatory toxicology (Thomas et al., 2019). We recently applied
HTTr to compare potencies of BPA alternative chemicals using
benchmark concentration (BMC) analysis in H9 human
embryonic stem cells exposed to BPA, BPS, 4,40-BPF, and 3,30,5,50-
tetrabromobisphenol A (TBBPA) at a range of concentrations
(1–700 mM) (Peshdary et al., 2021). BPA, BPS, and 4,40-BPF had over-
lapping median BMCs (approximately 100 mM), suggesting similar
potencies, whereas the BMC of TBBPA was much higher. These
data must be interpreted with caution as the absence of ERs in
human embryonic stem cells underestimates the toxicity of BPA
and similar compounds as the primary mode action of these
chemicals is via these nuclear hormone receptors.

An important approach to screening BPA alternatives is to spe-
cifically compare potencies based on interactions with the estrogen
receptor (ER). For example, Kitamura et al. (2005) tested the effects
of several BPA alternatives using MCF-7 cells transfected with an
estrogen response element (ERE) luciferase reporter. Bisphenol AF
(BPAF) and bisphenol B (BPB) stimulated estrogenic activity at lower
concentrations than BPA, whereas BPS and 4,40-BPF required higher
concentrations. Mesnage et al. (2017) found that BPAF, a hexafluoro
bisphenol, was the most potent BPA alternative in stimulating an
ERE-mediated luciferase reporter gene relative to BPB, bisphenol Z
(BPZ), BPA, 4,40-BPF, bisphenol AP (BPAP), and BPS (Mesnage et al.,
2017). The same study also conducted whole genome transcrip-
tomic analysis with MCF-7 cells and observed that these BPA alter-
natives perturbed Gene Ontology (GO) gene sets associated with
cell cycle, steroid hormone response, and breast cancer.

Nuclear receptor activation can also be assessed through the
analysis of transcriptomic biomarkers. Toward this, a 46 gene
transcriptomic biomarker associated with ERa was developed
and validated in ERa-positive MCF-7 cells (Ryan et al., 2016). This
biomarker predicts compounds that activate ERa through statisti-
cal comparisons of gene expression profiles of the 46 biomarker
genes, with a pattern of 32 upregulated and 14 downregulated
genes following ER activation. This approach can be used to iden-
tify the potential for the chemicals to interact with ERa and pro-
duce a composite BMC score for the biomarker to compare
potencies for this specific molecular initiating event.

In this study, we explored the use of HTTr to establish mecha-
nistic similarities, compare potencies, and increase understanding
of data-poor BPA alternative chemical toxicity in support of prioriti-
zation and assessment activities by Health Canada. MCF-7 cells
were exposed to BPA and 15 BPA alternatives (from 0.0005 to
100mM) or solvent (0.1% DMSO) for 48 h. The exposure concentra-
tions chosen span several orders of magnitude above and below the
mean concentration of BPA in human urine (approximately 10 nM;
Chen et al., 2016) and are consistent with those used in the US
Environmental Protection Agency’s ToxCast Program. HTTr using
Templated Oligo-Sequencing (TempO-Seq; BioSpyder Inc) was used
to examine both general toxicological transcriptional effects as well
as ERa-associated transcriptional changes induced by these chemi-
cals. BMC analysis was conducted to identify both general and ERa-
specific transcriptomic points of departure (tPODs). Finally, genes
showing robust concentration-response were subjected to
upstream regulator and pathway analysis to explore the toxicody-
namic similarity of these chemicals to each other and BPA.

Materials and methods
Chemicals

See Table 1 for a list of chemicals and suppliers. The final con-
centration of DMSO in media was 0.1% for all chemicals tested.

Cell culture

MCF-7 cells (ATCC, Manassas, Virginia) were routinely cultured in
100mm dishes (Corning Falcon No. 353003, Corning, New York)
with Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Thermo
Fisher Scientific, Waltham, Massachusetts) supplemented with 10%
fetal bovine serum (FBS; Wisent Bioproducts, Saint-Jean Baptiste,
QC, Canada), and 1% penicillin and streptomycin, at 37�C in a
humidified atmosphere with 5% CO2. For the experiment, cells were
seeded into clear 96-well plates (Corning Falcon No. 353075) at a
density of 2.5� 104 cells per well in phenol red-free DMEM supple-
mented with 5% charcoal-dextran stripped FBS (CD-FBS) to elimi-
nate estrogenic components (Atlas et al., 2003). The next day the
medium was removed and replaced with DMEM with 5% CD-FBS
containing the chemicals of interest. Cells were exposed, in quadru-
plicate, to BPA and 15 alternative chemicals at 10 different concen-
trations (0.0005, 0.001, 0.01, 0.1, 0.5, 1, 5, 10, 50, 100mM), the positive
control 17b-estradiol (E2; 0.0001, 0.001, 0.01, 0.1, 1, 10nM), dexame-
thasone (Dex—a non-ER-interacting control; 0.0001, 0.001, 0.01, 0.1,
1mM) alongside solvent controls (0.1% DMSO; 4 solvent controls per
plate) for 48 h. The exposure timeframe was selected based on pilot
studies to ensure that we observed a robust response in estrogen-
responsive genes (data not shown). After exposure, cells were
washed with phosphate-buffered saline (PBS) and lysed in situ using
2� TempO-Seq lysis buffer diluted with an equal amount of PBS for
Templated Oligo-Sequencing (TempO-Seq). We observed precipitate
in some chemicals and removed these concentrations from further
analysis: BPAF (50, 100mM), BPC (50, 100mM), BADGE (50, 100mM), D-
8 (50, 100mM), DDS (10, 50, 100mM), BTUM (50, 100mM), TGSA (50,
100mM), BPS-MPE (50, 100mM).

Cell viability assay

To determine cytotoxic concentrations, MCF-7 cells were seeded
into black 96-well plates (Corning Falcon No. 3603) and treated as
described above. Cell viability was measured using the CellTiter-
Blue Cell Viability Assay (Promega Corp, Madison, Wisconsin), as
per the manufacturer’s instructions. This is a fluorometric assay
that measures the metabolic capacity of cells based on resazurin
reduction. Briefly, after 48-h exposure, CellTiter-Blue reagent was
added to each well, and the plates were incubated at 37�C, 5%
CO2 for 2 h. The fluorescence was read using excitation wave-
length of 560 nm and emission of 590 nm using a SpectraMax M2
(Molecular Devices LLC, San Jose, California). Fluorescence read-
ings were expressed as a % ratio to that of the DMSO control
group. Cytotoxicity was defined as readings <50% of the control.

Cell proliferation assay

MCF-7 cells were seeded at a density of 2.0�104 cells/well in clear
24-well plates (Corning Costar No. 3524) and incubated at 37�C
and 5% CO2 in phenol red free DMEM/F12 media containing 5%
CD-FBS and 1% penicillin-streptomycin. Cells were treated with
BPA and 15 alternative chemicals at 11 different concentrations
(0.0001–100 mM) and positive control E2 (0.0001–10 nM). Chemical
concentrations that yielded precipitates were again excluded
from this analysis. Cells were counted on day 7 using a TC20
automated cell counter (Bio-Rad, Hercules, California) as per
manufacturer’s instructions after being treated every 2 days with
indicated concentration of the bisphenols and E2.

TempO-Seq library building and next-generation
sequencing

Gene expression was measured using the TempO-Seq Human
Whole Transcriptome v2.0 kit (BioSpyder Technologies Inc,

Matteo et al. | 267

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023



Carlsbad, California) as per manufacturer’s instructions and previ-
ously described (Yeakley et al., 2017). Cell lysates and positive tech-
nical controls (Human Universal Reference RNA—uhrRNA Agilent
Cat No. 740000, Santa Clara, California, and Human Brain Total
RNA brRNA—ThermoFisher AM7962, Waltham, Massachusetts),
as well as no-cell negative controls (1� TempO-Seq lysis buffer
alone) were hybridized to the Detector Oligo (DO) Pool in an
annealing kit for the whole human genome supplied by BioSpyder
for 10 min at 70�C followed by a temperature gradient with a ramp
rate of 0.5�C/min to 45�C over 50 min with a 16-h hold at 45�C and
then cooled to 25�C. Nuclease digestion was employed to remove
excess, unbound, or incorrectly bound DOs at 37�C for 90 min.
Amplification templates were generated by ligating DO pairs
bound to adjacent target sequences for 1 h at 37�C, followed by
enzyme denaturation for 15 min at 80�C. Amplification templates
(10mL) were pipetted into a 96-well PCR Pre-Mix and Primer plate
supplied by BioSpyder and amplified using a CFX96 Real-Time PCR
Detection System (Bio-Rad) to incorporate a unique sample index/
tag sequence and the sequencing adaptors for each sample. The
following PCR settings were used: 37�C for 10 min, 95�C for 1 min;
25 cycles of 95�C for 10 s, 65�C for 30 s, 68�C for 30 s (with optical
read for visual sample quality control [QC]); 68�C for 2 min; and
hold at 25�C prior to library pooling and purification. For a list of
attenuators used, see Supplementary File 1.

NucleoSpin Gel and PCR Clean-up kits were used to pool and
purify labeled amplicons. Next-generation sequencing (NGS)
libraries were sequenced using a NextSeq 500 High-Throughput
Sequencing System (Illumina, San Diego, California), using 50
cycles from a 75-cycle high-throughput flow cell to achieve a
median read depth of 2 million reads per sample. Data were proc-
essed as described below and reads were aligned to the BioSpyder
TempO-Seq Human Whole Transcriptome probe set (22 537
probes over 19 687 genes) using their purpose-built pipeline.

Data processing and QC

Reads were demultiplexed from the BCL files and processed into
FASTQ files using bcl2fastq v2.20.0.422. The FASTQ files were
processed with the TempO-SeqR analysis pipeline in R (version
3.1) supplied by BioSpyder. This script uses STAR v2.7.8a (Dobin
et al., 2013) to perform alignment of raw reads to the reference

sequence and the qCount function of the QuasR R package
(v1.30.0; Gaidatzis et al., 2015) to extract the feature counts from
the aligned reads (BAM files) using features specified in a GTF
file. The script generates a samples-by-probes count matrix.

Study-wide QC was performed on the count matrix using sev-
eral methods to measure consistency and remove low-quality
samples, using the methods in Harrill et al. (2021) as a guideline.
A cutoff of 0.1 for 1-Spearman’s q was used to remove samples
that were not correlated with others in the study. As in Harrill
et al. (2021), we also used a 10% cut-off of uniquely mapped reads
as the number of target sequences (eg, 100 000 reads to pass filter
when the target is 1 000 000 for TempO-Seq experiments). We
removed any samples outside of Tukey’s Outer Fence (3� inter-
quartile range) for: (1) the number of probes capturing the top
80% of the signal, and (2) the number of detected probes (those
with at least 5 mapped reads). Samples with a Gini coefficient
(which measures inequality in distributions) greater than 0.95
were excluded. Samples removed by these criteria are listed in
Supplementary File 2. We note that a technical error led to loss of
the 0.001 and 0.01 mM concentrations of BPA-exposed cells (see
Supplementary Table 2).

To create a matrix for biomarker analysis, individual pairwise
contrasts for each concentration and each chemical tested were
created to the respective 0.1% DMSO control samples for each plate.
Following the recommendations set out by the Omics Data Analysis
Frameworks for Regulatory application (R-ODAF) guidelines
(Verheijen et al., 2022), genes were filtered for each contrast tested
to include only those where 75% of at least one experimental group
were above 0.5 counts per million (CPM), and spurious spikes were
removed in which [max–median] of counts were less than [sum of
counts]/[number of replicatesþ1]. We used DESeq2 1.30.0 (Love
et al., 2014) to estimate fold changes and normalize for library size
within the TempO-Seq data. The ashr method was used to perform
log2FoldChange shrinkage (Stephens, 2017). The code used to per-
form processing of high-throughput sequencing data is available at
https://github.com/R-ODAF (last accessed December 2022).

Gene expression biomarker analysis

To determine if the BPA alternatives activated ER or stress path-
ways, the expression profile of each chemical-concentration

Table 1. List of chemicals and suppliers

Chemical name Abbreviation CASRN Purity (%) Supplier

2,2-Bis(4-hydroxyphenyl)propane BPA 80-05-7 >99 Sigma-Aldrich Inc
(St Louis, Missouri)Bis(4-hydroxyphenyl)methane 4,40-BPF 620-92-8 >98

1,1-Bis(4-hydroxyphenyl)-1-phenylethane BPAP 1571-75-1 99
4,40-Sulfonyldiphenol BPS 80-09-1 98
2,2-Bis[4-(glycidyloxy)phenyl]propane BADGE 1675-54-3 98
4,40-Dichlorodiphenyl sulfone DDS 80-07-9 98
4,40-(Hexafluoroisopropylidene)diphenol BPAF 1478-61-1 97
Dimethyl sulfoxide DMSO 67-68-5 >99
17b-estradiol E2 50-28-2 >98
Dexamethasone Dex 50-02-2 98
4-[[4-(2-Propen-1-yloxy)phenyl]sulfonyl]phenol BPS-MAE 97042-18-7 97 Toronto Research

Chemicals (Toronto,
ON, Canada)

4,400-Bis-(p-tolylsulfonyl ureido)-diphenylmethane BTUM 151882-81-4 98
4-((4-Isopropoxyphenyl)sulfonyl)phenol D-8 95235-30-6 98
4-Methyl-N-[[[3-[[(4-methylphenyl)sulfonyl]

oxy]phenyl]amino]carbonyl]benzenesulfonamide
Pergafast201 (P201) 232938-43-1 98

2,40-Dihydroxydiphenyl sulfone 2,40-BPS 5397-34-2 98
4-(4-Hydroxy-3-prop-2-enylphenyl)sulfonyl-2-prop-2-enylphenol TGSA 41481-66-7 98
2,2-Bis(4-hydroxy-3-methylphenyl)propane BPC 79-97-0 98 Accustandard (New

Haven, Connecticut)
2,40-Dihydroxydiphenylmethane 2,40-BPF 2467-03-0 >98 TCI America (Portland,

Oregon)
4-((4-(Benzyloxy)phenyl)sulfonyl)phenol BPS-MPE 63134-33-8 >98

268 | In Vitro Transcriptomic Analyses

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://github.com/R-ODAF


tested relative to solvent control (derived from the analysis of
biosets) was compared with a number of characterized bio-
markers as previously described (Kupershmidt et al., 2010).
Biosets are gene expression profiles of genes filtered by unad-
justed p value <.05 and absolute fold change �1.2. The bio-
markers included those that predict modulation of ER (Ryan et al.,
2016), Nrf2 (Rooney et al., 2020), heat shock factor 1 (Cervantes
and Corton, 2021), metal-induced transcription factor 1 (Jackson
et al., 2020), TGx-DDI (Li et al., 2019), TGx-HDACi (Cho et al., 2021),
and NF-kB (Korunes et al., 2022). The correlations between each
biomarker and the gene lists of BPA and alternatives were deter-
mined using the Running Fisher algorithm as described previ-
ously in the BaseSpace Correlation Engine (Ryan et al., 2016). The
Running Fisher algorithm provides an assessment of the statisti-
cal significance of the correlation of the overlapping genes
between the biomarker and each gene list providing a summary p
value. A complete description of the Running Fisher test is pro-
vided in Kupershmidt et al. (2010). The results were exported and
each p value was converted to a �log(p value). Negative values
were used to indicate negative correlation between the biomarker
and the gene list. Thresholds for significance were set at �log(p
value) �4 for activation or ��4 for suppression based on prior
studies.

Points of departure analysis

We prefiltered probes using a Williams trend test (p< .05) and
absolute fold change >1.5. Probes with the following criteria were
removed: (1) having a BMC greater than the concentration used
in the analysis after excluding cytotoxic/precipitating concentra-
tions; (2) mapping to more than one gene; (3) having a model fit p
value <.1 determined by a likelihood ratio test; and (4) having a
BMC upper (BMCU) to BMC lower (BMCL) ratio <40.

BMDExpress2.3 was used to derive BMCs for BPA alternative
chemicals as previously described (Farmahin et al., 2017).
Concentration-response modeling was done by fitting genes best
fit models, including polynomial 2�, linear, power (power term
constrained to �1), and exponential models (degrees 2–5); best fit
models were selected for each probe based on the lowest Akaike
Information Criterion (AIC). The BMDExpress parameters were:
250 maximum iterations, 0.95 confidence level, assume constant
variance, BMR type SD, BMR factor 1 SD, restrict power �1. The
BMDExpress model selection criteria were: compute and utilize in
best model selection, best poly model test Nested Chi Square, p
value cut-off .05. Probes that met all the BMC filtering criteria
were converted to their corresponding Entrez Identifiers. Data
were further analyzed using R-4.0.5 statistical software. Gene
accumulation plots were made by plotting the BMC of each gene
on the horizontal axis and representing the gene accumulation
number along the vertical axis.

Two approaches were used for tPOD generation intended for:
(1) general toxicity; and (2) ERa-specific response.

To establish general toxicity tPODs 2 methods were used:

a) Lowest median pathway BMC: We produced the lowest
median pathway BMC by aligning genes and their associ-
ated BMC/BMCL/BMCU values to the gene sets in the
Reactome Pathways database (https://reactome.org/; last
accessed December 2022), the Kyoto Encyclopedia of Genes
and Genomes (KEGG; https://www.genome.jp/kegg/; last
accessed December 2022), or the GO database (http://gen-
eontology.org; last accessed December 2022). Gene sets that
contained at least 3 gene BMCs (genes that pass all criteria
in the analysis) and were at least 5% populated (based on

total annotated gene number) were selected. The tPOD
chosen represents the gene set from the pathway database
with the lowest median BMC. This approach was based on
one recommended by the National Toxicology Program
(NTP) (National Toxicology Program, 2018).

b) We set the 25th rank-ordered gene as a threshold to identify
the concentration at which a robust change in the tran-
scriptome has occurred (eg, a concerted molecular
response; Johnson et al., 2022). Chemicals that do not have
at least 25 genes fitting BMCs were considered “inactive”
based on analysis of this tPOD.

To derive an ERa-specific tPOD, we calculated the median
BMC of the 46 genes from the published ERa biomarker (Ryan
et al., 2016). Confidence Intervals for gene sets were estimated
using a parametric bootstrap. For each gene in the gene set, a
uniform distribution was assumed with the lower limit repre-
sented by the gene BMCL and the upper limit represented by the
gene BMCU. A bootstrapped sample consisted of randomly gener-
ated BMCs for each gene that had a BMC. The median value from
each of the 2000 bootstrap samples was used to estimate the
bootstrap distribution for the median. The 95% confidence inter-
val was estimated using the 2.5th and 97.5th percentiles from the
bootstrap distribution.

Pathway and upstream regulator enrichment analysis

Ingenuity Pathway Analysis (IPA; QIAGEN, Redwood City,
California) was used to identify perturbed upstream regulators,
canonical pathways, and disease & functions. For each chemical
tested, we imported Excel files into IPA containing gene IDs
(Ensembl and Gene Symbol) for the genes that fit a BMC (ie,
passed all filtering criteria), as well as their Williams trend test p
values from the BMC analysis, and the largest fold-changes of the
gene relative to solvent controls (exported data from
BMDExpress). This approach allowed us to identify enrichment of
genes showing robust concentration-responses to the exposures.
IPA Core Analysis with a gene expression threshold of absolute
fold change �1.5 and FDR-adjusted p value �.05 was used with
the direct and indirect relationship settings based on experimen-
tal and highly predicted data (focusing on human sources from
breast cancer cell lines). Statistical significance of the overlap
(FDR-adjusted p value �.05) between the data set and known tar-
gets of upstream regulators in IPA were calculated using Fisher’s
exact tests. The z-score was calculated using Fisher’s exact test
based on the expected relationship for directions between
upstream regulators and target genes and those observed in the
data set. A z-score of >2 (activated) or <2 (inhibited) was consid-
ered statistically significant.

Results
Cell viability assay
Cytotoxicity was evaluated using the CellTiter-Blue Cell Viability
assay. There was no evidence of cytotoxicity using this assay, as
there was no decline in fluorescence for BPA or any of the alter-
native chemicals at any of the concentrations tested following
48-h exposure (Supplementary Figure 1).

Cell proliferation assay
Cell proliferation was evaluated by counting cells following expo-
sure to chemicals for 7 days and comparing with respective sol-
vent controls. BPA (50, 100 mM) and several alternative chemicals
decreased cell counts below 50% of solvent control, including:

Matteo et al. | 269

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://reactome.org/
https://www.genome.jp/kegg/
http://geneontology.org
http://geneontology.org
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data


4,40-BPF (100 mM), BPAP (10, 50, 100 mM), BPS-MAE (50, 100 mM),
2,40-BPF (100 mM) which suggests cell cytotoxicity (Supplementary
Figure 2).

General HTTr data summary
Transcriptomic data were processed using the R-ODAF pipeline
that assesses quality of preprocessed data, implements study-
wide alignment QCs (Harrill et al., 2021), and produces fold
changes and p values for each gene with the DESeq2 R package
(Love et al., 2014). The median read depth was approximately 2
million per sample with a median range of 96.6% mapped reads.
Of 960 samples, 59 were lost to quality assurance (QA)/QC analy-
sis. All QC reports are available in Supplementary File 2. The data
were further filtered at the gene or probe level to produce a read
count matrix for analysis in BMDExpress v2.3 to derive lists of
concentration-responsive genes for enrichment analysis and
tPOD derivation for toxicological potency ranking.

Stress response biomarker activation
We used previously published transcriptomic biomarkers to
explore general stress responses and toxicities following expo-
sure to BPA and its substitutes (Supplementary File 3). The bio-
markers were specific to: (1) nuclear factor erythroid 2-related
factor 2 (Nrf2) activation/suppression (Rooney et al., 2020); (2)
DNA damage induction (TGx-DDI; Li et al., 2019); (3) metal respon-
sive transcription factor 1 (MTF1) activation/suppression (Jackson
et al., 2020); (4) nuclear factor kappa beta (NFjB) activation/sup-
pression (Korunes et al., 2022) and an unpublished biomarker
gene set (5) for activation of heat shock transcription factor 1
(HSF1; Cervantes and Corton, 2021). We also investigated an epi-
genotoxicity biomarker for histone deacetylase inhibition
(HDACi; Cho et al., 2021). Coactivation of these biomarkers within
a condition generally occurs when there are broad, non-specific
changes in transcription that result from overt cellular stress; ie,
the cytotoxic burst phenomenon that occurs at concentrations
close to cell death (Escher et al., 2020). Supplementary Table 1
provides a summary of all stress response biomarkers activated
and inhibited.

BPAP activated most of the stress response biomarkers at the
top 2 concentrations: Nrf2 (50, 100 mM); TGx-DDI (100 mM); MTF1
(50, 100 mM); NFkB (100 mM); and HSF1 (50, 100 mM). BPA activated
both Nrf2 and MTF1 at 100 mM. 2,40-BPF activated Nrf2 (50,
100 mM) and NFkB (50 mM). BPS activated NFkB (100 mM) and MTF
(0.001 mM). BPS-MAE activated Nrf2 and MTF1 at 100 mM. P201
activated MTF1 (0.01, 0.1, 1, 10, 100 mM) and NFkB (100 mM). A few
other chemicals activated a single biomarker: 4,40-BPF activated
the Nrf2 biomarker at 100 mM; 2,40-BPS activated NFkB at 100 mM;
and BADGE (0.1, 0.5, 1 mM) and BTUM (5 mM) activated MTF1. The
HDACi biomarker was inhibited by BPA (5–50 mM) and the follow-
ing chemicals: 4,40-BPF (5–100 mM), BPAF (1 mM), BPAP (5–50 mM),
BPC (5–10 mM), and BPS (5 mM).

Final sample inclusion and exclusion for ERa activity,
potency comparison, and enrichment analysis
The stress response biomarker analysis was used to identify
those concentrations of the chemicals that were causing overt
cellular stress. Specifically, a chemical concentration was consid-
ered cytotoxic (overly stressed) if it activated 2 or more bio-
markers or if it decreased cell counts below 50% of controls;
these concentrations were subsequently removed from further
analyses. In general, the highest exposure concentrations of BPA
(50, 100 mM) and several BPA alternative chemicals met these
criteria, including: 2,40-BPF (50, 100 mM), 4,40-BPF (100 mM) BPAP

(10, 50, 100 mM), BPS-MAE (50, 100 mM), and P201 (100 mM).
Supplementary Table 2 provides a summary of the final concen-
trations retained for: (1) ERa biomarker analysis; (2) BMC analysis;
(3) tPOD derivation; (4) UR and pathway analysis. For some chem-
icals like BPA and BPAP, stress response biomarkers were acti-
vated at highest concentrations, whereas ERa biomarker
activation at these concentrations decreased. These data support
that the cellular generalized stress response was associated with
a loss of ERa signaling. Further, this stress response is in line with
a cytotoxic “burst” phenomenon associated with exposure to
high concentrations of certain chemicals (Judson et al., 2016).

Identification of ERa receptor agonists and potency ranking
for ERa effects
The ability of each BPA alternative chemical tested to activate
the ERa was predicted using a published 46 gene ERa biomarker
(Ryan et al., 2016). Nine chemicals including BPA were predicted
to be ERa agonists (Figure 1 and Supplementary File 3).

The median BMCs of genes that fit models using BMDExpress
from the ERa biomarker were also used to rank order the potency
of these chemicals for ERa activity. The rank from most potent to
least potent based on this was: BPAF, BPC, BPA, 4,40-BPF, BPAP,
BPS, 2,40-BPF, P201, 2,40-BPS (Figure 1).

Potency ranking using BMC analysis
The total number of genes fitting BMC models is shown in
Figure 1. The top 2 most transcriptionally active by this analysis
(ie, chemicals with the most genes fitting BMC models) were BPS
and 4,40-BPF, with 1040 and 861 genes fitting models, respectively.
In contrast, 3 BPA alternative chemicals (DDS, BTUM, BPS-MPE)
had fewer than 25 genes fitting BMC models (presented in more
detail below).

To rank order BPA alternative chemicals by toxicological
potency, we calculated 2 general systemic tPODs to represent a
“tipping point” for induction of transcriptomic perturbations: (1)
the 25th rank-ordered gene BMC, and (2) the lowest pathway
gene set median BMC as per the U.S. NTP’s recommendation
(National Toxicology Program, 2018).

Gene accumulation plot shows the ranking of BPA and its
alternatives potency according to their 25th ranked-ordered gene
BMC (Figure 2). The order of potency from the most to least
potent (lowest BMC to highest BMC) was: BPAF, BPA, BPC, 4,40-
BPF, BPAP, BPS, 2,40-BPF, BADGE, P201, D-8, 2,40-BPS, TGSA, BPS-
MAE. Three BPA alternative chemicals (DDS, BTUM, BPS-MPE) did
not have at least 25 genes with modeled BMC and thus were con-
sidered “not transcriptionally active” at the concentrations tested
(Figure 1). General system tPODs were used to compare potency
with our prototype agent: BPA. The halogenated bisphenol BPAF
was more potent than BPA (BMC¼ 0.00098 mM and 0.0084 mM for
BPAF and BPA, respectively). BPA is more potent than BPC
(BMC¼ 0.033 mM) by one order of magnitude, followed by 5 alter-
native chemicals (4,40-BPF, BPAP, BPS, 2,40-BPF, BADGE) that had
BMCs greater than one order magnitude higher (0.15–0.61 mM)
compared with BPC. Importantly, the 25th rank-ordered gene
BMC of BPA is within one or order of magnitude of the mean
human urinary BPA concentration (Chen et al., 2016). The remain-
ing 5 (ie, least potent) transcriptionally active BPA alternative
chemicals (P201, D-8, 2,40-BPS, TGSA, BPS-MAE) had BMCs greater
than the previous 5 compounds by one order of magnitude (2.9–
7.4 mM).

We then ranked BPA and its alternatives by potency based on
the lowest median BMC derived from pathways (Figure 1). The
order of potency from most to least potent was: TGSA (GO), BPAF

270 | In Vitro Transcriptomic Analyses

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data


(GO), 2,40-BPF (Reactome), 2,40-BPS (KEGG), BPA (GO), BPC (GO),

4,40-BPF (Reactome), BPAP (GO), BPS (GO), P201 (Reactome), D-8

(GO), BPS-MAE (GO). A full list of gene sets for each chemical

tested is included in Supplementary File 4.

Pathway and upstream regulator enrichment analysis
We uploaded the lists of genes with BMCs, along with the maxi-

mum fold change and Williams trend test p values, into IPA to

determine if there was enrichment of genes fitting models to spe-

cific canonical pathways and/or diseases, or associated with reg-

ulation by specific upstream molecules. This approach allowed

us to compare chemicals across all concentrations tested and

revealed that the positive control (E2), BPA, and 10 BPA alterna-

tive chemicals (BPAF, BPC, 4,40-BPF, BPAP, BPS, BPS-MAE, 2,40-BPF,

2,40-BPS, BPC, P201, and D-8) perturbed gene sets associated with

predicted upstream regulators and canonical pathways

(Figure 3). A full list of upstream regulators, canonical pathways,

and disease and functions is in Supplementary File 5. The gene

set enrichment analyses were remarkably similar across the ERa

active chemicals and E2, which indicates a high degree of con-
cordance in the transcriptional alterations that they induce.

Upstream regulator analysis (Figure 3) provides additional
support to the weight of evidence that many of the BPA alter-
native chemicals are ERa activators (see Figure 3) operating
through similar mechanisms. As expected, the estrogen recep-
tor 1 (ERa) was predicted to be activated by E2 and all the ERa

agonists identified in the biomarker analysis. Further, E2 and
many of the BPA alternative chemicals were also predicted to
activate the retinoic acid receptor alpha (RARa). From this anal-
ysis, D-8 and BPS-MAE were predicted to activate ESR1 (ie, ERa)
and RARa, albeit weakly, although these chemicals were not
predicted to be ERa active using the biomarker. E3 ubiquitin
ligase constitutive photomorphogenic 1 (COP1) and Paf1/RNA
polymerase II complex component (CTR9) were predicted to be
activated only by BPA, BPAF, 4,40-BPF, BPAP, and BPS and ring
finger protein (RNF181) activation was predicted for BPAF, 4,40-
BPF, BPAP, and BPS. The 2 upstream regulators most inhibited
by E2 and many of the BPA alternative chemicals were: lysine
demethylase 5B (KDM5B) and the tumor protein 53 (TP53).

Figure 1. Comparison of estrogen receptor alpha (ERa) bioactivity and transcriptomic points of departure (tPODs) derived from the whole transcriptome
analysis in MCF-7 cells (n¼3–4 per concentration) exposed to bisphenol A and 15 alternative chemicals at a range of concentrations (0.0005–100 mM) for
48 h. ERa prediction is based on at least one concentration yielding an agonist or antagonist call. For tPOD derivation, data were prefiltered using the
Williams trend test (p� .05) and a fold-change of �1.5 or �1.5. Data were also postfiltered with the following settings in BMDExpress v2.3: Best BMD/
BMDL <20, Best BMDU/BMD <20, Best BMDU/BMDL <40, and Best fitPvalue >.1. tPODs representing the 25th rank ordered gene benchmark
concentration (BMC) (shown in mM), the median gene BMC for the lowest pathway (at least 3 genes and 5% of pathway) as well as the median gene BMC
for the ERa biomarker gene set are shown in the table and in the top panel. The chemicals are shown in decreasing order of potency based on tPODS
from the 25th gene BMC.

Matteo et al. | 271

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data
https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data


Overall, BPA, BPAF, 4,40-BPF, BPAP, and BPS perturbed the most
similar upstream regulators.

In general, the most affected canonical pathways were asso-
ciated with sustained proliferative signaling and downregula-
tion of growth suppressors (Figure 3). For example, E2 and most
ERa active chemicals tested were predicted to activate the
kinetochore metaphase signaling pathway and inhibit the cell
cycle and DNA damage checkpoint control regulation pathway.
Some perturbed pathways, including aryl hydrocarbon receptor
and salvage pathways of pyrimidines ribonucleotides, likely
represent a non-specific transcriptomic response to xenobiotic
exposure. Full interpretation of altered canonical pathways is
beyond the scope of this study but serves to group BPA alterna-
tive chemicals based on possible mode of action. In sum, most

ERa active BPA alternative chemicals, including BPA, BPAF, 4,40-
BPF, BPAP, BPS, and 2,40-BPF perturbed similar canonical path-
ways.

Discussion
A major goal of this study was to provide information on the
mode of action and potency of data-poor BPA alternatives rela-
tive to more well-characterized replacements. We used a previ-
ously published transcriptomic ERa biomarker (Ryan et al., 2016)
to identify ERa active compounds and derive a tPOD for potency
comparison. The biomarker analysis identified 9 agonists among
the 16 chemicals tested, in line with previous data. Indeed, sev-
eral BPA alternative chemicals have been identified as ERa active

Figure 2. Comparison of ERa bioactivity and tPODs derived from the whole transcriptome analysis in MCF-7 cells (n¼ 3–4 per concentration) exposed to
BPA and 15 alternative chemicals at a range of concentrations (0.0005–100 lM) for 48 h. ERa prediction is based on at least one concentration yielding
an agonist or antagonist call. For tPOD derivation, data were prefiltered using the Williams trend test (p� .05) and a fold-change of �1.5 or �1.5. Data
were also post-filtered with the following settings in BMDExpress v2.3: Best BMD/BMDL <20, Best BMDU/BMD <20, Best BMDU/BMDL <40, and Best
fitPvalue >.1. tPODs representing the 25th rank ordered gene BMC (shown in lM), the median gene BMC for the lowest pathway (at least 3 genes and 5%
of pathway) as well as the median gene BMC for the ERa biomarker gene set are shown in the table and in the top panel. The chemicals are shown in
decreasing order of potency based on tPODS from the 25th gene BMC. 25th gene, ERa biomarker, lowest median pathway.

Figure 3. Ingenuity Pathway Analysis of genes fitting benchmark concentration models in the MCF-7 cells (n¼ 3–4 per concentration) exposed to
bisphenol A and 15 alternative chemicals for 48 h, as well as positive control 17b-estradiol. A list of genes fitting models, the Williams trend test p value
and the maximum fold-change were imported into IPA for this comparison. Z-score and adjusted p value filters were set to �2.0 and �.05, respectively.
A, Upstream regulators. B, Top 10 canonical pathways.

272 | In Vitro Transcriptomic Analyses

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023



using other methods, including BPC, 4,40-BPF, BPS, BPAF, and
BPAP (see Keminer et al., 2020; Pelch et al., 2019). BPAF had the
lowest ERa median BMC, followed by BPC. BPA was the third most
estrogenic compound tested, and had comparable estrogenicity
to 4,40-BPF, BPAP, BPS, and 2,40-BPF. Importantly, our data align
with Mesnage et al. (2017) who also applied this biomarker in
MCF-7 cells to identify the following BPA replacements with ERa

activity: BPAF, 4,40-BPF, BPS, and BPAP. Moreover, BPAF was iden-
tified as the most potent chemical tested among 20 BPA alterna-
tives using ToxCast data for gene sets associated with nuclear
hormone receptors, including the ERa (Pelch et al., 2017). Our data
add to the literature base to support that BPAF is more estrogenic
than BPA and other alternative chemicals in vitro indicating that
more information is required to scrutinize BPAF as an alternative
to BPA. This information could serve to support read-across
chemical safety evaluations.

Several data-poor BPA replacements demonstrated weak or no
estrogenicity based on our biomarker analysis. Specifically, 2,40-
BPF was less estrogenic than its isomer (4,40-BPF) by one order of
magnitude. To our knowledge, we are the first group to report on
the estrogenicity of 2,40-BPF in MCF-7 cells. Previous work sug-
gests that P201 is non-estrogenic (Goldinger et al., 2015; Keminer
et al., 2020); our results suggest that this alternative chemical is
weakly estrogenic and ERa agonism occurs at high concentra-
tions (�50 mM). 2,40-BPS was previously predicted to be ERa inac-
tive (Keminer et al., 2020; Pelch et al., 2019); herein, we predicted
ERa activity at high concentrations (�50 mM). Other data-poor
BPA alternatives, including TGSA, D-8, BPS-MAE, and BPS-MPE
were predicted to be ERa inactive, in line with previous data
(Goldinger et al., 2015; Keminer et al., 2020; Kuruto-Niwa et al.,
2005; Pelch et al., 2019). We were unable to find other published
reports on the estrogenicity of DDS or BTUM. Our data indicate
that further investigations using established methods (eg, com-
petitive binding assays) are warranted to confirm the estrogenic
activity of 2,40-BPF and P201.

We explored 2 approaches to identify a “general toxicological
response tPOD” that reflects a concentration at which the tran-
scriptome is perturbed based on (1) the 25th ranked gene BMC or
(2) the lowest median pathway BMC. We previously used the 5th
percentile as a tPOD for deriving bioactivity exposure ratios to be
conservative (Rowan-Carroll et al., 2021); we note that percentile
gene BMCs are influenced by the top concentration and total
number of responsive genes. The 5th percentile can also be
derived from very few genes and we sought to identify a method
that ensures some minimal level of transcriptional change. Thus,
the 25th rank-ordered gene was used for the first time as a fixed
point that represents a robust transcriptomic change relative to
matched solvent controls (Reardon et al., 2021); the 25th gene is a
more stable point for direct comparisons across chemicals that is
less influenced by inclusion/exclusion of top concentrations rela-
tive to the 5th percentile BMC. This was an arbitrary but conser-
vative selection based on our previous analysis that such an
approach provides reproducible potency rankings for per- and
poly-fluoroalkyated chemicals and represents approximately
0.1% of the genes in the genome. The lowest median pathway
was included because it is also less influenced by inclusion/
exclusion of top concentrations and has previously been pro-
posed as an acceptable approach by a panel of experts (National
Toxicology Program, 2018).

The general toxicological response tPODs yielded potency
rankings that largely aligned with the ERa biomarker approach.
In particular, the 25th gene BMC yielded very similar potency
rankings as the ERa BMC, but was generally more conservative by

one order of magnitude. This suggests that there is biological
activity impacting genes not directly regulated by ERa at lower
exposure concentrations than those regulated by the receptor.
BPAF had the lowest 25th gene BMC, followed by BPA, BPC; all
other chemicals had 25th gene BMCs within 2 orders of magni-
tude of BPA. The lowest pathway median BMC was also very simi-
lar to the 25th gene BMC in most cases. Examination of the
accumulation plots indicates the 25th gene BMC serves as a suit-
able representation of the concentration at which transcriptional
activity ramps up. Interestingly, TGSA, 2,40-BPF, and 2,40-BPS had
very low pathway median BMCs, whereas they were among the
least potent chemicals tested according to the other 2 tPODs.
This may represent a limitation to the NTP approach, which is
not encountered using the 25th gene BMC. Further investigation
is warranted into the transcriptomic effects of data-poor BPA
alternative chemicals. Overall, these data align with previous
data from our group that found considerable overlap between
pathway- and global-based tPODs from microarray data of in vivo
toxicological dose-response studies (Farmahin et al., 2017).

We used a suite of previously published transcriptomic bio-
markers to explore stress responses induced by the chemicals in
the MCF-7 cells (Jackson et al., 2020; Korunes et al., 2022). Most
cellular stress response biomarkers, as well as the DNA damage
biomarker, were activated by BPA and alternative chemicals at
the highest (�50 mM) concentrations tested. We also used a novel
HSF1 biomarker to assess cellular stress and observed a similar
trend. Interestingly, BPA and several alternatives inhibited the
HDACi biomarker at non-cytotoxic concentrations. Thus, these
BPA alternatives may decrease histone acetylation, perhaps indi-
rectly, adding to a growing body of evidence that supports that
in vitro and in vivo BPA exposure affects this pathway (Qin et al.,
2020). Overall, those alternative chemicals with the most struc-
tural homology to BPA had the most similar patterns of bio-
marker perturbation.

Genes fitting BMCs revealed toxicodynamic similarities across
the chemical sets supporting read-across for data-poor BPA alter-
native chemicals. All chemicals identified as ERa active by bio-
marker analysis were also predicted to activate ERa as an
upstream regulator, along with RARa. RARa is necessary for ERE-
mediated transcription and proliferation in breast cancer cells
(Ross-Innes et al., 2010) and has significantly overlapping target
genes with ERa (Hua et al., 2009). Most ERa active BPA alternative
chemicals were predicted to inhibit KDM5B. Previous work sup-
ports a role for KDM5B in cultured breast cancer cell proliferation
(Zhang et al., 2019), and its predicted inhibition in our dataset
may represent a compensatory mechanism in response to sus-
tained proliferative signaling by BPA alternative chemicals.
Furthermore, a subset of the ERa active chemicals activated
upstream regulators that facilitate ERa signaling including CTR9
(Zeng et al., 2016), RNF31 (Zhu et al., 2014), RNF181 (Zhu et al.,
2020), and COP1 (Tang et al., 2022). These data support our identi-
fication of ERa active BPA alternatives and suggest that they
cause robust activation of the ERa proliferative network in MCF-7
cells.

Our pathway analysis suggests that several BPA alternatives
induce transcriptional changes in gene sets associated with non-
specific proliferative signaling and dysregulated cell cycle
control. MAPK1 and CREB1 were predicted as activated upstream
regulators by a subset of chemicals tested. Another group identi-
fied CREB1 as an upstream regulator in hippocampi of BPA-
exposed mice using RNA-seq (Henriksen et al., 2020). Further, the
most perturbed canonical pathways in this study by BPA and
alternative chemicals were the kinetochore metaphase signaling

Matteo et al. | 273

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023



pathway, mitotic roles of polo-like kinase, and G2/M DNA dam-

age checkpoint. Importantly, a previous transcriptomic profiling

study using RNA-seq in MCF-7 cells found that BPA and 5 alterna-

tive chemicals induced GO pathways associated with cell cycle

and breast cancer (Mesnage et al., 2017). Collectively, these

enriched gene set data suggest that exposure to BPA and several

alternative chemicals affect the ability of MCF-7 cells to regulate

cell division.
In conclusion, our HTTr data analysis pipeline enabled us to

identify potential ERa agonists within our 15 BPA alternative

chemicals, several of which are supported by additional evidence

in the literature. Application of BMC modeling facilitated toxico-

logical potency ranking based both on ER-specific hazard as well

as more generalized transcriptomic responses. A novel approach

to gene set enrichment based on genes fitting BMC models

revealed a high degree of similarity across the gene expression

profiles induced by these chemicals, supporting their grouping

for this and future analyses. This study also supports previous

data indicating that BPAF is a potent ERa agonist that induces

transcriptional changes at low concentrations, indicating that

further research of halogenated BPA alternative chemicals is war-

ranted. Given that these data are derived from an in vitro cell

model that overexpresses the ERa, they are limited in their ability

to generalize to human toxicity. Thus, further investigation in

primary cell models would be beneficial. Together, our study pro-

vides key support for future read-across applications and sug-

gests that several of the proposed BPA alternative chemicals may

cause hazards at similar concentrations to BPA.

Supplementary data
Supplementary data are available at Toxicological Sciences online.

Acknowledgments
We thank Eunnara Cho, Nikolai Chepelev, and Solomon Appovoo

for helpful comments on this manuscript.

Funding
This work was supported by Health Canada and with some sti-

pend support to G.M. through the Natural Sciences and

Engineering Research Council of Canada. The research was

undertaken, in part, thanks to funding to C.L.Y. from the Canada

Research Chairs Program.

Declaration of conflicting interests
The authors declared no potential conflicts of interest with

respect to the research, authorship, and/or publication of this

article.

References
Atlas, E., Bojanowski, K., Mehmi, I., and Lupu, R. (2003). A deletion

mutant of heregulin increases the sensitivity of breast cancer

cells to chemotherapy without promoting tumorigenicity.

Oncogene 22, 3441–3451.

Björnsdotter, M. K., de Boer, J., and Ballesteros-G�omez, A. (2017).

Bisphenol A and replacements in thermal paper: a review.

Chemosphere 182, 691–706.

Cervantes, P. W., and Corton, J. C. (2021). A gene expression bio-

marker predicts heat shock factor 1 activation in a gene expres-

sion compendium. Chem. Res. Toxicol. 34, 1721–1737.

Chen, D., Kannan, K., Tan, H., Zheng, Z., Feng, Y. L., Wu, Y., and

Widelka, M. (2016). Bisphenol analogues other than BPA: environ-

mental occurrence, human exposure, and toxicity – a review.

Environ. Sci. Technol. 50, 5438–5453.

Cho, E., Rowan-Carroll, A., Williams, A., Corton, J. C., Li, H. H.,

Fornace, A. J., Hobbs, C. A., and Yauk, C. L. (2021). Development

and validation of the TGx-HDACi transcriptomic biomarker to

detect histone deacetylase inhibitors in human TK6 cells. Arch.

Toxicol. 95, 1631–1645.

Cimmino, I., Fiory, F., Perruolo, G., Miele, C., Beguinot, F., Formisano,

P., and Oriente, F. (2020). Potential mechanisms of bisphenol A

(BPA) contributing to human disease. Int. J. Mol. Sci. 21, 5761–5722.

den Braver-Sewradj, S. P., van Spronsen, R., and Hessel, E. V. S.

(2020). Substitution of bisphenol A: a review of the carcinogenic-

ity, reproductive toxicity, and endocrine disruption potential of

alternative substances. Crit. Rev. Toxicol. 50, 128–147.

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S.,

Batut, P., Chaisson, M., and Gingeras, T. R. (2013). STAR: ultrafast

universal RNA-seq aligner. Bioinformatics 29, 15–21.

Escher, B. I., Henneberger, L., König, M., Schlichting, R., and Fischer,

F. C. (2020). Cytotoxicity burst? Differentiating specific from non-

specific effects in tox21 in vitro reporter gene assays. Environ.

Health Perspect. 128, 077007–077010.

Farmahin, R., Williams, A., Kuo, B., Chepelev, N. L., Thomas, R. S.,

Barton-Maclaren, T. S., Curran, I. H., Nong, A., Wade, M. G., and

Yauk, C. L. (2017). Recommended approaches in the application

of toxicogenomics to derive points of departure for chemical risk

assessment. Arch. Toxicol. 91, 2045–2065.

Gaidatzis, D., Lerch, A., Hahne, F., and Stadler, M. B. (2015). QuasR:

quantification and annotation of short reads in R. Bioinformatics

31, 1130–1132.

Goldinger, D. M., Demierre, A., Zoller, O., Rupp, H., Reinhard, H.,

Magnin, R., Becker, T. W., and Bourqui-Pittet, M. (2015).

Endocrine activity of alternatives to BPA found in thermal paper

in Switzerland. Regul. Toxicol. Pharmacol. 71, 453–462.

Harrill, J. A., Everett, L. J., Haggard, D. E., Sheffield, T., Bundy, J. L.,

Willis, C. M., Thomas, R. S., Shah, I., and Judson, R. S. (2021).

High-throughput transcriptomics platform for screening envi-

ronmental chemicals. Toxicol. Sci. 181, 68–89.

Henriksen, A. D., Andrade, A., Harris, E. P., Rissman, E. F., and

Wolstenholme, J. T. (2020). Bisphenol A exposure in utero

disrupts hypothalamic gene expression particularly genes sus-

pected in autism spectrum disorders and neuron and hormone

signaling. Int. J. Mol. Sci. 21, 3129–3116.

Hua, S., Kittler, R., and White, K. P. (2009). Genomic antagonism

between retinoic acid and estrogen signaling in breast cancer. Cell

137, 1259–1271.

Jackson, A. C., Liu, J., Vallanat, B., Jones, C., Nelms, M. D., Patlewicz, G.,

and Corton, J. C. (2020). Identification of novel activators of the metal

responsive transcription factor (MTF-1) using a gene expression bio-

marker in a microarray compendium. Metallomics 12, 1400–1415.

Johnson, K. J., Auerbach, S. S., Stevens, T., Barton-Maclaren, T. S.,

Costa, E., Currie, R. A., Dalmas Wilk, D., Haq, S., Rager, J. E.,

Reardon, A. J. F., et al. (2022). A transformative vision for an

omics-based regulatory chemical testing paradigm. Toxicol. Sci.

190, 127–132.

Judson, R., Houck, K., Martin, M., Richard, A. M., Knudsen, T. B.,

Shah, I., Little, S., Wambaugh, J., Woodrow Setzer, R., Kothiya, P.,

et al. (2016). Analysis of the effects of cell stress and cytotoxicity

274 | In Vitro Transcriptomic Analyses

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023

https://academic.oup.com/toxsci/article-lookup/doi/10.1093/toxsci/kfac127#supplementary-data


on in vitro assay activity across a diverse chemical and assay

space. Toxicol. Sci. 152, 323–339.

Keminer, O., Teigeler, M., Kohler, M., Wenzel, A., Arning, J., Kaßner,

F., Windshügel, B., and Eilebrecht, E. (2020). A tiered high-

throughput screening approach for evaluation of estrogen and

androgen receptor modulation by environmentally relevant

bisphenol A substitutes. Sci. Total Environ. 717, 134743.

Kitamura, S., Suzuki, T., Sanoh, S., Kohta, R., Jinno, N., Sugihara, K.,

Yoshihara, S., Fujimoto, N., Watanabe, H., and Ohta, S. (2005).

Comparative study of the endocrine-disrupting activity of

bisphenol A and 19 related compounds. Toxicol. Sci. 84, 249–259.

Korunes, K. L., Liu, J., Huang, R., Xia, M., Houck, K. A., and Corton, J.

C. (2022). A gene expression biomarker for predictive toxicology

to identify chemical modulators of NF-jB. PLoS One 17, e0261854.

Kupershmidt, I., Su, Q. J., Grewal, A., Sundaresh, S., Halperin, I., Alag,

S., Akhtari, S., and Ronaghi, M. (2010). Ontology-based meta-

analysis of global collections of high-throughput public data.

PLoS One 5, e13066.

Kuruto-Niwa, R., Nozawa, R., Miyakoshi, T., Shiozawa, T., and Terao,

Y. (2005). Estrogenic activity of alkylphenols , bisphenol S, and

their chlorinated derivatives using a GFP expression system.

Environ. Toxicol. Pharmacol. 19, 121–130.

Li, H. H., Yauk, C. L., Chen, R., Hyduke, D. R., Williams, A., Frötschl,

R., Ellinger-Ziegelbauer, H., Pettit, S., Aubrecht, J., and Fornace, A.

J. (2019). TGx-DDI, a transcriptomic biomarker for genotoxicity

hazard assessment of pharmaceuticals and environmental

chemicals. Front. Big Data 2, 36–14.

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of

fold change and dispersion for RNA-seq data with DESeq2.

Genome Biol. 15, 550–521.

Ma, Y., Liu, H., Wu, J., Yuan, L., Wang, Y., Du, X., Wang, R., Marwa, P.

W., Petlulu, P., Chen, X., et al. (2019). The adverse health effects of

bisphenol A and related toxicity mechanisms. Environ. Res. 176,

108575.

Mesnage, R., Phedonos, A., Arno, M., Balu, S., Corton, J. C., and

Antoniou, M. N. (2017). Transcriptome profiling reveals bisphenol

A alternatives activate estrogen receptor alpha in human breast

cancer cells. Toxicol. Sci. 158, 431–443.

National Toxicology Program. (2018). NTP Research Report on National

Toxicology Program Approach to Genomic Dose-Response Modeling:

Research Report 5 [Internet]. National Toxicology Program, Durham,

NC.

Pelch, K. E., Li, Y., Perera, L., Thayer, K. A., and Korach, K. S. (2019).

Characterization of estrogenic and androgenic activities for

bisphenol A-like chemicals (BPs): in vitro estrogen and androgen

receptors transcriptional activation, gene regulation, and binding

profiles. Toxicol. Sci. 172, 23–37.

Pelch, K. E., Wignall, J. A., Goldstone, A. E., Ross, P. K., Blain, R. B.,

Shapiro, A. J., Holmgren, S. D., Hsieh, J. H., Svoboda, D., Auerbach,

S. S., et al. (2017). NTP Research Report on Biological Activity of

Bisphenol A (BPA) Structural Analogues and Functional Alternatives:

Research Report 4. National Toxicology Program, Research

Triangle Park, NC.

Peshdary, V., Hobbs, C. A., Maynor, T., Shepard, K., Gagn�e, R.,

Williams, A., Kuo, B., Chepelev, N., Recio, L., Yauk, C., et al. (2021).

Transcriptomic pathway and benchmark dose analysis of

Bisphenol A, Bisphenol S, Bisphenol F, and 3,3’,5,5’-

Tetrabromobisphenol A in H9 human embryonic stem cells.

Toxicol. In Vitro 72, 105097.

Qin, T., Zhang, X., Guo, T., Yang, T., Gao, Y., Hao, W., and Xiao, X. F.

(2020). Epigenetic alteration shaped by the environmental chemi-

cal bisphenol A. Front. Genet. 11, 618966.

Reardon, A. J. F., Rowan-Carroll, A., Ferguson, S. S., Leingartner, K.,

Gagne, R., Kuo, B., Williams, A., Lorusso, L., Bourdon-Lacombe, J.

A., Carrier, R., et al. (2021). Potency ranking of per- and polyfluor-

oalkyl substances using high-throughput transcriptomic analysis

of human liver spheroids. Toxicol. Sci. 184, 154–169.

Rochester, J. R., and Bolden, A. L. (2015). Bisphenol S and F: a system-

atic review and comparison of the hormonal activity of bisphenol

A substitutes. Environ. Health Perspect. 123, 643–650.

Rogers, L. D. (2021). What does CLARITY-BPA mean for Canadians?

Int. J. Environ. Res. Public Health 18, 7001–7009.

Rooney, J. P., Chorley, B., Hiemstra, S., Wink, S., Wang, X., Bell, D. A.,

van de Water, B., and Corton, J. C. (2020). Mining a human tran-

scriptome database for chemical modulators of NRF2. PLoS One

15, e0239367.

Ross-Innes, C. S., Stark, R., Holmes, K. A., Schmidt, D., Spyrou, C.,

Russell, R., Massie, C. E., Vowler, S. L., Eldridge, M., and Carroll, J.

S. (2010). Cooperative interaction between retinoic acid receptor-

a and estrogen receptor in breast cancer. Genes Dev. 24, 171–182.

Rowan-Carroll, A., Reardon, A., Leingartner, K., Gagn�e, R., Williams,

A., Meier, M. J., Kuo, B., Bourdon-Lacombe, J., Moffat, I., Carrier,

R., et al. (2021). High-throughput transcriptomic analysis of

human primary hepatocyte spheroids exposed to per- and poly-

fluoroalkyl substances as a platform for relative potency charac-

terization. Toxicol. Sci. 181, 199–214.

Ryan, N., Chorley, B., Tice, R. R., Judson, R., and Corton, J. C. (2016).

Moving toward integrating gene expression profiling into high-

throughput testing: a gene expression biomarker accurately

predicts estrogen receptor a modulation in a microarray compen-

dium. Toxicol. Sci. 151, 88–103.

Stephens, M. (2017). False discovery rates: a new deal. Biostatistics 18,

275–294.

Tang, S. C., Lion, Q., Peulen, O., Chariot, P., Lavergne, A., Mayer, A.,

Fuster, P. A., Close, P., Klein, S., Florin, A., et al. (2022). The E3

ligase COP1 promotes ERa signaling and suppresses EMT in

breast cancer. Oncogene 41, 173–190.

Thomas, R. S., Bahadori, T., Buckley, T. J., Cowden, J., Deisenroth, C.,

Dionisio, K. L., Frithsen, J. B., Grulke, C. M., Gwinn, M. R., Harrill, J.

A., et al. (2019). The next generation blueprint of computational

toxicology at the U.S. Environmental protection agency. Toxicol.

Sci. 169, 317–332.

Verheijen, M. C. T., Meier, M. J., Asensio, J. O., Gant, T. W., Tong, W., Yauk,

C. L., and Caiment, F. (2022). R-ODAF: omics data analysis framework

for regulatory application. Regul. Toxicol. Pharmacol. 131, 105143.

Yeakley, J. M., Shepard, P. J., Goyena, D. E., Vansteenhouse, H. C.,

Mccomb, D., and Seligmann, B. E. (2017). A trichostatin A expres-

sion signature identified by TempO-Seq targeted whole transcrip-

tome profiling. PLoS One 12, e0178302.

Zeng, H., Lu, L., Chan, N. T., Horswill, M., Ahlquist, P., Zhong, X., and

Xu, W. (2016). Systematic identification of Ctr9 regulome in ERa-

positive breast cancer. BMC Genomics 17, 1–14.

Zhang, Z. G., Zhang, H. S., Sun, H. L., Liu, H. Y., Liu, M. Y., and Zhou,

Z. (2019). KDM5B promotes breast cancer cell proliferation and

migration via AMPK-mediated lipid metabolism reprogramming.

Exp. Cell Res. 379, 182–190.

Zhu, J., Li, X., Su, P., Xue, M., Zang, Y., and Ding, Y. (2020). The ubiqui-

tin ligase RNF181 stabilizes ERa and modulates breast cancer

progression. Oncogene 39, 6776–6788.

Zhu, J., Zhao, C., Kharman-Biz, A., Zhuang, T., Jonsson, P., Liang, N.,

Williams, C., Lin, C. Y., Qiao, Y., Zendehdel, K., et al. (2014). The

atypical ubiquitin ligase RNF31 stabilizes estrogen receptor a and

modulates estrogen-stimulated breast cancer cell proliferation.

Oncogene 33, 4340–4351.

Matteo et al. | 275

D
ow

nloaded from
 https://academ

ic.oup.com
/toxsci/article/191/2/266/6935810 by H

ealth C
anada user on 16 N

ovem
ber 2023