Introduction

Estrogens, including the endogenous hormone estradiol (E2), are key regulators for growth, differentiation, and a variety of physiological functions. E2 acts through the two estrogen receptors (ERs), $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$ (Lubahn et al. 1993; McDonnell and Norris 2002; Nilsson et al. 2001; Pettersson and Gustafsson 2001). These two receptors belong to the nuclear receptor (NR) superfamily of ligand-inducible transcription factors (TFs) (Hall and McDonnell 2005). Both exhibit distinct tissue-specific expression patterns and biological roles (Deroo and Korach 2006). Exposure to estrogen results in a dramatic increase in the amount of ER binding to chromatin (Taylor and Al-Azzawi 2000). When ERs bind to chromatin, their interaction with the RNA polymerase II holocomplex and the surrounding chromatin environment is modified by various coregulators (Lonard and O'Malley 2012). The major mechanisms of ER-mediated transcriptional gene regulation are a) the classical genomic mechanism, in which the ligand-bound ER interacts directly with DNA estrogen response elements (EREs) that are located in the regulatory regions of the ER target genes; b) the nonclassical genomic tethered mechanism, in which the ER interacts with other TFs, such as AP-1 or Sp1, to regulate gene expression; and c) the rapid nongenomic action mechanism, which activates intracellular signaling cascades involving p44/42 MAKP, Src, and Akt (Burns et al. 2011; Hall and McDonnell 2005; O'Lone et al. 2004).

In the past decade, more than 400 transcriptional coregulators of NRs have been identified and characterized (Dasgupta et al. 2014). Coregulators, which are either coactivators or corepressors, interact with NRs and other TFs to alter chromatin and facilitate gene activity (Dasgupta et al. 2014). Steroid receptor coactivators (SRCs), SRC-1, SRC-2, and SRC-3, were the first gene family to be identified and characterized as coactivators (Lonard and O'Malley B 2007). Each member of the SRC family can enhance the transcriptional activity of NRs by acting as a bridging molecule, which assists protein–protein interactions between NRs and multiple components that facilitate the assembly of the transcriptome complex at the target gene promoter (Dasgupta et al. 2014). Peroxisome proliferator-activated receptor gamma coactivator 1s (PGC1s), including $\text{PGC}1\mathrm{\alpha }$ and $\text{PGC}1\mathrm{\beta }$, are a family of coactivators that are key in energy metabolism (Knutti and Kralli 2001; Puigserver and Spiegelman 2003), but they also increase the transcriptional activity of many other NRs, including the ERs (Tcherepanova et al. 2000).

A large number of natural and synthetic chemicals are reported to disrupt normal functions of the endocrine system (Henley and Korach 2010). Bisphenol chemicals (BPs), including bisphenol A (BPA), bisphenol AF (BPAF), and bisphenol S (BPS), are a group of synthetic chemicals widely used in plastics, food packaging, and other products (Wetherill et al. 2007). BPA is an endocrine-disrupting chemical (EDC) that might be harmful to human health due to its ability to regulate the actions of ERs (Li et al. 2012; Rochester 2013; Welshons et al. 2006). Because of its use in a number of consumer products, concerns regarding the toxicological effects of BPA grew, and BPA is gradually being replaced with analogs such as BPAF and BPS (Liao et al. 2012a). BPAF is a fluorinated derivative of BPA used in many consumer products (Akahori et al. 2008). Humans are exposed to both BPA and BPAF, and effects on reproductive and developmental systems have been demonstrated in vitro (Bay et al. 2004; Bermudez et al. 2010). Recent studies have demonstrated that the estrogenic activity of BPS is similar to BPA (Grignard et al. 2012; Rochester and Bolden 2015). Thermal receipt paper is a potential source of occupational exposure to BPS (Thayer et al. 2015). Both BPA and BPS have been implicated in obesity and steatosis processes, but findings from a recent in vitro study suggest that they may act through different metabolic pathways (Héliès-Toussaint et al. 2014). However, overall, the molecular mechanisms attributed to the BPs are still poorly understood.

In this study, we performed functional analysis using several cell model systems to compare the agonistic effects of BPA, BPAF, and BPS on $\text{ER}\mathrm{\alpha }$ responses. Using Chromatin Immunoprecipitation (ChIP) analysis, we determined direct $\text{ER}\mathrm{\alpha }$ binding occurs in the presence of BPA, BPAF, and BPS to the ERE of an ER target gene promoter. We profiled the recruitment of coregulator-derived peptides in the ligand (BPA, BPAF, or BPS)–$\text{ER}\mathrm{\alpha }$ complexes using a peptide chip, Microarray Assay for Real-Time Coregulator-Nuclear receptor Interaction (MARCoNI) assay. Furthermore, we evaluated differences or similarities among E2, BPA, BPAF, and BPS in the $\text{ER}\mathrm{\alpha }$ ligand complex through computational analysis.

Materials and Methods

Chemicals

The 17 $\mathrm{\beta }\text{-}\text{estradiol}$ (E2) was purchased from Sigma-Aldrich, and ICI 182,780 (ICI) was purchased from Tocris Bioscience. BPA, BPAF, and BPS used in this study were provided by the Midwest Research Institute via a contract with the National Toxicology Program (NTP). All $10\text{-}\text{mM}$ stock solutions were made in DMSO and kept at $-20°\mathrm{C}$.

Plasmids

The expression vectors pcDNA3 was purchased from Invitrogen. An internal control plasmid for transfection efficiency, pRL-TK renilla luciferase (pRL-TK Luc), was purchased from Promega. The luciferase reporter plasmid, 3xERE Luc (synthetic vitellogenin ERE-TATA box sequence fused to a luciferase reporter gene), $\text{pcDNA}/\text{hER}\mathrm{\alpha }$, and $\text{pcDNA}/\text{hER}\mathrm{\beta }$ (full-length human ER expression plasmids) have been described previously. The following plasmids were gifts from D. McDonnell (Duke University): pcDNA/SRC-1, -2, -3, and $\text{pcDNA}/\text{PGC}1\mathrm{\alpha }$, $\mathrm{\beta }$.

Cell Lines and Tissue Culture

The human endometrial adenocarcinoma stable cell lines Ishikawa/vector (Ishikawa/vec) and $\text{Ishikawa}/\text{WT}\text{ER}\mathrm{\alpha }$ ($\text{Ishikawa}/\text{ER}\mathrm{\alpha }$) have been described previously (Burns et al. 2011, 2014). The human breast cancer cell line MCF-7 and human hepatocellular cancer cell line HepG2 were purchased from American Type Culture Collection. The human ovarian cancer cell line BG-1 FR has been described previously (Li et al. 2014). The human bone osteosarcoma epithelial recombinant U2OS cells stably expressing enhanced green fluorescent protein (EGFP)–tagged human $\text{ER}\mathrm{\alpha }$ fused to the C-terminus of EGFP were purchased from Thermo Fisher Scientific. Ishikawa/vec and $\text{Ishikawa}/\text{ER}\mathrm{\alpha }$ cells were maintained in phenol red–free Dulbecco’s Modified Essential Medium (DMEM):F12 medium (Invitrogen) supplemented with 10% Fetal Bovine Serum (FBS; Gemini Bio Products) and $1\mathrm{mg}/\mathrm{mL}$ Geneticin® (G418; Invitrogen). HepG2 cells were maintained in phenol red–free Minimum Essential Medium (MEM; Invitrogen) supplemented with 10% FBS and $4\text{mM}$ L-glutamine (Invitrogen). MCF-7 and BG-1 FR cells were maintained in phenol red–free DMEM:F12 medium (Invitrogen) supplemented with 10% FBS and $4\text{mM}$ L-glutamine. Recombinant $\mathrm{U}2\text{OS}/\text{ER}\mathrm{\alpha }$ cells were maintained in phenol red–free DMEM supplemented with 10% FBS and $4\text{mM}$ L-glutamine. For cell-starving conditions, 5% or 10% charcoal/dextran stripped FBS (sFBS; Gemini Bio Products) was substituted for FBS in the medium.

Transient Transfection and Luciferase Reporter Analysis

Cells were seeded in 24 well plates with 10% sFBS medium overnight. Effectene transfection reagent (Qiagen) was used for transient transfection. For ER-negative HepG2 cells, a total of $0.5\mathrm{\mu }\mathrm{g}$ of DNA, including $0.2\mathrm{\mu }\mathrm{g}$ of $\text{ER}\mathrm{\alpha }$ or $\text{ER}\mathrm{\beta }$, $0.2\mathrm{\mu }\mathrm{g}$ of 3xERE Luc, and $0.1\mathrm{\mu }\mathrm{g}$ of pRL-TK Luc plasmids were transfected into the cells. For Ishikawa/vec, $\text{Ishikawa}/\text{ER}\mathrm{\alpha }$ stable cells, and ER-positive cells (MCF-7 and BG-1 FR), a total of $0.3\mathrm{\mu }\mathrm{g}$ of DNA, including $0.2\mathrm{\mu }\mathrm{g}$ of 3xERE Luc and $0.1\mathrm{\mu }\mathrm{g}$ of pRL-TK Luc plasmids, were transfected. After 8 h, cells were changed to fresh 10% sFBS medium overnight and then treated with vehicle control (DMSO, final concentration $<0.01\%$), E2, or BPAs (1, 10, 100, and $1,000\text{nM}$). For coregulator experiments in HepG2 cells, a total of $0.8\mathrm{\mu }\mathrm{g}$ of DNA, including $0.1\mathrm{\mu }\mathrm{g}$ of $\text{ER}\mathrm{\alpha }$, $0.4\mathrm{\mu }\mathrm{g}$ of either SRC-1, SRC-2, SRC-3, $\text{PGC}1\mathrm{\alpha }$, or $\text{PGC}1\mathrm{\beta }$, $0.2\mathrm{\mu }\mathrm{g}$ of 3xERE Luc, and $0.1\mathrm{\mu }\mathrm{g}$ of pRL-TK plasmids were transfected into the cells. After 8 h, cells were changed to fresh 10% sFBS medium overnight and then treated with vehicle control (DMSO, final concentration $<0.01\%$), E2 ($10\text{nM}$), or BPAs (100 and $1,000\text{nM}$) for 18 h.

Luciferase assays were performed using the Dual Luciferase Reporter Activity System (Promega). Transfection efficiency was normalized against the renilla luciferase. Fold changes were calculated relative to vehicle controls. All experiments were repeated at least three times. Data shown are the average of triplicate determinations in a representative experiment. Values were calculated relative to vehicle control and presented as $\pm \text{standard}\text{error}\text{of}\text{the}\text{mean}$ (SEM).

RNA Extraction and Real-Time Polymerase Chain Reaction Analysis

Cells were seeded in $60\times 15\mathrm{mm}$ dishes with 5% sFBS medium and starved for 2 d. After changing to fresh 5% sFBS medium, cells were treated with vehicle control (DMSO, final concentration $<0.01\%$), $10\text{nM}$ E2, or $100\text{nM}$ BPs for 18 h. Total RNA was extracted using RNeasy Mini Kit (Qiagen). First-strand cDNA synthesis was performed using SuperScript reverse transcriptase according to the manufacturer’s protocol (Invitrogen). The mRNA levels were measured using SYBR™ green assays (Applied Biosystems). The sequences of primers used in quantitative real-time polymerase chain reaction (qPCR) and ChIP qPCR are shown in Table S1. Cycle threshold values were obtained using the ABI PRISM™ 7900 Sequence Detection System and analysis software (Applied Biosystems). Each sample was normalized against $\mathrm{\beta }\text{-}\text{actin}$ expression. Experiments were repeated three times, and results are presented as $\text{mean}\pm \text{SEM}$. Changes of endogenous gene expression by E2 were calculated as fold change relative to the vehicle control group in each cell line.

Chromatin Preparation and ChIP–Quantitative PCR Analysis

Cells were seeded in $100\times 20\mathrm{mm}$ dish with 5% sFBS medium and starved for 2 d. After changing to fresh 5% sFBS medium, cells were treated with vehicle control (DMSO, final concentration $<0.01\%$), E2 ($10\text{nM}$), or BPs ($100\text{nM}$) for 40 min. Chromatin preparation and ChIP assay were carried out with $\text{anti}\text{-}\text{ER}\mathrm{\alpha }$, anti-p300, or anti-CBP antibody using the Chromatin Immunoprecipitation Assay Kit (Cat.# 17-295; Millipore/Upstate®). Purified chromatin DNA was used for qPCR analysis. The sequences of the ChIP primers used in real-time PCR were as follows: for amplicon $-699/-422$ ($\text{AP}\text{-}1+\text{ERE}$): the forward primer 5′-TGGCACATAGTAGGTCCAGC-3′, reverse primer 5′-GTTTCTGGCAGGCAGATTGG-3′; for amplicon $\mathrm{-593}/-422$ (ERE): the forward primer 5′-TGTTGTGCCTCCAGCTCCTG-3′, reverse primer 5′-GGTTTCTGGCAGGCAGATT-3′. Relative recruitment was normalized to signals obtained from input DNA. Experiments were repeated three times, and results are presented as $\text{mean}\pm \text{SEM}$.

Cell Lysis Preparation and Microarray Assay for Real-Time Coregulator-Nuclear Receptor Interaction Assay

Recombinant U2OS cells stably expressing human $\text{ER}\mathrm{\alpha }$ fused to the C-terminus of EGFP were maintained in phenol red–free DMEM in the presence of $2\text{mM}$ L-Glutamine, 1% penicillin–streptomycin, $0.5\mathrm{mg}/\mathrm{mL}$ G418, and 10% FBS at 37°C, 5% ${\text{CO}}_2$, 95% humidity incubator. Snap-frozen cell pellets were slowly thawed on ice and suspended in ice-cold lysis buffer (Perera et al. 2017) supplemented with Halt™ protease inhibitor cocktail (Pierce Biotechnology), using 150 uL per $1\times {10}^7$ cells. The cell suspension was further snap-frozen, lysed, and centrifuged to obtain clear supernatant for the analyses later. Cell lysates were isolated from the full-length human $\text{ER}\mathrm{\alpha }\text{-}\text{transfected}$ U2OS cells. For MARCoNI, $25\mathrm{\mu }\mathrm{L}$ assay mixtures that contain cell lysates, $25\text{nM}$ of Alexa Fluor 488–conjugated GFP antibody (Invitrogen), and $0.2\text{mM}$ ligand (E2, BPA, BPAF or BPS, prediluted in DMSO, final concentration, 2%) were prepared on ice. The ligand-modulated coregulator interactions with the EGFP-tagged $\text{ER}\mathrm{\alpha }$ was assessed using a PamChip® plate that contained the 154 coregulator-derived binding peptides (Excel Table S1), including the LXXLL coactivator motif or the LXXXIXXXL corepressor motif from 66 different coregulators (PamGene International B.V.) (Beekmann et al. 2015; Koppen et al. 2009). The reaction was performed in a fully automated microarray processing platform at 20°C for 80 cycles (two cycles per minute). After removal of the unbound receptor by washing each array with $25\mathrm{\mu }\mathrm{L}$ in Tris-buffered Saline (TBS) buffer, images were obtained by the Charge-couple device (CCD) camera of the PamStation®96, and then were analyzed for quantification of ER binding using BioNavigator software (PamGene International B.V.). For data and statistical analysis, BioNavigator, Microsoft Excel (version 14.0.7,106.5,003; Microsoft Corporation) IBM SPSS Statistics (version 22; International Business Machines Corporation), and R (version 2.15.2; R Development Core Team), a language and environment for statistical computing, were used. Modulation index (MI) is the ${\text{log}}_{10}$-transformed relative binding value, which is calculated by the compound’s binding value, relative to the vehicle control (DMSO) binding value. Positive interaction indicates that ligands increase $\text{ER}\mathrm{\alpha }–\text{coregulator}$ motif interactions, and the cutoff is MI, $>0$; the binding value, $>50$; the relative binding value, $>1$; p-value, $<0.05$. Negative interaction indicates that ligands decrease $\text{ER}\mathrm{\alpha }\text{-}\text{coregulator}$ motifs interactions and the cutoff is MI, $<0$; the binding value, $>50$; the relative binding value, $<1$; p-value, $<0.05$.

Molecular Modeling

Using the X-ray crystal structures of E2- (pdb ID 3UUD), BPA- (pdb ID 3UU7), and BPAF- (pdb ID 3UUA) bound to ER, starting structures of the ligand-binding domains of ER were created for molecular dynamics (MDs). Missing residues (462–464 in 3UUD, 462–466 in 3UU7, and 462–463 in 3UUA) were introduced using the program Modeller (version 9.14; Andrej Sali) (Sali and Blundell 1993) and missing side chains (305, 306, 397, 467–477, 492, 531, and 545 in 3UUD; 306, 330, 419, 467–472, 531–533, and 548 in 3UU7; and 306, 330, 334, 335, 337, 339, 342, 363, 413, 416, 417, 419, 423, 466–473, 477, 492, 523, 525, and 533 in 3UUA) were reconstructed using the program Coot (version 0.8.2; Paul Emsley) (Emsley et al. 2010). Initial configurations of E2–ER were directly constructed from the crystal structure 3UUD. Initial structures with the normal pose of the ligands BPA, BPAF, and BPS bound in agonist conformations were modeled using the crystal structure of 3UU7, and the alternate poses of BP compounds were based on the crystal structure of 3UUA. Starting conformations of the ligand-free ER came from 3UUD (in which ER was in the agonist conformation) and the antagonist conformation of ER found in BPC bound to ER (pdb ID 3UUC) after removing the ligands. Protons were added, followed by the addition of counter ions, and the resulting complex structures were then solvated in a box of water for each system. The amount of water molecules varied from 25,500 to 26,800 for these systems. Prior to equilibration, each system was subjected to the following steps: a) over a nanosecond of constrained MDs with the peptide and the ligand constrained to the original position with a force constant of $10\text{kcal}/\text{mol}/\text{nm}$, b) minimization, c) low-temperature, constant-pressure MD simulation to assure a reasonable starting density, d) minimization, e) stepwise heating MD at constant volume, and f) a constant-volume MD run for another 9 nanoseconds. All final unconstrained trajectories were calculated at 300 K under constant volume MD (300 ns total dynamics for ligand bound systems and 400 ns for ligand-free systems at 1 fs time step) using PMEMD from the Amber 14 program to accommodate long-range interactions (Case et al. 2014). The amino acid parameters were taken from the Amber SB14 force field. The charges of atomic positions of the ligands were derived using the ChelpG scheme with the 6–31g(d) basis, set at the B3LYP level using the program Gaussian 09.D01 (Gaussian Inc.). The other ligand force field parameters were generic parameters listed in the Amber force field. Additional variable length MD simulations were performed with the same protocol for each ligand in water to establish the solvation behavior and energies of each ligand. Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) calculations were performed for all the ligand-bound systems. Configurations saved from the final 250 ns of the dynamic trajectories were used for the MM/GBSA calculations. The Generalized Born model (Onufriev et al. 2004) ($\text{default parameter}\text{igb}=5$) was used to calculate polar contributions of desolation. The dielectric constants of 80 and 1 were used for the solvent and solutes. Using the normal mode analysis, conformational entropy contributions to the binding energy were calculated using Amber 14.

Statistical Analysis

One-way and two-way analysis of variance with Dunnett’s Multiple Comparison test (**$p<0.01$, ***$p<0.001$ or ****$p<0.0001$) were performed using GraphPad Prism (version 7.01; GraphPad Software, Inc.).

Results

Estrogenic Activation of $\mathit{\text{ER}}\mathrm{\alpha }$ and $\mathit{\text{ER}}\mathrm{\beta }$ by BPA, BPAF, and BPS, Alterations in Expression of ER-Target Genes, and Direct Binding of $\mathit{\text{ER}}\mathrm{\alpha }$ to an ERE Located in the Promoter of ER-Regulated Genes

To evaluate the ERE-mediated transcriptional activity of $\text{ER}\mathrm{\alpha }$ when treated with BPA, BPAF, and BPS, promoter activation was performed using several cell systems. The chemical structures of compounds used in this study are shown in Figure 1A. The ER dependence of estrogenic activation was confirmed using human endometrial cancer Ishikawa cells stably expressing $\text{ER}\mathrm{\alpha }$ cells (Burns et al. 2011, 2014). As shown in Figure 1B, $1\text{nM}$ E2 induced ERE-mediated activation by fourfold. BPA-induced ERE-mediated activation by threefold at $100\text{nM}$. BPAF induced ERE-mediated activation by threefold at $10\text{nM}$, and matched the efficacy of E2 at $1,000\text{nM}$. BPS only induced ERE-mediated activation by threefold at the highest concentration ($1,000\text{nM}$) in $\text{Ishikawa}/\text{ER}\mathrm{\alpha }$ cells. ERE-mediated activation was weak or absent in Ishikawa/vector cells, suggesting that the activation by BPA, BPAF, and BPS are predominantly $\text{ER}\mathrm{\alpha }\text{-}\text{dependent}$. To investigate the estrogenic activity in the different cell types through an endogenous receptor, human breast cancer MCF-7 ($\text{ER}\mathrm{\alpha }$ positive) and ovarian cancer BG-1 FR ($\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$ positive) cells (Burns et al. 2014; Li et al. 2014) were used in reporter assays. In MCF-7 cells, $1\text{nM}$ E2 induced ERE-mediated activation by fourfold. However, BPA, BPAF, and BPS induced the reporter activation weakly, even at the higher concentrations ($1,000\text{nM}$) (Figure 1C, left). In BG-1 FR cells, similar dose responses were observed after treatment with E2 or BPA and the two analogs ($1–1,000\text{nM}$, Figure 1C, right). These data indicate all transcriptional activations with BPA, BPAF, and BPS occur in a concentration-dependent manner, and the sensitivity of responses to the individual compounds are cell type–specific.

To compare the $\text{ER}\mathrm{\alpha }$- and $\text{ER}\mathrm{\beta }\text{-}\text{mediated}$ transcriptional activation, we performed reporter assays in human hepatocellular cancer HepG2 (ER negative) cells after transfecting with an $\text{ER}\mathrm{\alpha }$ or $\text{ER}\mathrm{\beta }$ expression plasmid. For $\text{ER}\mathrm{\alpha }$ ERE-mediated activation, cells were highly responsive to E2, with up to 35-fold increases at $1\text{nM}$ (Figure 1D, left). BPA induced this activation by 25-fold at $100\text{nM}$. BPAF induced activation by 17-fold at $10\text{nM}$ and 30-fold at $100\text{nM}$. In contrast, BPS only weakly activated the $\text{ER}\mathrm{\alpha }$ ERE at $100\text{nM}$ (eightfold). With the highest concentration of $1,000\text{nM}$, all three chemicals displayed strong activation of $\text{ER}\mathrm{\alpha }$ by greater than 30-fold. For $\text{ER}\mathrm{\beta }$, E2 induced ERE-mediated activation with a lower fivefold increase at $1\text{nM}$ (Figure 1D, right). BPA at $100\text{nM}$ and BPAF at $10\text{nM}$ had weaker activation when compared to E2 at the same concentrations. BPA at $1,000\text{nM}$ and BPAF at $100–1,000\text{nM}$ had strong activity for $\text{ER}\mathrm{\beta }$ ERE-mediated activation. However, BPS did not have $\text{ER}\mathrm{\beta }$ ERE-mediated activation in HepG2 cells. These results indicate that BPA and BPAF can activate ERE-mediated transcription via $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$, but BPS exhibited more ligand specificity for $\text{ER}\mathrm{\alpha }$.

To investigate the ER-dependent response of BPA, BPAF, and BPS, we examined their effects on three well characterized ER target genes, GREB1 (gene regulation by estrogen in breast cancer 1), PGR (progesterone receptor), and WISP-2/CNN5 (WNT1-inducible signaling pathway protein 2) using real-time PCR analysis in MCF-7 cells (Figure 1E). All three chemicals, with varying degrees, significantly induced GREB1, PGR, and WISP-2/CNN5. In addition, the pure ER antagonist, ICI, blocked expression of all three ER target genes when cotreated with E2, BPA, BPAF, and BPS, indicating the ER receptor specificity of the activity response. These results indicate that BPA, BPAF, and BPS can alter transcriptional regulation directly on endogenous ER target genes.

To determine if $\text{ER}\mathrm{\alpha }$ directly interacts with ERE binding sites in the presence of BPA, BPAF, and BPS, we performed ChIP–qPCR analysis in MCF-7 cells. Promoter analysis showed that the upstream region of WISP-2/CNN5 harbors two AP-1 sites ($-624/-612$ and $-654/-641$) and an ERE site ($-562/-548$) (Figure 1F, top). As expected, E2 increased $\text{ER}\mathrm{\alpha }$ enrichment significantly to both regions by 14-fold on the $-699/-422$ ($\text{AP}\text{-}1+\text{ERE}$) and fivefold on the $-593/-422$ (ERE) regions (Figure 1F, bottom). BPA and BPAF significantly increased $\text{ER}\mathrm{\alpha }$ enrichments in both $-699/-422$ ($\text{AP}\text{-}1+\text{ERE}$) and $-593/-422$ (ERE) regions. BPS weakly increased $\text{ER}\mathrm{\alpha }$ enrichments in the $-593/-422$ (ERE) region (Figure 1F, bottom). However, $\text{ER}\mathrm{\alpha }$ enrichment with BPS was not seen in the $-699/-422$ ($\text{AP}\text{-}1+\text{ERE}$) region of the WISP-2 gene promoter. These data demonstrate $\text{ER}\mathrm{\alpha }$ directly binds differentially to EREs of ER target gene promoters when treated with BPA or BPAF.

$\mathit{\text{ER}}\mathrm{\alpha }–\mathit{\text{Coregulator}}$ Interaction Profiling and Comparison of BPA, BPAF, and BPS Coregulator-Derived Binding Peptides on $\mathit{\text{ER}}\mathrm{\alpha }$

To further support the molecular basis of BPA actions, we investigated ligand-modulated coregulator recruitment using MARCoNI. A working model of MARCoNI is shown in Figure S1. The ligand (E2, BPA, BPAF, and BPS)-modulated interaction of $\text{ER}\mathrm{\alpha }$ with coregulators was quantified (Figure 2). In the heatmap, the red color depicts positive interactions, which suggests that the ligand increases peptide binding, and the blue color depicts negative interactions, which suggests that the ligand decreases peptide binding to $\text{ER}\mathrm{\alpha }$ (Figure 2A). As a parameter for modulation of the interaction of $\text{ER}\mathrm{\alpha }$ with each coregulator peptide motif on the array, MI was calculated (Figure 2B). The complete data set of the peptide array analysis for E2, BPA, BPAF, and BPS is summarized in Excel Table S2.

A summary of the positive interaction peptides in the $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes are shown in Figure 2C. E2 facilitated the binding of 74 peptides from 32 different coregulators. In contrast, BPA significantly enhanced the binding of 8 peptides from 7 different coregulators. BPAF significantly enhanced the binding of 9 peptides from 9 different coregulators, and BPS significantly enhanced the binding of 25 peptides from 14 different coregulators (Figure 2C, left). To further understand the differential effects of BPA, BPAF, and BPS on $\text{ER}\mathrm{\alpha }$ function, the recruitment of coregulator motifs, as induced by BPA, BPAF, and BPS, were compared to that observed with E2 stimulation of the $\text{ER}\mathrm{\alpha }$. BPA- and BPS-bound peptides overlapped 100% with E2-modulated peptides. In comparison, only 33% of BPAF-bound peptides (3 out of 9) were also recruited by E2 (Figure 2C, right). Comparing the profiles of positive interactions of BPA, BPAF, and BPS to E2 on the $\text{ER}\mathrm{\alpha }$ complex show highly differential interactions (Tables S2–S5). In certain instances, the interactions were similar, but with a lower binding value compared to E2 (e.g., NCOA1_620_643 was 2.0 for BPA; 45.1 for E2) and with another portion of the same coactivator, NCOA1, no binding occurred (NCOA1_1421_1441) (Tables S2–S5). Analysis of the BPAF profile showed overall low binding values, but showed interactions to different factors not bound by E2 [e.g., peroxisome proliferator-activated receptor gamma, coactivator-related 1 (PPRC1)] (Tables S2–S5). The $\text{BPS}–\text{ER}\mathrm{\alpha }$ complexes exhibited more interactions than with BPA or BPAF, but these binding activities were generally lower than with E2 in the $\text{ER}\mathrm{\alpha }$ complexes, and in certain examples, selective to the BPS complex only (e.g., JHD2C, LCOR, and certain portions of MED1_591_614) (Tables S2–S5). In general, BPS had weaker interactions with the SRC family members than BPA or BPAF. The potential coregulators from positive interaction peptides in the ligand–$\text{ER}\mathrm{\alpha }$ complexes are summarized in Figure 2C, bottom. Interestingly, the unique peptides that were derived from PR, TF65, MAPE, KIF11, and PRDP2 were recruited only by BPAF.

A summary of the negative interaction peptides in the $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes are shown in Figure 2D. E2 was found to significantly modulate $\text{ER}\mathrm{\alpha }$ binding with 35 peptides from 22 different coregulators, BPA bound 44 peptides (28 coregulators), BPAF bound 75 peptides (38 coregulators), and BPS bound 31 peptides (23 coregulators) in the $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes (Figure 2CD, right). Comparing the profiles of the negative interactions showed that all three BPs overlapped with the E2 profile, but to varying degrees (Tables S6–S9). Fifty-seven percent of BPA-bound peptides (25 out of 44) overlapped with E2-modulated peptides, and 77% of BPS peptides (24 out of 31) overlapped with E2 peptides. However, only 32% of BPAF peptides (24 out of 75) recruited were the same as E2 (Figure 2D, right). Interestingly, all three BPs showed unique peptide interactions with factors not bound by E2 (e.g., BLIS1, CBP, GNAQ) (Tables S6–S9). With NCOR1 and NCOR2, there were interactions with E2 (NCOR1_2039_2061 and NCOR2_2330_2352), but at different sites in the coactivator (Tables S6–S9). NELFB interacted with E2, BPAF, and BPS, but not with the $\text{BPA}–\text{ER}\mathrm{\alpha }$ bound complex (Tables S6–S9). The potential coregulators (from negative interaction peptides) in the $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes are summarized in Figure 2D, bottom. Interestingly, the peptides that were derived from N-CoR1, N-CoR2/SMART, OIP4, and TIF1-a were recruited by all ligands. Thus, these results together demonstrate some of the first evidence that the coregulators in the $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes are differentially recruited in the presence of BPA, BPAF, or BPS.

Differential Recruitments of Coactivators SRCs and PGCs in $\mathit{\text{Ligand}}–\mathit{\text{ER}}\mathrm{\alpha }$ Complexes

To evaluate the peptide screening results, we selected two well-characterized coactivator families, SRC and PGC1, for further experimentation. To determine the role of $\text{ER}\mathrm{\alpha }$ with coregulators in ERE-mediated transcriptional regulation using reporter assays, we co-transfected SRCs or PGC1s with $\text{ER}\mathrm{\alpha }$ in HepG2 cells. As expected, SRC-1, SRC-2, SRC-3, $\text{PGC}1\mathrm{\alpha }$, and $\text{PGC}1\mathrm{\beta }$ induced $\text{ER}\mathrm{\alpha }$ ERE-mediated activity when treated with $10\text{nM}$ E2 (Figure 3A). For SRCs, in the presence of 100 or $1,000\text{nM}$ BPA, SRC-1 significantly induced $\text{ER}\mathrm{\alpha }$ ERE-mediated activity at both concentrations, but SRC-2 had minimal effects at the lower concentration, and SRC-3 showed no effect at either concentration (Figure 3B). In contrast, BPAF at 100 or $1,000\text{nM}$ had further stimulation with SRC-3, but no additional effect was observed with SRC-1 or SRC-2 at either concentration. In the presence of $1,000\text{nM}$ BPS, all three SRCs significantly activated $\text{ER}\mathrm{\alpha }$ ERE-mediated transcription. However, there was no coactivation with either SRCs at $100\text{nM}$. $\text{PGC}1\mathrm{\alpha }$ significantly induced $\text{ER}\mathrm{\alpha }$ ERE-mediated activity at high concentrations of $1,000\text{nM}$ BPA, and $\text{PGC}1\mathrm{\beta }$ was active by BPA at the lower concentration with $100\text{nM}$ (Figure 3C). $\text{PGC}1\mathrm{\alpha }$ also significantly activated the $\text{ER}\mathrm{\alpha }$ ERE-mediated transcription in the presence of 100 or $1,000\text{nM}$ BPAF. $\text{PGC}1\mathrm{\beta }$, however, further induced activity only with $1,000\text{nM}$ BPAF. BPS showed a clear differential preference of activation with $\text{PGC}1\mathrm{\alpha }$ in the presence of $1,000\text{nM}$ BPS, but $\text{PGC}1\mathrm{\beta }$ did not show any additional coactivator activity with BPS. These results confirmed the coactivator peptide profile experimentally detected by the MARCoNI assay as shown in Table 1. These data also demonstrate that BPA, BPAF, and BPS induce differential recruitments and interactions with coactivators such as SRCs and PGC1 in the different $\text{ligand}–\text{ER}\mathrm{\alpha }$ complexes.

Simulating BPA-, BPAF-, and BPS-Bound $\mathit{\text{ER}}\mathrm{\alpha }$ by Molecular Dynamics

To investigate if any potential modifications to structure appear on the coregulator binding surface of $\text{ER}\mathrm{\alpha }–\text{LBD}$ in the presence of BPs at its natural ligand-binding site, we performed a series of MD simulations. In the X-ray crystal structure determined by Delfosse et al. in 2012, two distinct ligand conformations were identified for bisphenol compounds when bound to $\text{ER}\mathrm{\alpha }$ at the ligand-binding site (Delfosse et al. 2012). In that work, BPA was bound to $\text{ER}\mathrm{\alpha }$ in a conformation like the natural substrate, E2 (Figure S2A), and we refer to this mode of the ligand as the canonical conformation, or conformation 1. In this canonical mode of the ligand binding, one phenolic group of BPA interacted with H524 and M421, while the other phenolic oxygen interacted with E353 and R394 and helix 12 (h12) was in a similar conformation to $\mathrm{E}2–\text{ER}\mathrm{\alpha }$ (Figure S2B). In addition to the canonical conformation (Figure S2C), for $\text{BPAF}–\text{ER}\mathrm{\alpha }$, the ligand was found to be bound in an alternate conformation (Figure S2D), or conformation 2. Interactions of BPAF with ER residues in the canonical conformation resembled those of BPA, and the alternate conformation displayed a phenolic group interacting with T347 instead of H524 and M421 (Figure S2D). In the X-ray crystal structures of $\text{BPAF}–\text{ER}\mathrm{\alpha }$, even though BPAF conformations were different, h12 was observed in a similar conformation to BPA– or $\mathrm{E}2–\text{ER}\mathrm{\alpha }$ (i.e., $\text{ER}\mathrm{\alpha }$ was in the agonist form). Meanwhile, when another BPA analog, bisphenol C (BPC) was used in the crystal trials (Delfosse et al. 2012), only the alternate conformation of the ligand was observed in the X-ray structure of $\text{BPC}–\text{ER}\mathrm{\alpha }$. Interestingly, in the case of $\text{BPC}–\text{ER}\mathrm{\alpha }$, h12 was located elsewhere with a conformation labeled as an antagonist form of $\text{ER}\mathrm{\alpha }$. It is worth noting that the $\text{ER}\mathrm{\alpha }$ used in the X-ray crystallographic study (Delfosse et al. 2012) was a mutant (Y537S) that enabled successful crystal trials. In contrast, we used the published X-ray crystal coordinates as our starting structures in the MD simulations with Y537 (wide-type $\text{ER}\mathrm{\alpha }$) in the simulated systems.

In the series of MD simulations reported here, three simulations were devoted to establishing a reference system for comparison. In the first, E2 was used to establish the structure dynamic of $\mathrm{E}2–\text{ER}\mathrm{\alpha }$. Two additional MD simulations were performed without ligand presenting at the ligand-binding site to access the characteristics of tertiary structural elements of $\text{ER}\mathrm{\alpha }$, with one having h12 in the agonist conformation, while the other was in the antagonist conformation. Both simulations may be essential, as h12 is implied to be a key player in the determination of the $\text{ER}\mathrm{\alpha }$ agonist or antagonist conformation. For each ligand, parallel simulations were performed by placing BPs into the two conformations described earlier. In all cases, helix h12 of $\text{ER}\mathrm{\alpha }$ was in an agonist conformation to match with the agonist effects described in the experimental setup. Furthermore, additional simulations were performed for E2, BPA, BPAF, and BPS ligands in a pure solvent environment (i.e., without $\text{ER}\mathrm{\alpha }$ present) so that the ligand behavior could be characterized for these ligands prior to their binding at the binding site of $\text{ER}\mathrm{\alpha }$.

Convergence of Molecular Dynamic Simulations and Ligand Dynamics

MD simulations for ligand-bound systems were performed over 300 ns, and the ligand-free systems over 650 ns. Dynamic stabilities of the systems were analyzed using the root mean square deviations (RMSD) of the atomic positions of the heavy atoms of the receptor. RMSD values were calculated by superimposing the conformations at each time point with the starting reference conformation of each system while using only the protein heavy atoms for alignment. In general, a model system is said to be unstable if RMSD values continue to rise with time. The ligand-free system that started with the helix h12 in the agonist conformation showed relatively smaller deviations over the entire 650 ns segment of the simulation, indicating that the agonist conformation of $\text{ER}\mathrm{\alpha }$ is rather stable over a relatively lengthy period, even without the ligand present in the binding site. More dynamic, but still a stable $\text{ER}\mathrm{\alpha }$ conformation was observed in the RMSD of the ligand-free system with the h12 in an antagonist conformation. Comparison of the final structures with the corresponding starting $\text{ER}\mathrm{\alpha }$ configurations clearly demonstrate that once h12 is in the agonist conformation, the overall structural features of $\text{ER}\mathrm{\alpha }$ remain quite unaltered (Figure S4). Contrarily, significant changes in the structure were observed in the antagonist conformation (Figure S5). However, larger deviations in the backbone alignment of h4, a part of the coregulator binding surface, is noticeable in both cases.

Meanwhile, all ligand-bound systems showed relatively small RMSD values ($<3\AA$), and the fluctuations of RMSD values (shown in Figure S6) are well within the acceptable range for a stable peptide system of this size with normal fluctuations in atomic positions. In general, this is a reasonable test to confirm that all ligand-bound systems are dynamic but stable during the trajectory calculations over 300 ns. $\mathrm{E}2–\text{ER}\mathrm{\alpha }$ showed relatively smaller RMSD values ($<2\AA$) over the 300-ns trajectory, and the final conformation (Figure S6) did not exhibit much divergence from the starting X-ray crystal conformation. Both $\text{BPA}–\text{ER}\mathrm{\alpha }$ and $\text{BPAF}–\text{ER}\mathrm{\alpha }$ showed similar variations ($\text{RMSD}<2\AA$) to that found in $\mathrm{E}2–\text{ER}\mathrm{\alpha }$, irrespective of ligand conformation. Meanwhile, $\text{BPS}–\text{ER}\mathrm{\alpha }$ displayed slightly more mobility in the canonical conformation of the ligand as opposed to the alternate conformation. Direct comparison of the two final configurations from each ligand-bound $\text{ER}\mathrm{\alpha }$ displayed minor changes in the tertiary structural features in ER (Figures S8–S10). Relatively larger deviations were observed among the BPS-bound systems (Figure S10), specifically h12.

RMSD values calculated from the heavy atoms of the ligands may be an important parameter in analyzing the ligand dynamics, and these values for solvated ligands are given in Figure S11. E2 practically showed no deviations from the starting configuration for over 70 ns, but both BPA and BPS structures showed rapid variations due to faster ring-flipping dynamics, while BPAF conformations were rather stable over a long duration. Only three ring-flipping events were observed for BPAF during the longer 450-ns trajectory. Because we observed how ligand RMSDs correlate with the rigidity (or flexibility) of ligands in solution, one can use the RMSDs of ligands in complex with $\text{ER}\mathrm{\alpha }$ to obtain information on ligand dynamics at the ligand-binding site during a MD trajectory calculation. For E2, very little mobility was observed, as expected from a rigid ligand in its natural $\text{ER}\mathrm{\alpha }$ ligand conformation (Figure 4A). The positions of starting and final E2 conformations overlap quite well when $\text{ER}\mathrm{\alpha }$ structures were used in the alignment (Figure 4A, left).

Meanwhile, rearrangements in conformations are evident from the ligand heavy atom RMSDs in all three BPs, but with a lesser frequency when compared to the solution structures. BPA in the canonical conformation showed almost no variation in the positions except the ring flipping (Figure 4B), whereas the alternate conformation showed much larger heavy atom positional displacements. BPAF showed almost no ring flipping for both conformations, but showed a larger positional displacement during the early part of the simulation (during the first 50 ns). This final BPAF conformation was closer to the alternate conformation than the starting canonical conformation of the ligand (Figure 4C). BPS, in the canonical conformation, showed more ring flipping than displacement, but when started in the alternate conformation, there were positional displacements of BPS that brought the final conformation more towards the canonical conformation (Figure 4D).

Energetics of Ligand Binding

We observed the potential rearrangements of BPs within the binding site of $\text{ER}\mathrm{\alpha }$ during dynamics. Therefore, using the conformations selected from the final 250-ns production runs, we performed MM/GBSA calculations to obtain the binding free energies ($\mathrm{\Delta }{\mathrm{G}}_{\text{bind}}$) of BPs to estimate how strongly different ligands bind to $\text{ER}\mathrm{\alpha }$. The binding free energy results, along with their components, are summarized in Table 2. Enthalpic contributions to the binding free energy given by $\mathrm{\Delta }{\mathrm{H}}_{\text{total},\text{GB}}$ were the dominant component in the binding interaction. In all cases, the entropic contributions given by $\mathrm{T}\mathrm{\Delta }\mathrm{S}$ were unfavorable. However, all systems still showed favorable total binding free energies. Total ligand-binding free energies of all BPs were quite comparable to the endogenous E2. However, BPA seemed to show the weakest interactions with $\text{ER}\mathrm{\alpha }$. Comparing the enthalpy contributions to the binding free energy $\mathrm{\Delta }{\mathrm{H}}_{\text{total},\text{GB}}$, two conformations of BPA showed only $2.1\text{kcal}/\text{mol}$ difference, and that was further reduced to $1.3\text{kcal}/\text{mol}$ when the binding free energies ($\mathrm{\Delta }{\mathrm{G}}_{\text{bind}}$) were compared. Similar energy behavior with slightly larger magnitudes was observed for both BPAF (about $5.9\text{kcal}/\text{mol}$) and BPS (about $3.7\text{kcal}/\text{mol}$), showing greater enthalpic contributions ($\mathrm{\Delta }{\mathrm{H}}_{\text{total},\text{GB}}$) to the binding free energy when compared with the values for configuration 1. However, BPAF and BPS configurations were favored only by 4.6 and $0.9\text{kcal}/\text{mol}$ when the binding free energies were compared. At this point, looking at the final conformations from dynamics and taking the free energies into consideration, we confirmed the fact that BPA and BPS favor conformation 1 that resemble the X-ray crystal conformation of BPA in $\text{ER}\mathrm{\alpha }$, whereas BPAF prefers conformation 2 in the complex that differs from the former in its $\text{ER}\mathrm{\alpha }$ residue interactions.

Structures and Dynamics of Ligand Binding and a Close Inspection of the Coregulator Binding Surface

During the examination of ligand RMSDs, we concluded that there were rearrangements in ligand conformations. If the ligands undergo conformation changes, can they change the dynamics of interacting $\text{ER}\mathrm{\alpha }$ residues by searching for new interacting partners when accommodating changes in ligand conformations? Using the conformations of ligands and their ER interacting partners on the binding site in the X-ray crystal structures (Figure S2) as the reference, we systematically analyzed the ligand–receptor interactions during dynamics. Representative ligand conformations near the end of each trajectory are shown in Figure 5. In the $\mathrm{E}2–\text{ER}\mathrm{\alpha }$ conformation, the interactions between R394 and E353 of $\text{ER}\mathrm{\alpha }$ was enhanced during the simulation. The strengthening of the salt bridge between R394 and E353 may have weakened their interactions with the phenolic oxygen of E2 displayed in the X-ray crystal structure, but the phenolic oxygen (of E2) remained in the close proximity to these two residues during the entire time of the simulation. The other oxygen of E2 seemed to make more stable contacts with H524, G521, and L525 (Figure 5A), and these contacts remained intact during the entire simulation. It is noticeable (but not surprising, due to its bulkiness) that E2 displayed the largest number of interacting $\text{ER}\mathrm{\alpha }$ residues (18), and this number remained relatively un-changed (with residues T347 and M528 alternatively appearing in the coordinating shell). Between 13 and 17 $\text{ER}\mathrm{\alpha }$ residues were involved in direct interactions with BPA, BPAF, or BPS. Due to the flexibility and variable hydrophobicity of BPs, variations in $\text{ER}\mathrm{\alpha }$ residue participation in the stabilization of these compounds at the ligand-binding site were observed via close examination of the ligand interactions (Figure 5B–D). This variability of the bisphenol ligand interactions with $\text{ER}\mathrm{\alpha }$ may have led to the introduction of significant mobility into the assimilation of the coregulator binding surface formed by the helices h3, h4, and h12. In addition to the multiple BPAF configurations observed in the crystal structures (Delfosse et al. 2012), we found that the other two BPs were also capable of accommodating alternate binding conformations when they occupied the ligand-binding site of $\text{ER}\mathrm{\alpha }$.

Next, we focused on the variations in electrostatic structures of the coregulator binding surface defined by the helices h4, h12, and a part of h3. From the available X-ray crystal structures, the location of the cofactor binding surface was identified. Since we recognized the coregulator binding surface on $\text{ER}\mathrm{\alpha }$, changes in structures and other associated properties of this surface during dynamics were used to characterize the coregulator binding. Information gathered from electrostatic potentials (ESP) mapped onto the surface of the molecule was a reasonable tool for predicting what type of interactions the molecule may have with an incoming molecule. Differentiation of coregulator binding due to variable ligand occupancy at the ligand-binding site can be carried out by comparing the ESP of the above surface. ESP surfaces calculated using the Amber charge distribution for the atoms in the amino acid residues of $\text{ER}\mathrm{\alpha }$ systems with different ligands at the ligand-binding site are given in Figure 5E (on the right side of each ligand-binding site conformation). In the figures showing electrostatic potential energy surfaces (Figures 5A–D, right), coregulator interaction areas bracketed by helices h3, h4, and h12 are circled in white. For $\mathrm{E}2–\text{ER}\mathrm{\alpha }$ (shown in Figure 5A), the basic residues K362 and R363 identified in blue depicted in the top part of the circle interacted with some coregulator residues. The acidic residues E380, D538, E542, and D545 (represented in red in the bottom of the circle) have also been shown to interact with coregulator residues. The middle region stayed charge neutral to facilitate the $\text{ER}\mathrm{\alpha }$ interaction with the coregulator binding motif. All three bisphenol compounds preserved general features of their electrostatic potential of this area suitable for cofactor binding (Figures 5B–D). However, the general electrostatic potential features of the coregulator interacting surfaces of the ER with various BPs were comparable to those of $\mathrm{E}2–\text{ER}\mathrm{\alpha }$; some changes in the fine details were observed in the ESP surface of the relevant area, especially when the same ligand occupied the ligand-binding site in two different conformations.

The three main helices, h3, h4, and h12, are presented in Figure 5E that define the cofactor binding surface. The embedded leucine residues in the ${\mathrm{L}}_1{\mathrm{X}}_2{\mathrm{X}}_3{\mathrm{L}}_4{\mathrm{L}}_5$ motif of the SRC-1 peptide are shown using the space-filling model. The space for leucine binding is preserved in all systems in their representative structures, as in Figure 5E, and variations in the position of the backbone h4 and h12 can also be observed. Interactions of the ${\mathrm{X}}_2$ residue with helix h4 were altered due to the mobility observed in that helix region. Also, the residues in the N-terminal side from ${\mathrm{L}}_1$ of the SRC-1 peptide may potentially be in direct contact with h12 (specifically with the negatively charged residues marked in Figure 4), and even small variations in the position of h12 may lead to significant alteration in the binding strength of the coregulator peptide and $\text{ER}\mathrm{\alpha }$. Therefore, taken together, the information from the electrostatic potential surfaces and the variations in the side chain locations of potentially cofactor interacting residues in helices h3, h4, and h12 supports the observed variations in binding affinity of various coregulators, even when the same ligand is present in the ligand-binding site.

Discussion

The ER ERE–Mediated Transcriptional Activity and Alteration of ER Target Genes by BPA, BPAF, and BPS

A recent review concluded BPA is a reproductive toxicant based on a series of studies reporting that BPA impacts female reproduction and potentially affects the male reproductive system in humans and animals (Peretz et al. 2014). As BPA substitutes are replacing BPA in manufacturing, these analogs, like BPAF and BPS, are now more prevalent in consumer products (Liao et al. 2012a; Liao et al. 2012b; Liao and Kannan 2014). Our studies concluded that BPA and BPAF can function as EDCs by acting as cell type–specific agonists ($\ge 10\text{nM}$) or antagonists ($\le 10\text{nM}$) for mouse $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$ (Li et al. 2012). In the comparison of twelve EDCs, including BPA and BPAF, we found that EDCs with similar chemical structures tended to have comparable mouse $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$ ERE-mediated activities. However, similar chemical structures did not correlate with ligand-binding affinities of the EDCs to the ERs (Li et al. 2013). In the present study, we found that BPA and BPAF activate ERE-mediated transcription for human $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$, and the results were comparable due to the high homology between the DNA sequences of ERs in mouse and human. Interestingly, the estrogenic activity of BPS is more $\text{ER}\mathrm{\alpha }\text{-}\text{specific}$ in the cell systems we used. Human WISP-2/CCN5 is an important regulator involved in the maintenance of a differentiated phenotype in breast tumor epithelial cells (Fritah et al. 2008). Promoter analysis showed that endogenous CBP and p21 were recruited in an E2-dependent manner to the chromosomal human WISP-2/CNN5 promoter, which contains an ERE site in MCF-7 cells (Fritah et al. 2006). Using an ERE site in the human WISP-2/CNN5 gene promoter as a tool, we found that BPA, BPAF, and BPS induced $\text{ER}\mathrm{\alpha }$ direct binding to the ERE site of this gene promoter. Using Ishikawa cells stably expressing mouse $\text{ER}\mathrm{\alpha }$, we found that ER target genes, including GREB1, PGR, and WISP-2/CNN5, were increased by BPA and BPAF (Li et al. 2012, 2013). Here, we confirmed the effect of BPA, BPAF, and BPS on three ER target genes in MCF-7 cells, suggesting these target genes may be commonly upregulated by BPs.

Differential Coregulator Recruitment in the $Ligand–\mathit{\text{E}}\mathit{\text{R}}\mathrm{\alpha }$ Complex

Previous studies showed that ER mediates transcription using a common cohort of coregulators through many target genes and recruitment was based on ER confirmation (DeVilbiss et al. 2013). Coregulators are intrinsically disordered proteins that become structured upon interaction with other proteins. This flexibility allows a single coregulator to harbor many TF binding motifs that cause intrinsically disordered regions of the protein (Millard et al. 2013). In this study, we used a peptide chip array that contains the most well-characterized 154 NR interaction surfaces from 66 coregulators. These included coactivators with highly conserved helical LXXLL (NR box) and corepressors with highly conserved helical LXXXIXXXL (CoRNR box) motifs, and were used to investigate the modulation of coregulators in the $\text{BP}–\text{ER}\mathrm{\alpha }$ complex. Our findings indicated that coactivators such as SRCs and PGCs were differentially recruited in the presence of BPA, BPAF, and BPS. Interestingly, there are more corepressors, such as small heterodimer partner (SHP) and receptor-interacting protein 140 (RIP140), recruited by the $\text{BPS}–\text{ER}\mathrm{\alpha }$ complex; this may explain why BPS has weaker estrogenic activation than BPA and BPAF in our in vitro system. Of note, unlike other coregulators that are known either for the coactivating or co-repressing roles, RIP140 is capable of acting both as a coactivator and a co-repressor in the different biological systems (Nautiyal 2017). From our peptide studies, we noticed the differential bindings of nuclear receptor-interacting protein 1 (NRIP) peptides, which were deviated from RIP140 in the presence of BPA, BPAF, and BPS. Furthermore, a more recent study indicated that BPS had more than tenfold activity/toxicity in vivo compare to in vitro (Conley et al. 2016). Our new findings may also enhance our understanding of differential molecular mechanisms between the in vivo and in vitro studies and the biological activities in a ligand- and tissue-specific manner.

Dynamic Simulations of Bisphenol Chemicals Bound to $\mathit{\text{ER}}\mathrm{\alpha }$

Binding interactions of several bisphenol compounds and the endogenous ligand E2 with $\text{ER}\mathrm{\alpha }$ were characterized using the solution structure evaluations by MD simulations. E2 was recognized as a ligand with relatively strong binding interactions and was used as a reference system. Among the bisphenol compounds studied, BPA showed weaker interactions with $\text{ER}\mathrm{\alpha }$ at the binding site, while BPAF and BPS showed comparable interaction strengths to E2 in their binding conformations. However, energetics and structures of bisphenol compounds display variable binding modes with differing strengths at the ligand-binding site of ER.

Ligand-free simulations are quite instructive in estimating stabilities brought into $\text{ER}\mathrm{\alpha }$ due to ligand binding. In addition, they can provide valuable information on intrinsic stabilities of tertiary structural motifs of $\text{ER}\mathrm{\alpha }$. The major structural change that is said to be associated with the ligand-binding is the relocation of the h12 helix to create the coregulator binding surface. Our simulation, starting with the agonist form of $\text{ER}\mathrm{\alpha }$, showed a rather stable structure around the coregulator binding site bracketed by helices h3, h4, and h12 over the 650-ns simulations, even in the absence of any ligand in the $\text{ER}\mathrm{\alpha }$ ligand-binding site. Once positioned, h12 in the agonist conformation seemed to be in a rather stable conformation compared to the other conformations for this helix in $\text{ER}\mathrm{\alpha }$. However, some noticeable deviations of h4 can be seen in the backbone alignment during dynamics.

BPA, BPAF, and BPS show varying conformations in the ligand-binding site of ER. Dynamical trajectory calculations of bisphenol ligand-bound $\text{ER}\mathrm{\alpha }$ systems showed no major modifications in the coregulator-binding surface due to ligand mobility, even with a variable conformation space for ligands at the ligand-binding site. However, these minor alterations may be sufficient to introduce variability to the cofactor binding strengths specifically due to the repositioning of some key side chains in helices h4 and h12. Further studies are underway to characterize the antagonist features of BPs interacting with $\text{ER}\mathrm{\alpha }$.

Conclusions

In summary, our data indicate that BPA, BPAF, and BPS have varying estrogenic activity for $\text{ER}\mathrm{\alpha }$ and $\text{ER}\mathrm{\beta }$, except BPS, which did not show activation of $\text{ER}\mathrm{\beta }$. We also conclude that coregulators are differentially recruited in the ligand–$\text{ER}\mathrm{\alpha }$ complex when bound with BPA, BPAF, and BPS, which may help explain their differential biological actions. This is the first report to investigate coregulator recruitment of the $\text{bisphenol}–\text{ER}\mathrm{\alpha }$ complex and the subsequent discovery that BPs stimulate formation of more corepressor complexes than E2. Molecular modeling suggests BPs can bind to the ER–ligand-binding domain with differing energetics and conformations. Furthermore, we characterized the surface of coregulator binding on $\text{ER}\mathrm{\alpha }$ in complexes with bisphenol compounds using molecular modeling results. Taken together, these findings provide an important basis for the understanding of the molecular mechanisms of BP activity in ER-mediated activation and the differential biological activities of these EDCs that are ubiquitous in our everyday products.