The Effect of Structural Diversity on Ligand Specificity and Resulting

Subscriber access provided by BUFFALO STATE

Article

The Effect of Structural Diversity on Ligand Specificity and Resulting Signaling Differences of Estrogen Receptor # Qiao Xue, Xian Liu, Xiuchang Liu, Wenxiao Pan, Jianjie Fu, and Aiqian Zhang Chem. Res. Toxicol., Just Accepted Manuscript • DOI: 10.1021/acs.chemrestox.8b00338 • Publication Date (Web): 29 Mar 2019 Downloaded from http://pubs.acs.org on April 1, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemical Research in Toxicology

The Effect of Structural Diversity on Ligand Specificity and Resulting Signaling Differences of Estrogen Receptor α Qiao Xue1, Xian Liu1, Xiu-Chang Liu1,2, Wen-Xiao Pan1, Jian-Jie Fu1, Ai-Qian Zhang1,2,3* Address: 1State Key Laboratory of Environmental Chemistry and Ecotoxicology, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing 100085, P. R. China. 2College

of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, P. R. China

3

Institute of Environment and Health, Jianghan University, Wuhan 430056, P.R. China

ACS Paragon Plus Environment

Chemical Research in Toxicology 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Contents


Page 2 of 51

Page 3 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ABSTRACT

Numerous chemicals have been reported to exert estrogen-like endocrine disrupting effects via a receptor binding mechanism that directly interacts with the ligand binding domain (LBD) of estrogen receptor α (ERα). However, not only their binding affinities to ERα but also their interference in specific cell and tissue functions are clearly different. In this regard, significant regulation differences among three representative estrogenic chemicals (diethylstilbestrol (DES), bisphenol A (BPA) and diarylpropionitrile (DPN)), well-known ERα agonists with very similar structures, have been observed. Molecular dynamics (MD) simulation is used to explore the underlying mechanism of different regulation effect induced by the similar estrogen-like chemicals. DES induced 12 Å motion of the H9-H10 loop markedly expands the negative electrostatic potential surface of the AF-2 domain, which is consistent with the over-regulation effect of the agonist. In comparison, the 3 Å motion induced by BPA and DPN corresponds to the low-regulation effect of the chemicals. Cross-correlation analysis indicates that the different ERα motions and resulting surface feature of AF-2 domain is brought by the distinguished binding modes of the agonists. Moreover, only hydrophobic DES with estrogen-like size and flexibility has high binding affinity of -23.47 kcal/mol binding free energy. Both the hydrophilic group in DPN and the small molecular size of BPA dramatically decrease the agonist binding ability, and their binding free energies are only -12.43 kcal/mol and -11.82 kcal/mol, respectively. Our study demonstrates that similar chemicals interact



Page 4 of 51

differently with ERα and induce different allosteric effects, which explains the observed regulation diversity.

KEYWORDS Endocrine disrupting chemicals, estrogen receptor α, molecular dynamics simulations, regulation effect

INTRODUCTION Endocrine disrupting chemicals (EDCs) interfere with the hormone system in mammals and may adversely affect the development, reproduction, transport, and metabolism of hormones.1, 2 An important category of EDCs is environmental estrogens, which display estrogenic

activity

and

interfere

with

normal

estrogen

signaling

pathway.3

Diethylstilbestrol (DES) and bisphenol A (BPA) are well known environmental estrogens (EEs). DES has been used to prevent miscarriage, premature labor, and related complications of pregnancy. However, it has been proven to increase the risk of breast cancer and breast cancer mortality. Moreover, the daughters of the women who are exposed to DES during pregnancy are more likely to suffer vaginal/cervical clear cell adenocarcinoma.4 BPA is a common industrial material that is used to produce plastics and epoxy resins, and its adverse health effects in humans are now being recognized.5 BPA exhibits estrogenic properties and can affect the reproductive process by mimicking


Page 5 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


endogenous estrogens.6,

7

Previous research has identified that DES and BPA produce

health hazard effects by directly interact with estrogen receptor α (ERα).8, 9 Generally, ERs contains two subtypes, ERα and ERβ, which both have unique and complementary roles in estrogen signaling pathways.10,

11

These two ERs have similar

structures with three functional domains: the N-terminal transactivation domain, the central DNA binding domain, and the C-terminal ligand binding domain (LBD).12 The ER-LBD contains 12 α-helices (H1-H12) and a short two-stranded antiparallel β-sheet (s1 and s2) arranged in a three-layered sandwich structure (Figure S1). Such arrangement yields a hydrophobic pocket, comprised of H3, H6-H8, H11, H12, s1 and s2. Crystal structures show that most of the EEs initially bind with the pocket.13,

14

Notably, this

molecular-driven step initiates the entire signaling pathway, suggesting that the pocket binding mechanism is a crucial factor in the relevant endocrine disrupting process. When the ligands are captured in the binding pocket of the ER-LBD, the apo conformation of the unbound LBD adopts either an agonism conformation or antagonism conformation. Experiments have shown that conformational diversity triggered by either behavior can affect the signaling pathway and induce different gene expression.15, 16 Thus, identifying the structural influence of ligand binding will be helpful for studying the estrogenic effects of ER ligands. Agonists can induce allosteric regulation of H12, thereby forming of activation function 2 (AF-2), which would complete the recruitment of coactivators and activate the transcription.12 The remaining question is why the regulation effect of the agonists is different.



Both BPA and DES are ER agonists, but their ER-mediated effects are quite different. BPA has significant effect on expressions of mRNAs and the sex-determining gene in the mouse urogenital sinus, while DES has little effect on them.17 Another common ER agonist, diarylpropionitrile (DPN), which can also directly bind to ER-LBD to modulate the transcriptional activity of ER, also has its own regulation characteristics.18, 19 The chemical structures of the three agonists are very similar (Figure 1), but their ER binding affinities are significantly different. DES has the highest binding affinity, followed by DPN and BPA.20, 21 Unfortunately, the difference in binding affinity cannot be readily explained by the obtained static structures of the LBD.8, 12, 22 In fact, the cocrystalized structure indicates that the LBD conformations of the ER complexes with BPA, DES and 17β-estradiol (E2) share the common feature, and the position of H12, which is thought to be critical for regulation, is highly conservative. But the agonism effects of the three chemicals on ER vary and are different from those of endogenous hormone E2 as well. Exposure to DES leads to an increase in aggressive behavior and has little impact on male sexual activity.23 In contrast, the aggressiveness of rats is not influenced by BPA and it slightly impairs the male sexual activity of rats.24 The effect of DPN is more different than DES and BPA. BPA significantly affects the proliferation of MCF-7 cells, while DPN has little effect on it.25 Moreover, DES can induce vaginal epithelial cell proliferation which is not observed for DPN.26 Thus, what is baffling is significant differences among the transcriptional effects and binding affinities of the three agonists of similar structures. To date, the obtained information on both


Page 6 of 51

Page 7 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


structural biology and experimental data cannot support the elucidation of ER regulation diversity. The relationship among the ligand structure, binding affinity, protein structure and the final regulation effect is unclear. There must be an intricate regulatory mechanism, not just the positioning of H12, in controlling the agonist induced regulation of ERα. Detailed analysis of the binding mode of an individual ligand and its influence on LBD structure is vital to understanding how small ligands with similar structures govern different regulation effect. In recent years, molecular dynamics (MD) simulations combined with molecular mechanics – generalized born surface area (MM-GBSA) calculations have been successfully used to investigate recognition mechanism between nuclear receptors (NRs) and ligands.27-29 This method can provide not only a wealth of dynamic structural information about the protein-ligand interactions but also detailed energy information.30, 31

The dynamic atomic level information offers us the opportunity to explore the

allosteric mechanism of ERα. The previous study successfully revealed the structural basis of receptor selectivity for the contaminants such as BPA analogues, bisphenol S analogues and DDT analogues, while little effort has been made to investigate the molecular mechanism of regulation diversity induced by the agonist binding to the same receptor.32-34 Our work has been carried out to obtain the effect of structural diversity on ligand specificity and resulting signaling differences of ERα to fill in the gap. Thus, we construct the complexes of human ERα-LBD with E2, DES, BPA and DPN, and present the results of 160 ns all-atom MD simulations. The focus is to analyze the regulation



diversity induced by the different binding mechanism of DES, BPA and DPN. We will compare the differentiated motion of ERα-LBD for three EEs with that of the endogenous estrogen E2 to indicate the factors that influence the binding mode and regulation effect. METHODS Preparation of Molecular Systems The initial structures were taken from the Protein Data Bank (PDB). The crystal structures of the ERα-E2 (PDB ID: 1ERE12) and ERα-DES (PDB ID: 3ERD8) complexes were used as the initial structures for the simulations. Crystal structure of the ERα-BPA complex is a mutant Y537S which is not suitable for this study. The structure of ERα-DPN complexes is not available in PDB. Thus, we used a molecular docking method to construct them by exploiting the Autodock 4.2 program.35 The 1ERE crystal structure was chosen as the template for this procedure. The Kollman charge was added to the protein. Autogrid was used to generate grid maps fixed in the E2 position. The box size was set at 60 Å in all three dimensions, and the spacing between the grid points was set at 0.375 Å. The Lamarckian Genetic Algorithm (LGA) search function was used and the parameter was set as follows. The number of GA runs equaled 100. The number of individuals in the population equaled 300. The maximum number of energy evaluations equaled 25000000. The maximum number of generations was 27000. And other parameters were set as default. After docking, all the conformations of ligand were analyzed using cluster analysis and energy ranking. The conformation with the lowest energy was chosen as the starting structure for the subsequent MD simulations.


Page 8 of 51

Page 9 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hydrogen atoms were added to the crystal structures using the tLeap module of AMBER program.36 All crystallographic water molecules were included in the MD simulations. All of the missing and incomplete residues in the four complexes were repaired by SYBYL-X 1.2 software package from Tripos, Inc. Co.. The amino acid protonation states were calculated by PROPKA37 and set as the calculated results. The force field parameters of the four ligands were supplied using the B3LYP/6-31G* theory level with the Gaussian 09 program38. The RESP fitting procedure of AMBER was used to derive the charges based on the optimized geometries. MD simulations All MD simulations were performed using AMBER with the ff12SB force field, which includes several updates and modifications.39 Each system was solvated in a box containing TIP3P water with at least a 10 Å distance around the complex. To maintain system neutrality, Na+ counterions were added. The solvated system was subjected to 2000 steps of steepest decent minimization followed by 3000 steps of conjugate gradient minimization. Subsequently, the system was heated from 0 to 310K for 300 ps using Langevin dynamics with a collision frequency of 1 ps-1. Then, the system was subjected to 500 ps equilibrium MD simulations for relaxation. Finally, a total of 40 ns were simulated for each system under NPT ensemble conditions without any restraint. Periodic boundary conditions and particle-mesh Ewald (PME) for long-range electrostatics were used in the simulations.40 The cut-off was 12 Å, and bonds involving hydrogen were held fixed using the SHAKE algorithm. The time step was set to 2 fs. A hydrogen bond was



identified if the distance between the donor and acceptor was shorter than 3.5 Å and the angle was less than 60°. All of the graphics were created with GaussView and Chimera.41, 42

Principal component analysis Principal component analysis (PCA) was used to separate the protein's conformational space into essential subspaces with enough degrees of freedom to describe the motions of the proteins. PCA uses a covariance matrix constructed with structures sampled from the MD simulations. By diagonalization of the covariance matrix of the system’s coordinates, the motions of the structure along the trajectory can be determined. To obtain the dominant motion over the MD simulation, the trajectory is projected along the direction described by a selected eigenvector. Subsequently, by calculating the two largest projections on the average structure, the primary direction of the protein can be discerned. PCA was performed using ProDy software.43 Three-dimensional structural snapshots were visualized using VMD and its plugin NMWiz.43, 44 MM-GBSA Calculations A total of 2500 snapshots were extracted from each equilibrated trajectory to perform MM-GBSA45 calculations. Using this method, the binding free energy is computed as follows: ΔGbind = Gcomplex – Gprotein – Gligand Gcomplex, Gprotein, and Gligand are the free energies of the complex, protein, and ligand, respectively. The free energy, G, can be calculated using the following equations:


Page 10 of 51

Page 11 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


G = EMM + Gsol – TS EMM = Eele + Evdw Gsol = GGB + GSA GSA = 0.005 × SASA EMM, Gsol, and TS represent the gas-phase molecular mechanics components, the stabilization energy caused by solvation, and a vibrational entropy term, respectively. EMM is given as the sum of Eele and Evdw, which are the Coulomb and van der Waals interaction terms, respectively. Gsol is the solvation energy and can be separated into an electrostatic solvation free energy (GGB) and a nonpolar solvation free energy (GSA). GGB can be calculated with the generalized Born (GB) method.46 In this study, GB model developed by Onufriev (GBOBCII) was used, which presented good result in calculating the electrostatic part of the solvation free energy.47 The radius used for the calculation of the polar solvation free energy was mbondi. The salt concentration was set to 0.1M. Generally, GSA is proportional to the molecular solvent-accessible surface area (SASA) that is buried on binding.48 Normal-mode analysis was performed to estimate the change in conformational entropy upon ligand binding (the TS term) using the nmode module of AMBER.49, 50 To obtain a detailed view of the ERα-ligand interactions, the MM-GBSA method was also used to calculate the binding free energy of each residue at the interface between the protein and ligand. The snapshots used in the binding free-energy decomposition are identical to those used in the binding free-energy calculations.



RESULTS Stability of the Global Structure The structures and calculated RESP distributions of E2, DES, BPA, and DPN are shown in Figure 1. As seen, DES, BPA and DPN share similar structures. The two phenol groups of DES, BPA and DPN are linked by 1-2 sp3 or sp2 carbon atoms. Likewise, despite structurally distinct, E2 presents a similar overall size to those of DES and DPN. In contrast, BPA only has one linking carbon atom between the two phenol groups and its size is relatively smaller than those of E2, DES and DPN. The hydroxyl groups of the four chemicals represent a strongly charged region. The hydrogen atoms in the hydroxyl groups are positively charged, whereas the oxygen atoms are negatively charged. Thus, the hydroxyl groups may form electrostatic interactions with the surrounding ERα residues. The other structural fragments of E2, DES and BPA are essentially neutral. In addition, the cyano group of DPN is negatively charged, which may affect the interactions in the DPN-ERα complex. After the simulations, the root-mean-square deviation (RMSD) was calculated to monitor system equilibrium. Figure S2 indicates that all four systems reached equilibrium very quickly. Moreover, all complexes remained stable during the simulations. The mean RMSD values were 2.12 Å (standard deviation 0.34 Å), 1.55 Å (0.17 Å), 2.02 Å (0.16 Å), and 1.60 Å (0.18 Å) for ERα-E2, ERα-DES, ERα-BPA, ERα-DPN, respectively. The ERα-DES complex has the smallest RMSD relative to its initial structure.


Page 12 of 51

Page 13 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We also performed cluster analysis on the four complexes to identify structural changes and to choose a representative structure for the subsequent calculations. Each trajectory is divided into five clusters (termed C1, C2, C3, C4 and C5, separately) using the average linkage algorithm. The RMSD values of C1-C5 corresponding to the initial structure are rather small (data not shown). These results indicate the limited fluctuation of the four structures over the course of simulations. The largest cluster of each trajectory was then chosen to be the representative structure of each complex during the MD simulations. We also compared the structure of largest cluster with the initial one (Figure S3). The global structures of the proteins are highly conserved, though the secondary structure and some loops fluctuate during the simulations, which would be discussed in the following sections. Additionally, the binding modes of the four ligands are quite different. Notably, the orientations of E2 and DES in the pocket are similar, and both agonists are essentially located at the initial position with similar poses. On the contrary, the positions and orientations of BPA and DPN change significantly relative to their initial configurations. BPA moves downwards into the cavity bottom of binding pocket and the conformation of DPN twists. To measure the stability of the hydrophobic pocket, the root-mean-square fluctuation (RMSF) of the backbone atoms was determined (Figure 2). The RMSF values of the pocket, including H3, H6-H8, H11, s1 and s2, are rather small in all tested complex systems. Even the important helix H12 in the DES, BPA and DPN complexes fluctuates very little. These results indicate that the hydrophobic pocket remains rigid during the



simulations. Residues of the loop between H9 and H10 are flexible in the ERα-E2 and ERα-DES systems, but they fluctuate slightly in the ERα-BPA and ERα-DPN complexes. Furthermore, H12 in the ERα-E2 complex exhibits higher fluctuation than those in other three complexes. Correlation Analysis of the Residues The above analysis indicates that the protein structure remained conserved when it interacted with different agonists. To highlight the correlations among residues, we performed cross-correlation analysis on the trajectories of the four complexes. As seen in Figure 3, each residue is positively correlated with the adjacent residues in all four complexes. However, the residue correlations for the other regions differ among the four complexes. For instance, H4-H5 (residues 365-385) is positively correlated with H9-H10 (residues 460-485) and with H1-H2 (residues 305-325) in the ERα-E2 complex, respectively. Additionally, H11-H12 (residues 505-550) is positively correlated with many nearby residues. Compared with the significant correlations observed in the ERα-E2 complex, many disappear in the ERα-DES complex. Positive correlations between H4-H5 (residues 365-385) and H9-H10 (residues 460-485) become uncorrelated. The ERα-BPA complex has more positive correlations than the ERα-E2 system in many regions. H6 (residues near 385) is positively correlated with H3 (residues near 350) and H8 (residues near 430), respectively. The β-sheet s1 (residues near 405) is highly positively correlated with H8 (residues near 425). For the ERα-DPN complex, many positively correlations in the ERα-E2 complex become negative ones. Rare positive


Page 14 of 51

Page 15 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


correlations are evident in ERα-LBD when it is bound to DPN. In addition, the portion of H3, locating near the binding pocket, has different correlation patterns with H12 in the four complexes. Primary Motions of the Four Complexes

PCA was used to obtain insight into the primary motions of ERα-LBD that could be affected by ligand binding. To this end, PCA was performed on the backbone atoms from all trajectories. The largest two eigenvalues were projected to indicate the distribution of the complex conformations (Figure S4). The conformation distributions of the four complexes are quite different, indicating their diverse motion patterns. Overall, the ERα-E2 complex has a broader conformational distribution than the other three systems. Essentially, its conformations mainly locate in the region with positive PC1 and negative PC2. The ERα-DES complex shows three distinct conformational clusters, which are distributed in both directions of PC1 and the positive direction of PC2. The ERα-BPA complex has a major conformational cluster that locates in the positive direction of PC1. The conformational distribution of the ERα-DPN complex is the most concentrated one and primarily locates in the negative direction of PC1 and the positive direction of PC2. Figure 4 further visualizes the motions presented by PC1, the largest PC mode, to monitor the ERα mobility. For the ERα-E2 complex, the loop connecting H9 and H10 is a flexible region, as is the H12 terminal. In comparison with that in the ERα-E2 complex (23 Å), less mobile H9-H10 loop (12 Å) is observed in the ERα-DES complex.(Figure S5) The mobility of the ERα-BPA system differs from that of the ERα complexes with E2



and DES. For example, in the ERα-BPA system, the H9-H10 loop is rigid while H1, the H2-H3 loop, H7 and H8 present mobility. For the ERα-DPN complex, the H1-H2 region is mobile, whereas other regions moved little. Overall, the PC1 defined motion mode corresponds to the breathing movement of ERα-LBD, and demonstrates how the receptor responds uniquely after binding different ligands. Different receptor motility leads to the corresponding change in the surface electrostatic potential.(Figure 5) In the structure of the ERα, the AF-2 region is generally negatively charged. After DES binding, the movement of the H9-H10 loop helps the exposure of the nearby negatively charged region, thereby increasing the negatively charged surface of the AF-2 region. This is very similar to the endogenous hormone E2. The situation is totally different for BPA and DPN. The area of the negatively charged surface decreased remarkably in the AF-2 region of the ERα-BPA and ERα-DPN complexes instead, which can be attributed to the weak mobility of the H9-H10 loop. Inter-molecular Interaction Analysis Previous research has shown that the ligands bind to ERα via a hydrophobic pocket. The CASTp program51 calculated the volume of the cavity to be 594 Å3, which was sufficient to accommodate the small ligands such as E2, DES, BPA and DPN. To investigate the hydrophobic interactions between the ligands and ERα, we illustrated the hydrophobic environment around each ligand as shown in Figure S6. Almost the entire pocket is of high hydrophobicity, whereas the region near Glu353 and Arg394 is hydrophilic. DES, DPN and E2 occupy the entire hydrophobic pocket, while


Page 16 of 51

Page 17 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the small-size BPA only partially occupies the pocket. Clearly, the nonpolar structural fragments of the ligands play an important role in ligand-receptor interaction. Among the four ligands, E2 has the most hydrophobic skeleton, steroid nucleus, while BPA had the least one. E2, DES and BPA all have two polar hydroxyl groups at each end, and the remaining portions of the molecules are nonpolar. One polar hydroxyl group in E2 or DES is located near the hydrophilic region of the pocket and the other in the hydrophobic region. But both hydroxyl groups in BPA locate in an area of less hydrophilicity. DPN has three polar groups, including two hydroxyl groups at each end and one cyano group at the center. Although the two hydroxyl groups and the nonpolar skeleton of DPN locate in the same environment as those of E2 and DES, a hydrophobic environment surrounds the polar cyano group. For this reason, DES, E2 and BPA were more compatible with the hydrophobic pocket than DPN. The formation of π-π interactions between ERα and each ligand was also examined and is shown in Figure 6. DES presents a very strong π-π interaction with Phe404. The phenyl ring of Phe404 is located just above the phenyl ring in DES. The two benzene rings are perpendicular to each other, which is the strongest type of π-π interactions.52 E2 is found to align in a similar manner as DES does inside the ERα-LBD. It also processes strong π-π interactions to Phe404. However, the phenyl groups of BPA and DPN exhibit weak π-π interactions with ERα-LBD. Structurally, BPA is shorter than DES, which creates a longer distance between its phenyl group and the key residue Phe404. As for DPN, there is repulsion between the hydrophobic pocket and its cyano group. The distances between



the phenyl group of each ligand and the phenyl group of Phe404 were also examined. According with the criterion of π-π interaction52, E2 and DES are approximately 5 Å away from Phe404 (Figure S7). DPN is larger than E2 and DES for about 0.4 Å. BPA is about 6 Å away from Phe404, which suggests weaker or non-existent π-π interactions. Hydrogen bonds (HBs) also play very important roles in the ERα-ligand interactions. Both phenolic hydroxyl groups of DES form HBs with the key residues Glu353 and Glu521 in ERα (Figure 7). Figure 8 clearly shows that both HBs of DES are conserved during the MD simulations. Three HBs involving E2, His524, Gly521 and Glu353 are also evident. One cyclopentanol group of E2 forms two HBs with His524 and Gly521, both are highly conserved during the MD simulations. The phenolic hydroxyl group of E2 unstably forms one hydrogen bond with Glu353, the cleavage of which occurs after approximately 7 ns. Thus, only one side of E2 establishes HBs with ERα in the stable ERα-E2 complex. BPA forms only one hydrogen bond with Thr347, which is not stable at the end of simulation. Notably, this residue is different from the residues participating in HBs formation in the ERα interaction with E2 or DES, which is decided by the different orientation of BPA within the binding pocket. DPN does not form a hydrogen bond with ERα. Binding Free Energy and the Driving Force for Binding To investigate the driving force of the ligand binding, binding free energies were computed using the MM-GBSA method. To minimize any potential bias from the initial structure, the first 30 ns of each simulation was excluded from the energy calculations.


Page 18 of 51

Page 19 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The overall binding free energy is divided into contributions from electrostatic interactions, van der Waals interactions, electrostatic solvation interactions, and nonpolar solvation interactions. The contribution of entropy was also considered. For all computed energies, the results are presented as ensemble-averaged values over 2500 snapshots. Table 1 lists the binding free energies of the four complexes. DES presents the highest binding affinity with ERα, followed by E2, DPN and BPA, in turn. These results are consistent with the experimental datas.20,

21

We also calculated the relative binding

affinity (RBA) in silico to establish a relationship between the RBA in silico and the RBA in vitro. A statistically significant correlation is obtained. (R2=0.99, p E2 (-1.81 kcal/mol) > DPN (-1.58 kcal/mol) > BPA (-1.13 kcal/mol). The Phe404 contribution to the binding affinity is larger for DES (-1.58 kcal/mol) and E2 (-1.80 kcal/mol) than that for BPA (-1.10 kcal/mol) and DPN (-1.25 kcal/mol). This result



illustrates the differences in both ligand binding sites and binding poses among the four complexes. DISCUSSION MD simulations, binding free energy calculations and PCA were adopted to investigate the underlying basis of the different ERα allosteric effect induced by structurally similar ligands, which are endogenous estrogen E2 and exogenous EEs including DES, BPA and DPN. The study provides many further insights into the essence of the ERα-EEs interactions and subsequent diverse regulatory effects. The Chemicals Present Estrogen-Like Activity All the ligands bind to ERα in the binding pocket stably with negative binding free energies, which shows that they are all agonists for ERα. The results from MM-GBSA calculations demonstrate that the major stability of the ligand-receptor complexes is attributed to the nonpolar interactions (mainly van der Waals interactions) between the protein and the ligands. Clearly, the hydrophobic pocket binds hydrophobic molecules with high affinity. The hydrophilic groups in the ligand structure interact unfavorably with the hydrophobic pocket, but polar residues that can establish stable HBs with them assist in the complex formation. All of four ligands have two hydroxyl groups (one on each end), and the potential repulsion of them by the binding pocket is indicated by the unfavorable electrostatic solvation free energy (EGB). Another feature of the pocket is its rigidity. The structural analysis indicates that the pocket-forming residues tent to maintain their positions during the MD simulations, and the ligand binding fails to induce


Page 22 of 51

Page 23 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


clear structural changes in the pocket. The chemicals with the ERα agonism must be accommodated in this relatively rigid pocket. MM-GBSA also helps us identify the key residues for the ERα-ligand interactions. The residues of the binding pocket form a hydrophobic network that interacts strongly with the hydrophobic groups of the ligands. Most of them are leucine, while phenylalanine offers hydrophobic binding to all ligands and can also establish π-π interactions with E2 and DES. Leu346, Leu387, Leu391, Phe404 and Ala350 play significant roles in receptor binding for the four ligands. Additionally, several polar residues located on the either side of the binding pocket are potential participants in forming conserved hydrogen bonds, such as Glu353, Gly521 and His524. These residues provide favorable binding free energies with the ligands that can fully occupy the pocket. We also established the relationship between the in silico and in vitro RBA. Generally, the MM-GBSA method is not able to reproduce the absolute value of the experimental results precisely. However, there is an excellent linear correlation between the calculated RBA and the experimental results with a significance level lower than 0.05. Our aim is to use MM-GBSA calculations to explore the sources of binding affinity rather than to calculate absolute binding free energies. Thus, we are more concerned with identifying the contribution of various interaction types and key residues responsible for ligand binding and transcription activation. In this study, MM-GBSA is proven to be a reliable approach to investigate ERα-ligand interactions, and it may be effectively extended to



calculate the binding free energies of other EDCs upon receptor-mediated effect to understand their molecular toxicological mechanisms. Similar chemicals with different binding affinities The hydrophobicity of the binding pocket in ERα determines that the number and categories of hydrophobic and hydrophilic groups in each chemical can significantly affect its ERα binding affinity. Through a structural comparison of E2, DES, BPA and DPN, we have observed that the number of hydrophobic groups decreases in the order of E2 > DES > BPA, DPN, which is consistent with the descending order of their van der Waals interaction energies. A dramatically small electrostatic interaction energy is found for DPN, which, relative to the other three ligands, has an additional hydrophilic cyano group repelled by the hydrophobic environment. This finding may be one of the reasons why DPN has low binding affinity for ERα.53 We can conclude that the agonists with high binding affinity for ERα should avoid hydrophilic groups along their backbone interior. The rigidity of the binding pocket determines that ERα tends to form stable inter-molecular interactions with the ligands of proper size. Hydrogen bonds play very important roles in the ligand-receptor interactions. E2, DES can form HBs with ERα-LBD at both ends. In contrast, BPA forms just one hydrogen bond at one end of the molecule because its relatively smaller size makes the establishment of an additional hydrogen bond at the other end impossible. Only the ligands completely occupying the ERα pocket possess HBs with the receptor at both ends. Moreover, the stabilities of the


Page 24 of 51

Page 25 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


HBs are rather different for the 3 chemicals (Figure 7 and 8). The phenolic hydroxyl groups at the either side of DES are able to establish conserved HBs with Glu353 and Gly521, respectively. For E2, the cyclopentanol group forms two HBs with ERα-LBD, via Gly521 and His524, but the HB between Glu353 and the phenolic hydroxyl group on the other side cannot be sustained during the simulation. In addition, the key residues responsible for hydrogen bonding, are differentiated for the good ligands DES and E2. The hydrogen bond between E2 and His524 is not observed in the DES-ERα complex, while His524 is thought to be essential for some normal regulation effects of ERα, indicating the different regulatory effects of DES relative to E2.54 Figure 1 shows that DES, BPA and DPN are relatively flexible molecules when the rigid steroid ring in E2 is considered. High ligand plasticity helps to initiate an effective induce-fit process to make the chemical adopt proper conformation and pose to fit in the rigid pocket of ERα in this regard. In contrast, it is difficult for the structurally rigid E2 to fit in the pocket and results in the unstable hydrogen bonding at the rigid phenyl side of E2. Obviously, DPN is less compatible with the ERα-LBD because of its cyano group, repelled by the hydrophobic pocket. This repulsion causes a break of HB in both ends of DPN. The MM-GBSA calculations support the above analysis and confirm our reasoning (Table 1). DES, which can form HBs on both sides of the ligands with the ERα-LBD, have inter-molecular electrostatic interaction energies of approximately -21 kcal/mol. As expected, E2 and BPA, establishing just one stable HB with the ERα-LBD, have inter-molecular electrostatic interaction energies of approximately -12 kcal/mol. DPN,



which cannot form HB with ERα-LBD, only presents the inter-molecular electrostatic interaction energy of approximately -3 kcal/mol. As for another kind of notable inter-molecular interaction, π-π interaction between ERα and the ligands, only E2 and DES can form well-defined π-π interactions with Phe404, respectively. The small size of BPA makes it difficult to form HBs and π-π interactions simultaneously. Although DPN occupies the binding pocket as completely as DES, its hydrophilic cyano group is exposed to the hydrophobic pocket. The resulting repulsion negatively influences the position and orientation of DPN in the pocket, pushing it away from Phe404. Compared with DES, BPA and DPN have weak π-π interactions with Phe404. For this reason, BPA and DPN, in comparison with E2 and DES, have diminished hydrophobic interactions with the binding pocket, further verified by their corresponding van der Waals energies calculated by the MM-GBSA method. In conclusion, as a result of the rigid cavity in ERα-LBD, only ligands with proper size, functional groups and flexibility can bind tightly with ERα. Different binding modes induce diverse regulatory effects Previous research has shown that E2, DES, BPA and DPN are agonists for ERα, though each ligand affects transcription regulation differently. This phenomenon is difficult to explain based on the current experimental data. Many important residues for ligand binding of ERα are different for the four agonists. Table 2 lists all the 16 important residues. Clearly, each ligand only interacts with some of these residues for complex formation. Although some residues play important roles in the stability of two or more


Page 26 of 51

Page 27 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


complexes, their contributions differ. The hydrogen bonds involved in the ligand-receptor interaction are also different for the four complexes. His524 has a stable hydrogen bond with E2, which cannot be formed with DES, BPA and DPN. The His524-involved HB is shown to be critical for the normal regulation of ERα-mediated transcription. Its disruption has obvious effects on gene expression.54 Overall, these results demonstrate the diverse binding modes of ERα with the different ligands. Although the entire protein is relatively rigid, the different binding modes still influence the orientation of the inner residues of the pocket and thus change their relations with external residues. More precisely, even though the four ERα complexes all possess a similar agonism conformation with a conserved H12 position, the ligand binding mode may affect the exposed surface of LBD associated AF-2, known as the coactivator binding site. Subtle modulation of ER-ligand interactions has been proven to be sufficient to induce highly divergent regulation results. Resveratrol was ever thought to act as a superagonist and provoke a super-induction effect initially, but many subsequent reports indicate that it is a partial agonist with non-proliferative activity for ER.55, 56 Nwachukwu et al. showed that the adoption of multiple orientations/poses in the binding pocket, which dramatically influence the exposed surface of AF-2 region, is the main reason for the regulatory diversity of resveratrol.57 In this study, the four complexes show obviously different correlation patterns between residues in H4-H5 and residues in H9-H10, as well as H3 to H12, and thus present how the subtle allosteric effect of the pocket affects the AF-2 region. PCA provides further proofs. The motions of ERα-LBD, especially the H9-H10



loop, adjacent to AF-2, vary when different ligands bind to the pocket of ERα. AF-2 is the binding site of cofactors and plays a very important role in determining the final regulatory effect. The negative electrostatic potential of the AF-2 region surface has been proven to be critical to the normal function of ERα.58 Regarding the ERα-DES interaction, the high mobility of H9-H10 loop increase the negative charged surface area of the AF-2 region, which leads to the readily recruitment of diverse co-activators and gene overexpression.59 The different AF-2 surface electrostatic potential induced by the H9-H10 loop motion in the ERα-E2 and ERα-DES complexes seems correlated with the different aggressive behaviors.23 However, the relatively decreased mobility of the H9-H10 loop in the ERα-BPA and ERα-DPN complexes leads to insufficient exposure of the coactivator binding site, and the following coactivator recruitment via protein-protein interaction becomes less active relative to that for E2 and DES. Correspondingly, in comparison with E2 and DES, BPA and DPN fail to cause significant changes in transcription activation and gene expression, which is indicative of their inability to recruit co-activators.60 In addition, DES has been reported to promote hypermethylation of the homeobox gene, Hoxa10, while exposure to BPA could result in hypomethylation of the Hoxa10 promoter.61 Such high selectivity in regulation effects of the different ligands can be partly attributed to both differentiated area and feature of the exposed AF-2 surface with negative charged potential in the ERα-DES and ERα-BPA complexes. Another possible explanation is that the mobility diversity of the H9-H10 loop also influence subsequent ERα dimerization since it lying at the interface of ERα dimerization.


Page 28 of 51

Page 29 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Moreover, it should be noted there are more regions affected by the ligand binding. The H2-H3 loop exhibits significant mobility in the ERα complexes with DES, BPA and DPN, but its conformation is relatively stable in the ERα-E2 complex. Since there is no complete structure of ERα, it is difficult to directly determine the function of this region. The detailed effects of these regions on regulation pathways require further investigation. CONCLUSION MD simulations combined with PCA and MM-GBSA calculations were performed to investigate ligand-specific responses of three similar chemicals and an endogenous estrogen. The results indicate that small structural changes in the ligands can cause distinct different in binding affinities and regulation effects on ERα. All four ligands bind within a rigid binding pocket, which is incompatible with hydrophilic ligands. Van der Waals interaction between the ERα pocket and the ligand is the main driving force to the ligand binding and it is directly influenced by the number of hydrophobic groups in a ligand. Hydrogen bonds established between the ligands and ERα-LBD also play very important roles in stabilizing the complexes. Among the four chemicals, DES with proper size and suitable flexibility is the most stable ligand within the rigid pocket of ERα. The MM-GBSA results indicate that Leu346, Leu387, Leu391, Phe404 and Ala350 are essential residues in ligand binding for all four chemicals. Other differentiated residues contribute to the distinct binding affinities and poses of the four ligands, and further influence the correlations among residues. Furthermore, different binding modes result in unique mobility for each complex, particularly in the region of the H9-H10 loop. The



distinct mobility of this region dramatically changes the electrostatic potential surface of AF-2, which in turn clearly affect the coactivator recruitment and the subsequent gene expression. Interestingly, high mobility of the H9-H10 loop induced by DES is correlated with the significant up-regulation of specific gene expression while its low mobility induced by BPA and DPN is correlated with the little up-regulation effect. Understanding how the binding mechanisms affect estrogen signaling pathways will enrich our knowledge of the effect of EEs on human health. ASSOCIATED CONTENT Supporting Information Representative structure of ERα-LBD; Time-dependent RMS deviation; Cluster analysis; Projection of PCA; Mobility of residues; Surface of the binding pocket; The distances between the aromatic rings of Phe404 and the ligands; the relationship of RBA in silico versus RBA in vitro; Binding free energy of each residue. AUTHOR INFORMATION Corresponding Author *State Key Laboratory of Environmental Chemistry and Ecotoxicology, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing 100085, China. E-mail: [email protected] Funding Information


Page 30 of 51

Page 31 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This study was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB14030500), and the National Natural Science Foundation (91743204, 21507152, 21607168, 21777181) . Notes The authors declare no competing financial interests. REFERENCES (1)

Braun, J. M. (2017) Early-life exposure to EDCs: role in childhood obesity and

neurodevelopment. Nat. Rev. Endocrinol. 13, 161-173. (2)

Gore, A. C., Chappell, V. A., Fenton, S. E., Flaws, J. A., Nadal, A., Prins, G. S.,

Toppari, J., and Zoeller, R. T. (2015) EDC-2: The endocrine society's second scientific statement on endocrine-disrupting chemicals. Endocr. Rev. 36, E1-E150. (3)

Adeel, M., Song, X. M., Wang, Y. Y., Francis, D., and Yang, Y. S. (2017)

Environmental impact of estrogens on human, animal and plant life: A critical review. Environ. Int. 99, 107-119. (4)

Vom Saal, F. S., Timms, B. G., Montano, M. M., Palanza, P., Thayer, K. A.,

Nagel, S. C., Dhar, M. D., Ganjam, V., Parmigiani, S., and Welshons, W. V. (1997) Prostate enlargement in mice due to fetal exposure to low doses of estradiol or diethylstilbestrol and opposite effects at high doses. Proc. Natl. Acad. Sci. 94, 2056-2061.



(5)

Usman, A., and Ahmad, M. (2016) From BPA to its analogues: Is it a safe

journey? Chemosphere 158, 131-142. (6)

Lejonklou, M. H., Dunder, L., Bladin, E., Pettersson, V., Ronn, M., Lind, L.,

Walden, T. B., and Lind, P. M. (2017) Effects of low-dose developmental bisphenol A exposure on metabolic parameters and gene expression in male and female fischer 344 rat offspring. Environ. Health. Perspect. 125. (7)

Rahman, M. S., Kwon, W.-S., Karmakar, P. C., Yoon, S.-J., Ryu, B.-Y., and

Pang, M.-G. (2017) Gestational exposure to bisphenol A affects the function and proteome profile of F1 spermatozoa in adult mice. Environ. Health. Perspect. 125, 238-245. (8)

Shiau, A. K., Barstad, D., Loria, P. M., Cheng, L., Kushner, P. J., Agard, D. A.,

and Greene, G. L. (1998) The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 95, 927-937. (9)

Delfosse, V., Grimaldi, M., Pons, J.-L., Boulahtouf, A., le Maire, A., Cavailles,

V., Labesse, G., Bourguet, W., and Balaguer, P. (2012) Structural and mechanistic insights into bisphenols action provide guidelines for risk assessment and discovery of bisphenol A substitutes. Proc. Natl. Acad. Sci. 109, 14930-14935. (10)

Enmark, E., Pelto-Huikko, M., Grandien, K., Lagercrantz, S., Lagercrantz, J.,

Fried, G., Nordenskjöld, M., and Gustafsson, J.-A. k. (1997) Human estrogen receptor


Page 32 of 51

Page 33 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


β-gene structure, chromosomal localization, and expression pattern 1. J. Clin. Endocrinol. Metab. 82, 4258-4265. (11)

Hall, J. M., and McDonnell, D. P. (1999) The estrogen receptor β-isoform (ERβ)

of the human estrogen receptor modulates ERα transcriptional activity and is a key regulator of the cellular response to estrogens and antiestrogens 1. Endocrinology 140, 5566-5578. (12)

Brzozowski, A. M., Pike, A. C., Dauter, Z., Hubbard, R. E., Bonn, T., Engström,

O., Öhman, L., Greene, G. L., Gustafsson, J.-Å., and Carlquist, M. (1997) Molecular basis of agonism and antagonism in the oestrogen receptor. Nature 389, 753-758. (13)

Nwachukwu, J. C., Srinivasan, S., Bruno, N. E., Nowak, J., Wright, N. J.,

Minutolo, F., Rangarajan, E. S., Izard, T., Yao, X. Q., and Grant, B. J. (2016) Systems structural biology analysis of ligand effects on ERα predicts cellular response to environmental estrogens and anti-hormone therapies. Cell. Chem. Biol. 24, 35. (14)

Delfosse, V., Grimaldi, M., Cavaillès, V., Balaguer, P., and Bourguet, W. (2014)

Structural and functional profiling of environmental ligands for estrogen receptors. Environ. Health Perspect. 122, 1306-1313. (15)

Shang, Y. F., and Brown, M. (2002) Molecular determinants for the tissue

specificity of SERMs. Science 295, 2465-2468.



(16)

Smith, C. L., and O'Malley, B. W. (2004) Coregulator function: A key to

understanding tissue specificity of selective receptor modulators. Endocr. Rev. 25, 45-71. (17)

Arase, S., Ishii, K., Igarashi, K., Aisaki, K., Yoshio, Y., Matsushima, A.,

Shimohigashi, Y., Arima, K., Kanno, J., and Sugimura, Y. (2011) Endocrine Disrupter Bisphenol A Increases In Situ Estrogen Production in the Mouse Urogenital Sinus. Biol. Reprod. 84, 734-742. (18)

Benmansour, S., Adeniji, O. S., Privratsky, A. A., and Frazer, A. (2016) Effects

of long-term treatment with estradiol and estrogen receptor subtype agonists on serotonergic function in ovariectomized rats. Neuroendocrinology 103, 269-281. (19)

Song, X. Z., and Pan, Z. Z. (2012) Estrogen receptor-beta agonist

diarylpropionitrile counteracts the estrogenic activity of estrogen receptor-alpha agonist propylpyrazole-triol in the mammary gland of ovariectomized Sprague Dawley rats. J. Steroid Biochem. Mol. Biol. 130, 26-35. (20)

Kuiper, G. G., Lemmen, J. G., Carlsson, B., Corton, J. C., Safe, S. H., van der

Saag, P. T., van der Burg, B., and Gustafsson, J.-A. k. (1998) Interaction of estrogenic chemicals and phytoestrogens with estrogen receptor β. Endocrinology 139, 4252-4263. (21)

Shanle, E. K., and Xu, W. (2010) Selectively targeting estrogen receptors for

cancer treatment. Adv. Drug Delivery Rev. 62, 1265-1276.


Page 34 of 51

Page 35 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22)

Delfosse, V., Grimaldi, M., Pons, J. L., Boulahtouf, A., le Maire, A., Cavailles,

V., Labesse, G., Bourguet, W., and Balaguer, P. (2012) Structural and mechanistic insights into bisphenols action provide guidelines for risk assessment and discovery of bisphenol A substitutes. Proc. Natl. Acad. Sci. 109, 14930-14935. (23)

Tanaka, M., Ohtani-Kaneko, R., Yokosuka, M., and Watanabe, C. (2004)

Low-dose perinatal diethylstilbestrol exposure affected behaviors and hypothalamic estrogen receptor-alpha-positive cells in the mouse. Neurotoxicol. Teratol. 26, 261-269. (24)

Farabollini, F., Porrini, S., Seta, D. D., Bianchi, F., and Dessifulgheri, F. (2002)

Effects of perinatal exposure to bisphenol A on sociosexual behavior of female and male rats. Environ. Health. Perspect. 110, 409-414. (25)

Lee, H. R., Hwang, K. A., Park, M. A., Yi, B. R., Jeung, E. B., and Choi, K. C.

(2012) Treatment with bisphenol A and methoxychlor results in the growth of human breast cancer cells and alteration of the expression of cell cycle-related genes, cyclin D1 and p21, via an estrogen receptor-dependent signaling pathway. Int. J. Mol. Med. 29, 883-890. (26)

Nakamura, T., Katsu, Y., Watanabe, H., and Iguchi, T. (2008) Estrogen receptor

subtypes selectively mediate female mouse reproductive abnormalities induced by neonatal exposure to estrogenic chemicals. Toxicology 253, 117-124.



(27)

Zhang, J., Li, Y., Gupta, A. A., Nam, K., and Andersson, P. L. (2016)

Identification and Molecular Interaction Studies of Thyroid Hormone Receptor Disruptors among Household Dust Contaminants. Chem. Res. Toxicol. 29, 1345-1354. (28)

Li, Y., Perera, L., Coons, L. A., Burns, K. A., Tyler, R. J., Pelch, K. E.,

Houtman, R., Van, B. R., Teng, C. T., and Korach, K. S. (2018) Differential in vitro biological action, coregulator interactions, and molecular dynamic analysis of bisphenol A (BPA), BPAF, and BPS ligand-ERα complexes. Environ. Health. Perspect. 126, 017012. (29)

Ruiz-Rodriguez, M. A., Vedani, A., Flores-Mireles, A. L., Chairez-Ramirez, M.

H., Gallegos-Infante, J. A., and Gonzalez-Laredo, R. F. (2017) In Silico Prediction of the Toxic Potential of Lupeol. Chem. Res. Toxicol.30, 1562-1571. (30)

Xue, Q., Zheng, Q. C., Zhang, J. L., Cui, Y. L., and Zhang, H. X. (2014)

Exploring the mechanism how Marburg virus VP35 recognizes and binds dsRNA by molecular dynamics simulations and free energy calculations. Biopolymers 101, 849-860. (31)

Xue, Q., Zhang, J.-L., Zheng, Q.-C., Cui, Y.-L., Chen, L., Chu, W.-T., and

Zhang, H.-X. (2013) Exploring the molecular basis of dsRNA recognition by Mss116p using molecular dynamics simulations and free-energy calculations. Langmuir 29, 11135-11144.


Page 36 of 51

Page 37 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(32)

Zhuang, S. L., Zhang, J., Wen, Y. Z., Zhang, C. L., and Liu, W. P. (2012)

Distinct mechanisms of endocrine disruption of DDT-related pesticides toward estrogen receptor alpha and estrogen-related receptor gamma. Environ. Toxicol. Chem. 31, 2597-2605. (33)

Lu, L. P., Zhan, T. J., Ma, M., Xu, C., Wang, J. P., Zhang, C. L., Liu, W. P., and

Zhuang, S. L. (2018) Thyroid disruption by Bisphenol S analogues via thyroid hormone receptor beta: in Vitro, in Vivo, and molecular dynamics simulation study. Environ. Sci. Technol. 52, 6617-6625. (34)

Zhuang, S. L., Zhang, C. L., and Liu, W. P. (2014) Atomic insights into distinct

hormonal activities of bisphenol A analogues toward PPAR gamma and ER alpha receptors. Chem. Res. Toxicol. 27, 1769-1779. (35)

Morris, G. M., Huey, R., Lindstrom, W., Sanner, M. F., Belew, R. K., Goodsell,

D. S., and Olson, A. J. (2009) AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 30, 2785-2791. (36)

Case, D., Darden, T., Cheatham III, T., Simmerling, C., Wang, J., Duke, R., Luo,

R., Walker, R., Zhang, W., and Merz, K. (2012) AMBER 12. University of California, San Francisco 1, 3.



(37)

Olsson, M. H., Søndergaard, C. R., Rostkowski, M., and Jensen, J. H. (2011)

PROPKA3: consistent treatment of internal and surface residues in empirical p K a predictions. J. Chem. Theory Comput. 7, 525-537. (38)

Frisch, A. (2009) Gaussian 09: User's Reference. Gaussian.

(39)

Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C.

(2006) Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins. 65, 712-725. (40)

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L.

G. (1995) A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577-8593. (41)

Dennington, R., Keith, T., and Millam, J. (2009) GaussView, version 5.

Semichem Inc., Shawnee Mission KS. (42)

Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M.,

Meng, E. C., and Ferrin, T. E. (2004) UCSF Chimera—a visualization system for exploratory research and analysis. J. Comput. Chem. 25, 1605-1612. (43)

Bakan, A., Meireles, L. M., and Bahar, I. (2011) ProDy: Protein dynamics

inferred from theory and experiments. Bioinformatics 27, 1575-1577. (44)

Humphrey, W., Dalke, A., and Schulten, K. (1996) VMD: Visual molecular

dynamics. J. Mol. Graph. Model. 14, 33-38.


Page 38 of 51

Page 39 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(45)

Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., Lee, M.,

Lee, T., Duan, Y., and Wang, W. (2000) Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889-897. (46)

Hawkins, G. D., Cramer, C. J., and Truhlar, D. G. (1996) Parametrized models

of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem. 100, 19824-19839. (47)

Feig, M., Onufriev, A., Lee, M. S., Im, W., Case, D. A., and Brooks, C. L.

(2004) Performance comparison of generalized born and Poisson methods in the calculation of electrostatic solvation energies for protein structures. J. Comput. Chem. 25, 265-284. (48)

Weiser, J., Shenkin, P. S., and Still, W. C. (1999) Approximate atomic surfaces

from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 20, 217-230. (49)

Hawkins, G. D., Cramer, C. J., and Truhlar, D. G. (1995) Pairwise solute

descreening of solute charges from a dielectric medium. Chem. Phys. Lett. 246, 122-129. (50)

Xu, B., Shen, H., Zhu, X., and Li, G. (2011) Fast and accurate computation

schemes for evaluating vibrational entropy of proteins. J. Comput. Chem. 32, 3188-3193. (51)

Dundas, J., Ouyang, Z., Tseng, J., Binkowski, A., Turpaz, Y., and Liang, J.

(2006) CASTp: computed atlas of surface topography of proteins with structural and



topographical mapping of functionally annotated residues. Nucleic. Acids. Res. 34, W116-W118. (52)

Tsuzuki, S., Honda, K., Uchimaru, T., Mikami, M., and Tanabe, K. (2002)

Origin of attraction and directionality of the π/π interaction: model chemistry calculations of benzene dimer interaction. J. Am. Chem. Soc. 124, 104-112. (53)

Sun, J., Baudry, J., Katzenellenbogen, J. A., and Katzenellenbogen, B. S. (2008)

Molecular basis for the subtype discrimination of the estrogen receptor-β-selective ligand, diarylpropionitrile. Mol. Endocrinol. 17, 247-258. (54)

Gao, L., Tu, Y., Agren, H., and Eriksson, L. A. (2012) Characterization of

agonist binding to His524 in the estrogen receptor α ligand binding domain. J. Phys. Chem. B 116, 4823-4830. (55)

Gehm, B. D., Mcandrews, J. M., Chien, P., and Jameson, J. L. (1997)

Resveratrol, a polyphenolic compound found in grapes and wine, is an agonist for the estrogen receptor. Proc. Natl. Acad. Sci. 94, 14138-14143. (56)

Bowers, J. L., Tyulmenkov, V. V., Jernigan, S. C., and Klinge, C. M. (2000)

Resveratrol acts as a mixed agonist/antagonist for estrogen receptors α and β. Endocrinology 141, 3657-3667. (57)

Nwachukwu, J. C., Srinivasan, S., Bruno, N. E., Parent, A. A., Hughes, T. S.,

Pollock, J. A., Gjyshi, O., Cavett, V., Nowak, J., and Garciaordonez, R. D. (2014)


Page 40 of 51

Page 41 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Resveratrol modulates the inflammatory response via an estrogen receptor-signal integration network. eLife 3. (58)

Jeselsohn, R., Yelensky, R., Buchwalter, G., Frampton, G., Meric-Bernstam, F.,

Gonzalez-Angulo, A. M., Ferrer-Lozano, J., Perez-Fidalgo, J. A., Cristofanilli, M., Gomez, H., Arteaga, C. L., Giltnane, J., Balko, J. M., Cronin, M. T., Jarosz, M., Sun, J., Hawryluk, M., Lipson, D., Otto, G., Ross, J. S., Dvir, A., Soussan-Gutman, L., Wolf, I., Rubinek, T., Gilmore, L., Schnitt, S., Come, S. E., Pusztai, L., Stephens, P., Brown, M., and Miller, V. A. (2014) Emergence of constitutively active estrogen receptor-alpha mutations in pretreated advanced estrogen receptor-positive breast cancer. Clin. Cancer. Res. 20, 1757-1767. (59)

Nair, R., and Shaha, C. (2003) Diethylstilbestrol induces rat spermatogenic cell

apoptosis in vivo through increased expression of spermatogenic cell Fas/FasL system. J. Biol. Chem. 278, 6470-6481. (60)

Gertz, J., Reddy, T. E., Varley, K. E., Garabedian, M. J., and Myers, R. M.

(2012) Genistein and bisphenol A exposure cause estrogen receptor 1 to bind thousands of sites in a cell type-specific manner. Genome. Res. 22, 2153-2162. (61)

Gibson, D. A., and Saunders, P. T. K. (2014) Endocrine disruption of oestrogen

action and female reproductive tract cancers. Endocr. Relat. Cancer 21.



Figure 1. Top: structures of E2, DES, BPA and DPN. Bottom: The electrostatic potential surface of E2, DES, BPA and DPN. The electron density ranges from -0.06 a.u. to 0.06 a.u. The negatively charged region is in red and the positively charged region is in blue.


Page 42 of 51

Page 43 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. RMSF of the backbone atoms relative to the initial structures for the ERα-E2, ERα-DES, ERα-BPA and ERα-DPN complexes.



Figure 3. Cross-correlation matrices of the coordinate fluctuations for Cα atoms of ERα-E2, ERα-DES, ERα-BPA and ERα-DPN complexes around their mean positions in the 40 ns simulations. Positive correlations in descending order are shown in red, orange, and yellow. Negative correlations are shown in blue. The main differences among the four complexes are labeled in red box.


Page 44 of 51

Page 45 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Porcupine representation of the PC1 mode of the ERα-E2, ERα-DES, ERα-BPA, and ERα-DPN complexes. Arrows in cyan show the direction of motion along PC1. Mobility ranges from high (blue) to low (red).



Figure 5. Electrostatic potential surfaces of the ERα-E2, ERα-DES, ERα-BPA, and ERα-DPN complexes. The positive area is colored in blue. The negative area is colored in red.


Page 46 of 51

Page 47 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. π-π Interactions for Phe404-E2, Phe404-DES, Phe404-BPA and Phe404-DPN. Phe404 is colored in pink, and each ligand is colored in orange. The aromatic ring is also shown.



Figure 7. Comparison of hydrogen bonds involving E2, DES, BPA and DPN with the receptor residues inside the binding pocket.

Figure 8. Hydrogen bond occupancies from the MD simulations for the ligand-receptor complexes.


Page 48 of 51

Page 49 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9. The superimposition of the representative conformation of each complex (E2 in orange, DES in yellow, BPA in green, and DPN in blue). Distributions of important residues (< -1 kcal/mol) around E2, DES, BPA and DPN are shown in the four side windows.



Page 50 of 51

Table 1. Binding free energy and its components for the ERα-E2, ERα-DES, ERα-BPA, ERα-DPN complexes (kcal/mol) Iterms

E2

DES

BPA

DPN

ΔEvdw

-43.91

-39.70

-35.22

-35.11

ΔEele

-12.41

-21.40

-11.36

-3.67

ΔGGB

23.89

24.01

21.36

12.58

ΔGSA

-5.27

-5.78

-4.87

-5.18

ΔHa

-37.71

-42.87

-30.09

-31.38

TΔS

-21.17

-19.40

-18.27

-18.95

ΔGb

-16.54

-23.47

-11.82

-12.43

RBA in silico

100

142

71

75

RBA in vitro

100

236

0.05

0.25

[a] ΔH = ΔEvdw + ΔEele + ΔGGB + ΔGSA [b] ΔG = ΔH – TΔS


Page 51 of 51 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2. Residues with binding affinities greater than -1 kcal/mol (kcal/mol) ERα-E2

ERα-DES

ERα-BPA

ERα-DPN

Residue Energy

Residue Energy

Residue Energy

Residue Energy

Met343

-1.16

Leu346

-1.57

Leu346

-2.03

Leu525

-1.98

Leu346

-1.55

Leu349

-1.25

Thr347

-1.12

Leu387

-1.58

Ala350

-1.00

Ala350

-1.22

Leu349

-1.05

Leu346

-1.47

Leu384

-1.26

Glu353

-1.66

Ala350

-1.03

Ala350

-1.42

Leu387

-1.81

Leu384

-1.18

Leu387

-1.13

Met421

-1.42

Met388

-1.06

Leu387

-2.09

Leu391

-1.40

Leu391

-1.30

Leu391

-1.03

Leu391

-1.04

Phe404

-1.10

Phe404

-1.25

Phe404

-1.80

Phe404

-1.58

Met421

-1.19

Thr347

-1.05

Met421

-1.06

Ile424

-1.24

Ile424

-1.30

His524

-2.82

Gly521

-2.23

Leu525

-1.21

Leu525

-1.37


The Effect of Structural Diversity on Ligand Specificity and Resulting

Recommend Documents