Ligand Binding Domain of Estrogen Receptor Alpha Preserve a Conserved Structural Architecture Similar to Bacterial Taxis Receptors

It remains a mystery why estrogen hormone receptors (ERs), which are highly specific toward its endogenous hormones, are responsive to chemically distinct exogenous agents. Does it indicate that ERs are environmentally regulated? Here, we speculate that ERs would have some common structural features with prokaryotic taxis receptor responsive toward environmental signals. This study addresses the low specificity and high responsiveness of ERs toward chemically distinct exogenous substances, from an evolutionary point of view. Here, we compared the ligand binding domain (LBD) of ER alpha (α) with the LBDs of prokaryotic taxis receptors to check if LBDs share any structural similarity. Interestingly, a high degree of similarity in the domain structural fold architecture of ERα and bacterial taxis receptors was observed. The pharmacophore modeling focused on ligand molecules of both receptors suggest that these ligands share common pharmacophore features. The molecular docking studies suggest that the natural ligands of bacterial chemotaxis receptors exhibit strong interaction with human ER as well. Although phylogenetic analysis proved that these proteins are unrelated, they would have evolved independently, suggesting a possibility of convergent molecular evolution. Nevertheless, a remarkable sequence divergence was seen between these proteins even when they shared common domain structural folds and common ligand-based pharmacophore features, suggesting that the protein architecture remains conserved within the structure for a specific function irrespective of sequence identity.

It remains a mystery why estrogen hormone receptors (ERs), which are highly specific toward its endogenous hormones, are responsive to chemically distinct exogenous agents. Does it indicate that ERs are environmentally regulated? Here, we speculate that ERs would have some common structural features with prokaryotic taxis receptor responsive toward environmental signals. This study addresses the low specificity and high responsiveness of ERs toward chemically distinct exogenous substances, from an evolutionary point of view. Here, we compared the ligand binding domain (LBD) of ER alpha (α) with the LBDs of prokaryotic taxis receptors to check if LBDs share any structural similarity. Interestingly, a high degree of similarity in the domain structural fold architecture of ERα and bacterial taxis receptors was observed. The pharmacophore modeling focused on ligand molecules of both receptors suggest that these ligands share common pharmacophore features. The molecular docking studies suggest that the natural ligands of bacterial chemotaxis receptors exhibit strong interaction with human ER as well. Although phylogenetic analysis proved that these proteins are unrelated, they would have evolved independently, suggesting a possibility of convergent molecular evolution. Nevertheless, a remarkable sequence divergence was seen between these proteins even when they shared common domain structural folds and common ligand-based pharmacophore features, suggesting that the protein architecture remains conserved within the structure for a specific function irrespective of sequence identity.

INTRODUCTION
At present, the evolution of nuclear hormone receptors (NRs) is traced back up to protostomes such as plathelminths and mollusks (Kohler et al., 2007;Vogeler et al., 2014;Baker and Lathe, 2018;Wu and LoVerde, 2019). According to the present perspective, Amphioxus, a prochodate deuterostome, contains the most primitive estrogen receptor (ER) (Baker and Chandsawangbhuwana, 2008;Lecroisey et al., 2012). However, the full transcriptional function of the ER evolved hundreds of millions of years later after the evolution of early deuterostomes is reported (Kao et al., 2000;Baker and Chandsawangbhuwana, 2008;Callard et al., 2011). ERα is believed to be the most primitive among all the NRs reported so far. Sequence analysis indicates that the other NRs such as other ER subtypes, progesterone receptors, androgen receptors, glucocorticoid receptors, and minerocorticoid receptors evolved from ERα through gene duplication and sequence divergence (Elliston and Katzenellenbogen, 1988;Thornton, 2001;Guerriero, 2009;Baker et al., 2015;Baker, 2019).

LBD OF ESTROGEN HORMONE RECEPTORS ARE SPECIFIC TOWARD ITS ENDOGENOUS LIGAND
All the known NRs discovered so far adopts common structural and functional characteristics and consists of evolutionarily conserved, structurally and functionally distinct three to six basic domains (Lakshmanan and Shaheer, 2020). ERs are made up of an N-terminal domain, central DNA binding domain, C-terminal LBD, and two distinct, conformationally active regions designated as activation function 1 (AF-1) and activation function 2 (AF-2) (Lakshmanan and Sadasivan, 2014). DNA-binding domain (DBD) is the most conserved among the domains, and the N-terminal domain is most variable in sequence and length. ER signaling depends on the ligand/hormone and begins with the binding of ligand to LBD (Lakshmanan and Shaheer, 2020). LBDs of nuclear hormones are specific toward its endogenous ligand and normally do not cross-interact with other non-specific endogenous hormones (Sasson and Notides, 1983;Klinge, 2001;Razandi et al., 2004;Yang et al., 2004;Morissette et al., 2008). On the other hand, they are responsive to metabolites and precursors of steroid hormones and wide variety of chemically distinct exogenous substances grouped as xenohormones. Both ERα and ERβ exhibits similar affinities for the endogenous hormones estradiol, estriol, and estrone (Sasson and Notides, 1983;Morissette et al., 2008).

WHY ERs ARE RESPONSIVE TO AN ARRAY OF DIVERSE EXOGENOUS SUBSTANCES?
Estrogen hormone receptors exhibits least substrate specificity and strong binding affinities for most of the steroid/phenolic and few of the non-phenolic/non-steroid agents classified as xenoestrogens (Cauley, 2015;Gore et al., 2015;Shafei et al., 2018;Pacyga et al., 2019;Rosenfeld and Cooke, 2019;Basak et al., 2020;Lakshmanan and Shaheer, 2020;Park et al., 2020). It is interesting to explore why a receptor that shows high specificity to its endogenous ligand can be transactivated by other hormone metabolites or chemically distinct exogenous agents. The present study was initiated to address this question through bioinformatics tools.

LBD Structural Architecture Is Conserved Between ERs and Bacterial Taxis Receptors Irrespective of Sequence Divergence
Emergence of proteins with novel domain and/or domain combinations, generated by either homologous or nonhomologous DNA repair or inserted into the genomes by transposition, is an important evolutionary mechanism, as it confers potentially diverse biological functions for the organisms (Itoh et al., 2007;Peisajovich et al., 2010;Forslund et al., 2019). Such domain accretion would have paved way in the emergence of novel proteins with specialized functions in the due course of eukaryotic evolution. While ERs are highly specific toward their endogenous hormone, they exhibit least substrate specificity and strong binding affinities toward array of exogenous ligand. This prompted us to speculate that LBD of ER alpha, the NR, would have some similarity with prokaryotic receptors responsive toward environmental signals (chemotaxis receptors). To prove this speculation, we downloaded the amino acid sequences and structures of LBDs of ERα and compared it with bacterial taxis receptors and sequences from other lower life forms such as protozoa and protostoma ( Table 1) using Constraint-Based Multiple Alignment Tool (COBALT) and sequence-independent structural alignment (TM-align).
Analysis based on sequence similarity such as Basic Local Alignment Search Tool (BLAST), BLAST-Like Alignment Tool (BLAT), and CLUSTAL multiple sequence alignment tools are not reliable in detecting the homology of distantly related species if similarity of the protein sequence is <30% (Bhagwat et al., 2012;Moreno-Hagelsieb and Hudy-Yuffa, 2014;Ward and Moreno-Hagelsieb, 2014). Nevertheless, protein domain composition is anticipated to be conserved throughout the course of evolution due to functional constraints (Forslund et al., 2019;Yu et al., 2019). The function of a particular protein, to a great degree, is decided by the orderly arrangement of protein sequence to specific domains that constitutes domain architecture. Proteins with the same domain architecture may probably have similar structures and hence related cellular function (Koonin et al., 2002).
Constraint-based multiple alignment tool has advantage over other protein multiple sequence alignment tools in that it finds a set of pairwise constraints derived from conserved domain database and protein motif database and sequence similarity, using RPS-BLAST, BLASTP, and PHI-BLAST. Pairwise constraints were then integrated into a progressive multiple alignment (Papadopoulos and Agarwala, 2007). All the protein sequences were downloaded from protein sequence database from NCBI. Protein structural files were downloaded from RCSB, protein databank, and COBALT and were used for multiple sequence alignment and prediction of common domain architecture. Molecular Evolutionary Genetics Analysis (MEGA) version X was used for constructing phylogenetic trees based on molecular evolution.
Interestingly, using COBALT, we found that LBDs of ER showed high conservation in domain architecture with all the selected sequences ( Figure 1A and Supplementary Figure 1). Furthermore, we selected only the first four proteins shown in Table 1 and performed the COBALT alignment. All the four sequences, that is, LBDs of ER, nuclear receptor of Ciona sp., and two bacterial taxis receptors, exhibited high degree of conservation in domain architecture ( Figure 1B and Supplementary Figure 2: first four sequence). As a validation to our hypothesis, we randomly selected some protein sequences (Supplementary Table 1) and repeated COBALT multiple sequence alignment. The result (Supplementary Figure 3) showed that the selected sequences have less conservation and more gaps.
The phylogenetic tree built for all the six sequences with mega gave an unrooted gene trees ( Figure 1C) with clustering of genes into two separate groups. Human ER was clustered into a common clade with Ciona sp., and Mytilus sp. was an outgroup from this clade. Two bacterial chemoreceptors were clustered together. Tom40 of Amoeba was an outgroup from the clade of bacterial receptors. Phylogenetic analysis showed that ERs and bacterial chemotactic receptors are evolutionarily unrelated.
Irrespective of their sequence and structural and functional diversity, ERs and other bacterial chemoreceptors such as taxis to serine and repellents (Tsr) and broad ligand-specific (BLS) receptor of Pseudomonas putida used in this study have certain common features such as a ligand binding domain (LBD) and homodimerization of the LBDs to exert their downstream signaling.
TM-align has advantage over other pairwise sequence alignment tools, as it confers sequence-independent protein structure comparison (Zhang and Skolnick, 2005). For any given two protein structures of unidentified similarity, the program uses heuristic dynamic programming iterations and initially generates optimized amino acid-to-amino acid alignment based on structural similarity. The tool then returns an optimal superposition of the two structures with a TM-score that can be used to scale the similarity. We performed TM-align to find a sequence-independent alignment based on local backbone similarity using heuristic dynamic programming iterations. The structures of LBDs of the serine receptor (Tsr) (PDB ID: 3ATP) of Escherichia coli with the LBD of ER-alpha (PDB ID: 1ERE) and LBDs of the chemoreceptor (BLS) (PDB ID: 6S33) of P. putida with the LBD of ER-alpha (PDB ID: 1ERE) were selected for TM-align.
Structural alignments with TM-align confirmed similar structural folds between LBDs of ER and LBDs of bacterial  Frontiers in Ecology and Evolution | www.frontiersin.org chemoreceptors (Figures 1D,E). The superposition of two proteins using T-align showed a clear picture on the similarity of Domain Structural architecture between Tsr and ER and between BLS and ER, substantiating our hypothesis on preservation of common fold architecture between the LBDs between bacterial chemoreceptor and human ERα. Frontiers in Ecology and Evolution | www.frontiersin.org

Ligands Specific for the Estrogen Receptor and Bacterial Chemotaxis Receptors Shares Common Pharmacophore Features
Pharmacophore models provide description of optimal supramolecular interactions that typically occur between small-molecule ligands and their respective protein receptor. The pharmacophores of ERs ligands and bacterial chemotaxis receptor ligands were compared to understand if these ligands shared any common pharmacophore features. Two pharmacophore hypotheses were developed in which the first hypothesis involves 17β-estradiol (human estrogen) as the reference ligand and the second hypothesis was focused on coumestrol (Phytoestrogen) as the reference ligand. The hypotheses were developed using the "Develop pharmacophore hypothesis" module of Maestro 11.0 Schrodinger suite. To develop the hypotheses using the reference ligands, hydrogen bond (HB) donors (D), HB acceptors (A), aromatic ring (R), hydrophobic/non-polar groups (H), and positive/negative ionizable groups were selected as the pharmacophore features. Furthermore, the ligands to be screened for matching pharmacophore features were subjected to "ligand-database pharmacophore screening" using the same software. Interestingly, pharmacophore modeling uncovered a set of common features among the ligands specific for the ER and bacterial chemotaxis receptors (Figures 1F-H). All the bacterial chemotaxis receptor ligands possessed pharmacophore features that are common for 17β-estradiol and coumestrol (Figures 1F-H). The bacterial chemotaxis receptor ligands such as salicylic acid, benzoic acid, and serine exhibited matching pharmacophores with the phytoestrogen "coumestrol, " a known ER agonist.

Natural Ligands of Bacterial Chemotaxis Receptors Exhibited Favorable Binding Affinity Toward Human Estrogen Receptor and Vice Versa
Our finding that the structural architecture was conserved between the LBD of ERs and bacterial Taxis receptors along with the existence of common pharmacophore feature between the ligands of both receptors prompted us to check for the interactions of the ligands with both receptors. The molecular docking studies were performed using Schrodinger suite 11.0. The three dimensional structures of ER (PDB ID 1ERE), E. coli chemotaxis receptor Tsr (PDB ID 3ATP), and P. putida PcaY_PP LBD (PDB ID 6S33) were downloaded from Protein Data Bank. The downloaded protein structures were processed and energy minimized using the Protein Preparation Wizard. The protein preparation involved assigning bond orders, adding hydrogens, creating disulfide bonds, and creating zero-order bonds to metals. The het states of the proteins were generated using "Epik" at pH 7 ± 2.0, and restrained minimization was performed using OPLS3 forcefield. Furthermore, a receptor grid was generated from each structure by selecting the cocrystal ligand as the centroid of the receptor. The ligand molecules (Figure 2) were downloaded from NCBI-PubChem database and prepared for docking with the LigPrep module. The ligand molecules were structurally optimized at near neutral pH (7 ± 1) and subjected to energy minimization using OPLS3 force field.
Furthermore, the prepared ligand sets were docked against the receptor grid of the above-mentioned target proteins. Extra precision (XP) flexible docking was performed using the GLIDE module, and the affinity of the ligands toward the target proteins were ascertained in terms of negative glide score (kcal/mol). Fifty docking poses were collected for each ligand. The binding poses were generated using Pymol software (free license). The glide score is an indirect measure of binding free energy of the ligand-receptor interaction. The more negative the glide score, the stronger the interaction.
For the molecular docking studies against human ER, 17β-estradiol was used as the reference ligand, and other ER ligands were used for validation of the experiment (Figure 2A). Similarly, quinic acid was used as the standard reference ligand for bacterial chemotaxis receptors, and other ligands specific for bacterial receptor was used for cross-validation of the experiment. The docking studies strongly support the findings of structural similarity between human ER and bacterial chemotaxis receptors that were presented in the structure alignment section. When discussing the target specificity of ER ligands against bacterial chemotaxis receptors, all the ligands exhibited thermodynamically favorable interaction ( Figure 2B). Phytoestrogens such as coumestrol, diadzeine, and naringenin exhibited stronger interaction with high binding affinity toward bacterial chemotaxis receptors than the estradiol or synthetic ER ligands (Figures 2C,D). On the other hand, the natural ligands of bacterial chemotaxis receptors exhibited strong interaction with human ER (Figures 2A,B). The favorable interactions of all ER ligands toward bacterial receptors irrespective of their binding strength and vice versa were the key findings of this experiment.

CONCLUSION AND PERSPECTIVE
We observed a remarkable similarity in the domain structural fold architecture of ERα and bacterial taxis receptors. The ligand molecules of both receptors also shared common pharmacophore features. The natural ligands of bacterial chemotaxis receptors exhibit favorable binding affinity with human ER. This holds true for ER-ligand interaction with bacterial receptors as well. However, phylogenetic analysis proves that these proteins are unrelated. Together, these results suggest that these receptors have evolved independently to respond against certain environmental agents, pointing toward a possibility of convergent molecular evolution. However, during the course of evolution, even when these sequences exhibited high divergence and acquired novel functions, the domain structural fold remained greatly preserved. Domain structural architecture for a particular acquired or specialized functions maybe conserved across species in due course of evolution even when the sequence composition varies. This hypothesis also explains the non-specific binding of the nuclear hormones toward an array of chemically distinct exogenous ligands and at the same time maintaining its high specificity toward the respective hormones.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material."

AUTHOR CONTRIBUTIONS
DL conceived the idea, developed the concept, designed the experiments, performed the experiments, and wrote the manuscript. SH contributed in doing revision experiments in the revised version of the manuscript. Both authors contributed to the article and approved the submitted version.

FUNDING
This research was funded with Seed grant (YU/Seed grant/065-2018), from Yenepoya (Deemed to be University).

ACKNOWLEDGMENTS
We would like to acknowledge Ranajith Das, Yenepoya Research Centre, for his comments and suggestions on the phylogenetic tree construction and analysis.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021. 681913/full#supplementary-material Supplementary Table 1 | Sequence ID and their Function for the random sequences selected for COBALT analysis.