Xenoestrogen and Its Interaction with Human Genes and Cellular Proteins: An In-Silico Study

Background: Breast cancer represents one of the leading causes of death worldwide. Apart from genetic factors, the sex hormone estrogen plays a pivotal role in breast cancer development. We are exposed to a plethora of estrogen mimics on a daily basis via various routes. Nevertheless, how xenoestrogens, the exogenous estrogen mimics, modulate cancer-associated signaling pathways and interact with specific genes is still underexplored. Hence, this study aims to explore the direct or indirect binding partners of xenoestrogens and their expression upon exposure to these estrogenic compounds. Methods: The collection of genes linked to the xenoestrogens Octylphenol, Nonylphenol, Bisphenol-A, and 2,2-bis(4-hydroxyphenyl)-1,1,1-trichloroethane were gathered from the Comparative Toxicogenomics Database. Venny 2.1 was utilized to pinpoint the genes shared by these xenoestrogens. Subsequently, the shared genes underwent Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis using the Database for Annotation, Visualization, and Integrated Discovery bioinformatics resource. A xenoestrogen-protein interaction network was constructed using Search Tool for Interactions of Chemicals. The expressions of common genes were studied with the microarray dataset GSE5200 from the Gene Expression Omnibus database. Also, the expression of a common gene set within different breast cancer subtypes was identified using the University of California, Santa Cruz Xena. Results: The genes linked to xenoestrogens were identified, and 13 genes were found to interact with all four xenoestrogens. Through DAVID analysis, the genes chosen are found to be enriched for various functions and pathways, including pathways in cancer, chemical carcinogenesis-receptor activation, and estrogen signaling pathways. The results of the Comparative Toxicogenomics Database and the chemical-protein interaction network derived from STITCH were similar. Microarray data analysis showed significantly high expression of all 13 genes in another study, with Bisphenol-A and Nonylphenol treated MCF-7 cells, most of the genes are expressed in luminal A or basal breast cancer subtype. Conclusion: In summary, the genes associated with the four xenoestrogens were mostly linked to pathways related to tumorigenesis, and the expression of these genes was found to be higher in breast cancer.


Introduction
There is a growing concern about the menace caused by xenoestrogen exposure to various diseases.The estrogen-mimicking compounds, which are also known as endocrine-disrupting chemicals, play multiple roles in hormone biosynthesis, metabolism, and action that result in the alteration of normal hormonal homeostasis [1].Since these compounds are lipophilic in nature, which is an insecticide and its application on crops results in accumulation in the human body [8].Several studies show that natural estrogen plays a significant role in the development and spread of hormone-sensitive cancers.Since xenoestrogen possesses estrogen-like characteristics, it also accelerates the growth of estrogendependent malignancies like breast, ovarian, and uterine cancers [9].
Most of the studies involving xenoestrogens are focused mainly on the effect of a single endocrine disruptor on specific cancers [10].However, it has not been explored whether multiple xenoestrogens can interact with the same target genes or influence similar biological and molecular processes in carcinogenesis.Thus, in the present study, we are trying to identify the genes commonly interacting either directly/ indirectly with OP, NP, BPA, and HPTE through various bioinformatic platforms and also to unravel their correlation with tumorigenesis by understanding the pathways and molecular functions influenced by them.

Screening of common genes directly/indirectly interacting with the xenoestrogens
Venn diagram was constructed by entering the gene list obtained from CTD into Venny 2.1(https://bioinfogp.cnb.csic.es/tools/venny/index.html)online tool to obtain the common genes associated with OP, NP, BPA, and HPTE.

Gene Ontology (GO) Enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis
The genes interacting with at least three of the four xenoestrogens studied were selected to perform functional enrichment analysis using The Database for Annotation, Visualization, and Integrated Discovery (DAVID) (version 2021; david.ncifcrf.gov)[12,13] database.GO function enrichment and KEGG pathway analysis were carried out in David 2021 database by selecting the "Homo sapiens" species.The GO results of significant terms for cellular component (CC), biological process (BP), and molecular function (MF) were ranked by p-value (Benjamini method) and exhibited as bar charts.p-value < 0.05 was considered as statistically significant.

Chemical-Protein interaction
The direct/indirect interactions among xenoestrogens were analyzed, and proteins were elucidated using the Search Tool for Interactions of Chemicals (STITCH) 5.0 database (http://stitch.embl.de/).The species was selected as 'Homo sapiens.The optional setting for network construction was shown as follows: number of interactors = not more than 50; minimum required interaction score = 0.700

Gene expression analysis using the Gene Expression Omnibus (GEO) database
The genes commonly interacting with the four xenoestrogens were analyzed for their expression status in a microarray dataset (GSE5200) obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/).GEO is an open-access functional genomics data repository of gene expression datasets and sequence-based data.GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is a network analysis tool that allows users to analyze deregulated genes under different conditions and is used for data processing.p-value < 0.05 (Benjamini method) was considered as statistically significant.

Gene expression analysis in different cancers
The mRNA expression levels of commonly interacting 13 genes in various human cancers were obtained from cBioPortal (www.cbioportal.org).The expression of genes in various breast cancer subtypes was studied using data from 1218 breast cancer samples (BRCA-RNAseq) taken from the TCGA PanCancer database using the University of California, Santa Cruz (UCSC) Xena.The UCSC Xena (http://xena.ucsc.edu/) is a functional genomics browser that offers visualization and integration for examining and viewing public data repositories.p-value < 0.05 (One-way ANOVA) was considered statistically significant.

Genes induced by xenoestrogens
The CTD was used to retrieve all chemical-gene interactions induced by the xenoestrogens OP, NP, BPA, and HPTE.A total of 86 genes were identified to have direct/indirect interactions with OP, 625 genes with NP, 25,521 genes with BPA, and 229 genes with HPTE.This gene set was analyzed using Venny 2.1 to find the common interacting genes that are induced by the xenoestrogens.Thirteen genes were found to be at the intersection that interacts with all the four xenoestrogens in the Venn diagram (Figure 1).
A set of 68 genes was found to interact with at least three of the four xenoestrogens studied (Table 1).This gene set was used for GO function enrichment analysis and KEGG analysis using DAVID 2021 database.Biological processes, cell composition, and molecular function are the three categories addressed by GO   Amidst all the cancers listed in Figure 2, breast cancer showed the highest association with the gene set.

Identification of xenoestrogen-protein interaction
The biological relevance of xenoestrogens was estimated by evaluating the chemical-protein interactions and molecular networks by STITCH database.STITCH is a database of known and predicted interactions between chemicals and proteins.The interactions are based on computational prediction, knowledge transfer across species, and interactions gathered from other (primary) databases; they comprise direct (physical) and indirect (functional) correlations.Xenoestrogen-protein interaction was analyzed, and the selected protein targets with a probabilistic confidence score of 0.700 were plotted as an interaction network (Figure 3).The nodes and the edges represent protein/gene targets and their interactions, respectively.OP produced a network consisting of 2 nodes, 1 edge, an average node degree of 1, a clustering coefficient of 1, an expected edge number of 2, and a protein-protein interaction (PPI) enrichment p-value of 0.67.BPA produced a network consisting of 50 nodes, 159 edges, an average node degree of 6.36, a clustering coefficient of 0.751, an expected edge number of 90, and a PPI enrichment p-value of 3.53e-11.NP produced a network consisting of 19 nodes, 42 edges, an average node degree of 4.42, a clustering coefficient of 0.544, an expected edge number of 23, and a PPI enrichment p-value of 0.000334.HPTE produced a network consisting of 9 nodes, 13 edges, an average node degree of 2.89, a clustering coefficient of 0.87, an expected edge number of 10, and a PPI Enrichment p-value of 0.234 (Figure 3).The gene set obtained from CTD was comparable to the outcomes obtained from STITCH database.

Gene expression analysis using GEO
The microarray dataset GSE5200 of MCF7 cells exposed to 1 nM BPA and 10 nM NP for 48 hours was retrieved from GEO.The thirteen genes commonly interacting with all four xenoestrogens obtained from CTD were analyzed for their expression pattern in this dataset using the GEO2R tool.The expression of ESR1, ESR2, MAPK3, TGFB3, BCL2, DDIT3, HSD3B1, CYP19A1, CLU, FOS, and PGR genes was found to be significantly upregulated in both the BPA and NP treated cells when compared to the control (vehicle -0.1% ethanol) (Figure 4).

Differential expression of the common gene set in breast cancer subtypes
The genes commonly interacting with all the xenoestrogens studied were analyzed for their expression in multiple malignancies as well as in breast cancer subtypes using cBioportal and UCSC Xena, respectively.RNA-seq data obtained from cBioportal indicated higher expression status for all the genes analyzed in multiple malignancies.Of note, breast cancer had the highest ESR1 expression among all solid malignancies.In addition, most of our genes of interest showed moderate to high levels of expression in breast cancer (Figure 5).Given that the majority of the genes exhibit increased expression in both breast cancer and in the MCF-7 cell line, we further examined the expression of the 68 genes, which were commonly expressed by at least three xenoestrogens on different subtypes of breast cancer using the UCSC Xena platform.Out of the 68 genes, 9 genes were excluded as it showed a median value of 0. Among the remaining 59 genes, 40.67% of the genes were overexpressed in the Luminal A subtype, 33.8% in the basal type, 22.03% in the Her2, and 11.86% in the Luminal B types (Table 2).

Discussion
Breast cancer represents one of the leading causes of death worldwide.More than 80% of breast cancer cells grow in response to the hormone estrogen [14].Apart from genetic factors, the sex hormone estrogen plays a pivotal role in breast cancer development and progression.Increased lifetime estrogen exposure brought on by early menarche, late menopause, long-term menopausal estrogen therapy, and exposure to estrogen mimics are linked to an increased incidence of breast cancer.The most prevalent explanation behind estrogen's role in tumorigenesis is estrogen, upon binding with estrogen receptor α (Erα), boosts cell proliferation and causes mutations, which occur as a result of errors occurring during DNA replication.The proliferation of cells carrying    In the present study, the genes directly or indirectly associated with the xenoestrogens OP, NP, BPA, and HPTE were retrieved from CTD. ESR1, ESR2, MAPK1, MAPK3, PGR, TGFB3, TNF, DDIT3, CYP19A1, CLU, FOS, HSD3B1, and BCL2 were the common set of genes obtained, which interact with all the four xenoestrogens analyzed.ESR1 and ESR2 encode estrogen receptors, which are important for sexual development and reproductive function [17].As these cancer cells were treated with xenoestrogens, the expression of estrogen receptor were expected.ERα, which is encoded by the ESR1 gene, is well correlated with breast cancer, whereas, ERβ encoded by ESR2, is mostly expressed by normal breast epithelial cells [18].MAPK1 and MAPK3 are two genes coming under the mitogen-activated protein kinase family [19].They strictly regulate cellular processes like cell proliferation, differentiation, and survival.Any dysfunction in the MAPK pathway triggers tumorigenesis [20].PGR is a member of the steroid receptor family, which codes for progesterone receptor mutations and polymorphisms, which has been reported to cause tumorigenesis and risk for cancers, including ovarian, breast, and endometrial cancers [21].Transforming growth factor β (TGFβ) has three isoforms; interestingly, previous reports showed Figure 5. Graphical Representation of the Expression Levels of 13 Genes Commonly Interacting with Xenoestrogens in Human Cancers.Data obtained from cBioPortal For Cancer Genomics.TCGA The Cancer Genome Atlas Network higher TGFβ1 and TGFβ3 levels in breast cancer tissues when compared to normal tissues.Moreover, they also found increased TGFb1 and TGFb3 and decreased TGFb2 expression in advanced lymph node (LN)-positive and metastatic tumors [22].TNF (Tumor necrosis factor) is a double-edged sword that can either act as a pro-or anti-tumorigenic factor.The anti-tumorigenic property of TNF causes the death of cancer cells; however, in cancer cells that are resistant to TNF-induced death, it promotes angiogenesis, migration, and cell proliferation [23].A study by Lin et al., 2020 demonstrated that gastric cancer showed elevated levels of DDIT3 (DNA damage-inducible transcript 3) and activation of DDIT3-mediated pathways [24].High levels of CYP19A1 mRNA expression are reported to be present in women with malignant tumors of the breast, endometrial, and ovary.Additionally, increased CYP19A1 mRNA levels were strongly linked to the likelihood of metastases, local recurrence, and breast cancer-related deaths [25].The cell activities involved in the development of tumors and carcinogenesis have all been linked to clusterin (CLU).Expression of CLU has been linked to the development of a number of malignancies, including cancers of the prostate, colon, and breast [26].Muhammad et al. observed upregulation of c-FOS in Head and neck squamous cell carcinoma (HNSCC) sphere-forming cells than in parental cells.In their study, the tumorigenic phenotype in immunodeficient mice is displayed by exogenous expression of c-Fos in non-tumorigenic HNSCC MDA1386Tu cells.Overall, in non-tumorigenic cell lines, they have demonstrated that overexpression of c-Fos makes the cells tumorigenic and increases the expression of EMT/CSC marker genes [27].A study by, Hearn et al., 2020 proposed the germline HSD3B1 genotype as a genetic biomarker of resistance to Androgen Deprivation Therapy (ADT) for prostate cancer [28].Human tumor tissues typically express high levels of BCL-2 family members that prevent apoptosis, such as BCL-2 or BCL-XL, which prevents tumor cells from dying and causes them to spread quickly [29].All the thirteen genes interacting with the xenoestrogens play a pivotal role in the development and progression of multiple malignancies.
GO enrichment and KEGG pathway analysis were carried out using the DAVID online tool to identify BP, CC, and MF and pathways involving the 68 commonly expressed genes.With regard to BP, these genes are mainly enriched in signal transduction and positively regulating transcription, gene expression, and cell proliferation.Genes enriched in the CC category are the nucleus, cytosol, extracellular space, and nucleoplasm.KEGG pathway analysis indicated the genes are mainly associated with pathways in cancer, chemical carcinogenesis-receptor activation, MAPK pathway, etc.These genes are found to be mostly enriched in breast cancer among all the other cancers.This could be due to the fact that estrogenic compounds have adverse effects on hormone-sensitive organs like the breast [30].Further, we looked into the gene expression pattern shown by MCF-7 cells treated with BPA and NP using a microarray dataset obtained from GEO.The in-silico analysis demonstrated higher expression of a set of genes in xenoestrogenexposed cells when compared to untreated counterparts, except for TNF and HSD3B1 expression upon treatment with BPA.A probable cause for this could be the lower concentration of BPA, which would not be enough to trigger the pro-inflammatory role of TNF.Besides, the reduced expression of HSD3B1, which is an enzyme that converts pregnenolone to progesterone, was observed probably due to a decrease in progesterone secretion at this concentration of BPA [31].
Since most of the genes showed increased expression in various cancer types, especially in breast cancer, the gene expression patterns that commonly interact with at least three xenoestrogens were examined in breast cancer subtypes using UCSC XENA.The gene expression analysis shown by the BRCA-RNAseq data taken from the TCGA PanCancer database disclosed a notable gene expression pattern in four breast cancer subtypes.A higher percentage of genes falls into the luminal A subtype of breast cancer, followed by basal, her2, and lastly, luminal B. Since many of these xenoestrogens have been shown to function through the estrogen receptor at the molecular level, this could account for the higher percentage of luminal A subtype of breast cancer [32].
Furthermore, our data point to the need for greater attention to be paid to the functions of xenoestrogens in the onset and progression of Triple Negative Breast Cancer (TNBC).Besides Erα, it has been demonstrated that xenoestrogen functions via ERRγ, G protein-coupled estrogen receptor (GPER), ER-β, and ER-α36 (a variant of ER-α) in TNBC [33].This implies that xenoestrogen exposure could be an important cause of high and increasing rates of hormone receptor-positive as well as hormone receptor-negative breast cancer.
In conclusion, the present study used bioinformatics tools to identify the genes and proteins interacting with xenoestrogens.We selected a set of thirteen genes that were found to commonly interact with OP, NP, BPA, and HPTE.These genes were enriched in cancer pathways, chemical carcinogenesis-receptor activation, MAPK signaling pathway, and estrogen signaling pathway.Also, higher expression of genes that interact with xenoestrogens was observed in breast cancer cell lines as well as in multiple cancer types; the expression was found to be high in the breast cancer Luminal A subtype.However, further studies are required to prove this experimentally.In the future, a better understanding of the genes that interact with xenoestrogens and their association with different cancer types enable us to identify the causative agent, improve cancer treatment, and also help us in cancer prevention to an extent.

Figure 2 .
Figure 2. GO Function Enrichment and KEGG Pathway Analysis of Genes Interacting with at Least Three Xenoestrogens as Obtained from the DAVID 2021 Database.A. BP: biological process; B. CC: cellular component; C. MF: molecular function; D. KEGG: Kyoto Encyclopedia of Genes and Genomes; E. Cancers

Figure 3 .
Figure 3. Schematic Illustration of the Chemical-Protein Networks of A. OP, B. BPA, C. NP and D. HPTE and Its Interacting Entities, Acquired from STITCH Database.It was generated according to the known and predicted interactions for Homo sapiens with the confidence score set to 0.700 with the maximum of 50 interactions.Thicker lines represent the stronger linkages.Gray and green lines show the protein-protein interaction.Individual nodes represent gene products.Small nodes represent proteins of unknown 3D structure, and large nodes represent proteins of known or predicted 3D structure.

Figure 4 .
Figure 4. mRNA Expression of the Thirteen Genes Commonly Interacting with Xenoestrogens.Data was obtained from the gene expression dataset (GSE5200) in NCBI, GEO database and analysed using GEO2R tool.

Figure 5 .
Figure 5. Graphical Representation of the Expression Levels of 13 Genes Commonly Interacting with Xenoestrogens in Human Cancers.Data obtained from cBioPortal For Cancer Genomics.TCGA The Cancer Genome Atlas Network

Table 1 .
The most enriched GO terms in the BP category were signal transduction, positive regulation of transcription, gene expression, cell proliferation, and apoptosis.MF was mainly enriched in protein binding, DNA binding, List of Intersecting Genes Obtained from Venny 2.
1 Figure 1.Venn Diagram Representing the Number of Genes that are Uniquely Expressed by each Xenoestro Gens OP, NP, BPA and HPTE.The overlapping regions shows the number of genes that are expressed in two or more samples.The diagram was built using Venny 2.1 online tool (http://bioinfogp.cnb.csic.es/tools/venny/)according to the data from the CTD database.

Table 2 .
Differential Expression of Genes Commonly Interacting with Xenoestrogens in Various Breast Cancer Subtypes.Based on the data of 1218 breast cancer samples (BRCA-RNAseq) taken from TCGA PanCancer database using UCSC Xena.Log2 normalized values are represented here.

Table 2 .
Continued macromolecular complex binding & RNA polymerase II transcription factor activity, and sequence-specific DNA binding.The most enriched GO terms in the CC category were nucleus, cytosol, nucleoplasm, and extracellular space.KEGG pathway analysis revealed that the common gene set was associated with pathways related to cancer, chemical carcinogenesis-receptor activation, MAPK signaling pathway, and estrogen signaling pathway.