Differential in Vitro Biological Action, Coregulator Interactions, and Molecular Dynamic Analysis of Bisphenol A (BPA), BPAF, and BPS Ligand–ERα Complexes

Background: Bisphenol A (BPA) is an endocrine-disrupting chemical (EDC) that might be harmful to human health. Recently, there has been widespread usage of bisphenol chemicals (BPs), such as bisphenol AF (BPAF) and bisphenol S (BPS), as replacements for BPA. However, the potential biological actions, toxicity, and the molecular mechanism of these compounds are still poorly understood. Objectives: Our objective was to examine the estrogenic effects of BPA, BPAF, and BPS and the molecular mechanisms of action in the estrogen receptor alpha (ERα) complex. Methods: In vitro cell models were used to compare the estrogenic effects of BPA, BPAF, and BPS to estrogen. Microarray Assay for Real-Time Coregulator-Nuclear receptor Interaction (MARCoNI) analysis was used to identify coregulators of BPA, BPAF, and BPS, and molecular dynamic (MD) simulations were used to determine the compounds binding in the ERα complex. Results: We demonstrated that BPA and BPAF have agonistic activity for both ERα and ERβ, but BPS has ERα-selective specificity. We concluded that coregulators were differentially recruited in the presence of BPA, BPAF, or BPS. Interestingly, BPS recruited more corepressors when compared to BPA and BPAF. From a series of MD analysis, we concluded that BPA, BPAF, and BPS can bind to the ER–ligand-binding domain with differing energetics and conformations. In addition, the binding surface of coregulator interactions on ERα was characterized for the BPA, BPAF, and BPS complexes. Conclusion: These findings further our understanding of the molecular mechanisms of EDCs, such as BPs, in ER-mediated transcriptional activation, biological activity, and their effects on physiological functions in human health. https://doi.org/10.1289/EHP2505


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), ERa and ERb (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 ligandbound 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 proliferatoractivated receptor gamma coactivator 1s (PGC1s), including PGC1a and PGC1b, 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 ERa responses. Using Chromatin Immunoprecipitation (ChIP) analysis, we determined direct ERa binding occurs in the presence of BPA, BPAF, and BPS to the ERE of an ER target gene promoter. We profiled the recruitment of coregulatorderived peptides in the ligand (BPA, BPAF, or BPS)-ERa 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 ERa ligand complex through computational analysis.

Chemicals
The 17 b-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-mM stock solutions were made in DMSO and kept at −20 C.

Cell Lines and Tissue Culture
The human endometrial adenocarcinoma stable cell lines Ishikawa/ vector (Ishikawa/vec) and Ishikawa=WT ERa (Ishikawa=ERa) have been described previously (Burns et al. 2011(Burns et al. , 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 . The human bone osteosarcoma epithelial recombinant U2OS cells stably expressing enhanced green fluorescent protein (EGFP)-tagged human ERa fused to the C-terminus of EGFP were purchased from Thermo Fisher Scientific. Ishikawa/vec and Ishikawa=ERa 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 mg=mL Geneticin ® (G418; Invitrogen). HepG2 cells were maintained in phenol red-free Minimum Essential Medium (MEM; Invitrogen) supplemented with 10% FBS and 4 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 mM L-glutamine. Recombinant U2OS=ERa cells were maintained in phenol red-free DMEM supplemented with 10% FBS and 4 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.
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 ± standard error of the mean (SEM).

RNA Extraction and Real-Time Polymerase Chain Reaction Analysis
Cells were seeded in 60 × 15 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 nM E2, or 100 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 b-actin expression. Experiments were repeated three times, and results are presented as mean ± SEM. Changes of endogenous gene expression by E2 were calculated Environmental Health Perspectives 017012-2 as fold change relative to the vehicle control group in each cell line.

Cell Lysis Preparation and Microarray Assay for Real-Time Coregulator-Nuclear Receptor Interaction Assay
Recombinant U2OS cells stably expressing human ERa fused to the C-terminus of EGFP were maintained in phenol red-free DMEM in the presence of 2 mM L-Glutamine, 1% penicillinstreptomycin, 0:5 mg=mL G418, and 10% FBS at 37°C, 5% 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 × 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 ERa-transfected U2OS cells. For MARCoNI, 25 lL assay mixtures that contain cell lysates, 25 nM of Alexa Fluor 488-conjugated GFP antibody (Invitrogen), and 0:2 mM ligand (E2, BPA, BPAF or BPS, prediluted in DMSO, final concentration, 2%) were prepared on ice. The ligand-modulated coregulator interactions with the EGFPtagged ERa 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

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)(468)(469)(470)(471)(472)(473)(474)(475)(476)(477)492,531,and 545 in 3UUD;306,330,419,(467)(468)(469)(470)(471)(472)(531)(532)(533)and 548 in 3UU7;and 306,330,334,335,337,339,342,363,413,416,417,419,423,(466)(467)(468)(469)(470)(471)(472)(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 kcal=mol=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) (default parameter 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.

Estrogenic Activation of ERa and ERb by BPA, BPAF, and BPS, Alterations in Expression of ER-Target Genes, and Direct Binding of ERa to an ERE Located in the Promoter of ER-Regulated Genes
To evaluate the ERE-mediated transcriptional activity of ERa 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 ERa cells (Burns et al. 2011(Burns et al. , 2014. As shown in Figure 1B, 1 nM E2 induced ERE-mediated activation by fourfold. BPA-induced ERE-mediated activation by threefold at 100 nM. BPAF induced ERE-mediated activation by threefold at 10 nM, and matched the efficacy of E2 at 1,000 nM. BPS only induced ERE-mediated activation by threefold at the highest concentration (1,000 nM) in Ishikawa=ERa cells. ERE-mediated activation was weak or absent in Ishikawa/vector cells, suggesting that the activation by BPA, BPAF, and BPS are predominantly ERa-dependent. To investigate the estrogenic activity in the different cell types through an endogenous receptor, human breast cancer MCF-7 (ERa positive) and ovarian cancer BG-1 FR (ERa and ERb positive) cells (Burns et al. 2014;Li et al. 2014) were used in reporter assays. In MCF-7 cells, 1 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 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 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 ERa-and ERb-mediated transcriptional activation, we performed reporter assays in human hepatocellular cancer HepG2 (ER negative) cells after transfecting with an ERa or ERb expression plasmid. For ERa ERE-mediated activation, cells were highly responsive to E2, with up to 35-fold increases at 1 nM ( Figure 1D, left). BPA induced this activation by 25-fold at 100 nM. BPAF induced activation by 17-fold at 10 nM and 30fold at 100 nM. In contrast, BPS only weakly activated the ERa ERE at 100 nM (eightfold). With the highest concentration of 1,000 nM, all three chemicals displayed strong activation of ERa by greater than 30-fold. For ERb, E2 induced ERE-mediated activation with a lower fivefold increase at 1 nM ( Figure 1D, right). BPA at 100 nM and BPAF at 10 nM had weaker activation when compared to E2 at the same concentrations. BPA at 1,000 nM and BPAF at 100-1,000 nM had strong activity for ERb EREmediated activation. However, BPS did not have ERb EREmediated activation in HepG2 cells. These results indicate that BPA and BPAF can activate ERE-mediated transcription via ERa and ERb, but BPS exhibited more ligand specificity for ERa.
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.

ERa-Coregulator Interaction Profiling and Comparison of BPA, BPAF, and BPS Coregulator-Derived Binding Peptides on ERa
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 ERa 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 ERa (Figure 2A). As a parameter for modulation of the interaction of ERa 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 ligand-ERa 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 ERa function, the recruitment of coregulator motifs, as induced by BPA, BPAF, and BPS, were compared to that observed with E2 stimulation of the ERa. 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 ERa 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 BPS-ERa complexes exhibited more interactions than with BPA or BPAF, but these binding activities were generally lower than with E2 in the ERa complexes, and in certain examples, selective to the BPS . Co-transfection with ERa or ERb in HepG2 cells. Cells were transfected with 3xERE-luc and pRL-TK plasmids in MCF-7, BG-1 FR, or co-transfection with the ERa or ERb expression plasmid in HepG2 cells overnight. After changing to fresh starved medium, cells were treated with vehicle control, 1-1,000 nM estradiol (E2), BPA, BPAF, and BPS for 18 h. The ER ERE-mediated activity was detected by the luciferase reporter assay as described in the materials and methods. Data shown are representative of triplicates as fold increases that were calculated relative to the vehicle ðcontrolÞ ± standard error of the mean ðSEMÞ, ***p < 0:001, **p < 0:01, *p < 0:05 compared to vehicle control. (E) Changes of the ER target genes in MCF-7 cells. Cells were treated with vehicle control, 10 nM E2, or 100 nM BPA, BPAF, and BPS, with or without 10 lM ICI 182,780 (ICI) for 18 h, respectively. Total RNA was extracted and used as a template for cDNA synthesis. Gene expression was measured by quantitative PCR (qPCR). Experiments were repeated three times, and results are presented as mean ± SEM, **p < 0:01, ***p < 0:001, ****p < 0:0001. (F) BPA-, BPAF-, and BPS-induced ERa binding at an ERE site of the human WISP-2 gene promoter. Cells were treated with vehicle control, 10 nM E2, and 100 or 1,000 nM BPA, BPAF, or BPS for 40 min. ChIP assay was carried out with anti-ERa antibody using a ChIP Assay Kit as described in the materials and methods. Purified chromatin DNA was used for qPCR analysis by the special ChIP primers for the amplicons −699= − 422 (AP-1 + ERE) or −593= − 422 (ERE). Relative recruitments of ERa were normalized to signals obtained from input DNA. Experiments were repeated three times, and results are presented as mean ± SEM, **p < 0:01, ***p < 0:001, ****p < 0:0001. 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-ERa 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 ligand-ERa complexes are shown in Figure 2D. E2 was found to significantly modulate ERa 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 ligand-ERa 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 BPA-ERa bound complex (Tables S6-S9). The potential coregulators (from negative interaction peptides) in the ligand-ERa 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

Differential Recruitments of Coactivators SRCs and PGCs in Ligand-ERa Complexes
To evaluate the peptide screening results, we selected two wellcharacterized coactivator families, SRC and PGC1, for further experimentation. To determine the role of ERa with coregulators in ERE-mediated transcriptional regulation using reporter assays, we co-transfected SRCs or PGC1s with ERa in HepG2 cells. As expected, SRC-1, SRC-2, SRC-3, PGC1a, and PGC1b induced ERa ERE-mediated activity when treated with 10 nM E2 ( Figure  3A). For SRCs, in the presence of 100 or 1,000 nM BPA, SRC-1 significantly induced ERa 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 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 nM BPS, all three SRCs significantly activated ERa EREmediated transcription. However, there was no coactivation with either SRCs at 100 nM. PGC1a significantly induced ERa EREmediated activity at high concentrations of 1,000 nM BPA, and PGC1b was active by BPA at the lower concentration with 100 nM ( Figure 3C). PGC1a also significantly activated the ERa ERE-mediated transcription in the presence of 100 or 1,000 nM BPAF. PGC1b, however, further induced activity only with 1,000 nM BPAF. BPS showed a clear differential preference of activation with PGC1a in the presence of 1,000 nM BPS, but PGC1b 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 ligand-ERa complexes.

Simulating BPA-, BPAF-, and BPS-Bound ERa by Molecular Dynamics
To investigate if any potential modifications to structure appear on the coregulator binding surface of ERa-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 ERa at the ligand-binding site (Delfosse et al. 2012). In that work, BPA was bound to ERa 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 E2-ERa ( Figure S2B). In addition to the canonical conformation ( Figure S2C), for BPAF-ERa, 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 BPAF-ERa, even though BPAF conformations were different, h12 was observed in a similar conformation to BPA-or E2-ERa (i.e., ERa 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 BPC-ERa. Interestingly, in the case of BPC-ERa, h12 was located elsewhere with a conformation labeled as an antagonist form of ERa. It is worth noting that the ERa 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 ERa) 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 E2-ERa. Two additional MD simulations were performed without ligand presenting at the ligand-binding site to access the characteristics of tertiary structural elements of ERa, 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 ERa 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 ERa 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 ERa present) so that the ligand behavior could be characterized for these ligands prior to their binding at the binding site of ERa.

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  BPAF, or BPS. Cells were transfected with 3xERE-luc, pRL-TK, pcDNA=ERa, and pcDNA/SRC1, pcDNA/SRC2, pcDNA/SRC3, pcDNA=PGC1a, or pcDNA=PGC1b plasmids overnight. After changing to fresh starved medium, cells were treated with the vehicle control, 10 nM E2, 100 or 1,000 nM BPA, BPAF, and BPS for 18 h. The ERa ERE-mediated activity was detected by a luciferase reporter assay. Data are representative of triplicates as fold increase calculated relative to the vehicle ðcontrolÞ ± standard error of the mean ðSEMÞ, ***p < 0:05 compared to vehicle (control). unstable if RMSD values continue to rise with time. The ligandfree 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 ERa is rather stable over a relatively lengthy period, even without the ligand present in the binding site. More dynamic, but still a stable ERa 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 ERa configurations clearly demonstrate that once h12 is in the agonist conformation, the overall structural features of ERa 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 A), 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. E2-ERa showed relatively smaller RMSD values (<2 A) over the 300-ns trajectory, and the final conformation ( Figure S6) did not exhibit much divergence from the starting X-ray crystal conformation. Both BPA-ERa and BPAF-ERa showed similar variations (RMSD <2 A) to that found in E2-ERa, irrespective of ligand conformation. Meanwhile, BPS-ERa 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 ligandbound ERa displayed minor changes in the tertiary structural 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 ERa 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 ERa ligand conformation ( Figure 4A). The positions of starting and final E2 conformations overlap quite well when ERa 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 ERa 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 (DG bind ) of BPs to estimate how strongly different ligands bind to ERa. The binding free energy results, along with their components, are summarized in Table 2. Enthalpic contributions to the binding free energy given by DH total,GB were the dominant component in the binding interaction. In all cases, the entropic contributions given by TDS 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 ERa. Comparing the enthalpy contributions to the binding free energy DH total,GB , two conformations of BPA showed only 2:1 kcal=mol difference, and that was further reduced to 1:3 kcal=mol when the binding free energies (DG bind ) were compared. Similar energy behavior with slightly larger magnitudes was observed for both BPAF (about 5:9 kcal=mol) and BPS (about 3:7 kcal=mol), showing greater enthalpic contributions (DH total,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 kcal=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 ERa, whereas BPAF prefers conformation 2 in the complex that differs from the former in its ERa 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 ERa 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 E2-ERa conformation, the interactions between R394 and E353 of ERa was enhanced during the simulation. The Table 1. Coactivator steroid receptor coactivators (SRCs) and peroxisome proliferator-activated receptor gamma coactivator 1 (PGC1) in the ligand-ERa complex by Microarray Assay for Real-Time Coregulator-Nuclear receptor Interaction (MARCoNI) assay. Environmental Health Perspectives 017012-11 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 ERa residues (18), and this number remained relatively unchanged (with residues T347 and M528 alternatively appearing in the coordinating shell). Between 13 and 17 ERa residues were involved in direct interactions with BPA, BPAF, or BPS. Due to the flexibility and variable hydrophobicity of BPs, variations in ERa 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 ERa 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 ERa.
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 ERa, 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 ERa 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 E2-ERa (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 ERa 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 E2-ERa; some changes The black curve represents the system starting with the ligand in its canonical conformation (conformation 1), and the red curve represents the system started with the ligand in its alternate conformation (conformation 2). In each case, the estrogen receptor (ERa) is in the agonist form. The final ligand structures are also shown (in dark ball and stick model) with their starting conformations (in lighter stick model). Table 2. Binding free energies (DG Bind ) and their enthalphic (DH total,GB ) and conformational entropic (TDS) components (in kcal/mol) of ligands in ERa predicted by MM/GBSA method.
in the fine details were observed in the ESP surface of the relevant area, especially when the same ligand occupied the ligandbinding 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 L 1 X 2 X 3 L 4 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 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 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 ERa. 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 (≥10 nM) or antagonists (≤10 nM) for mouse ERa and ERb (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 ERa and ERb 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 ERa and ERb, 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 ERa-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 ERa direct binding to the ERE site of this gene promoter. Using Ishikawa cells stably expressing mouse ERa, we found that ER target genes, including GREB1, PGR, and WISP-2/CNN5, were increased by BPA and BPAF (Li et al. 2012(Li et al. , 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-ERa 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 BP-ERa 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 BPS-ERa 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 ERa
Binding interactions of several bisphenol compounds and the endogenous ligand E2 with ERa 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 ERa 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 ERa due to ligand binding. In addition, they can provide valuable information on intrinsic stabilities of tertiary structural motifs of ERa. 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 ERa, 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 ERa 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 ERa. 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 ERa 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 ERa.

Conclusions
In summary, our data indicate that BPA, BPAF, and BPS have varying estrogenic activity for ERa and ERb, except BPS, which did not show activation of ERb. We also conclude that coregulators are differentially recruited in the ligand-ERa 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 bisphenol-ERa 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 ERa 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 ERmediated activation and the differential biological activities of these EDCs that are ubiquitous in our everyday products.