Bisphenol S Promotes the Progression of Prostate Cancer by Regulating the Expression of COL1A1 and COL1A2.

In recent decades, Bisphenol S (BPS), which have been considered as alternatives for Bisphenol A (BPA), have become widely used in personal care products, paper products and food. Clarifying the relationship between bisphenol and tumors is of great signi�cance for the treatment and prevention of diseases. In this work, we discovered a new method to predict the correlation between bisphenol interactive genes and tumors. The transcriptome pro�le and interactive genes of bisphenol were obtained from the Cancer Genome Atlas and Genotype-Tissue Expression, Comparative Toxicology Genomics and PharmMapper database. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis showed that interactive genes were mainly enriched in prostate cancer. Gene targetd prediction and gene set variation analysis proved that bisphenol exert potential effects on prostate cancer. The operating characteristic curves and survival analysis showed that role of COL1A1 and COL1A2 in predicting the prognosis of prostate cancer. Cell counting kit-8 assay demonstrated that the cells with BPS-treated could remarkably promote the cell proliferation ability in both PC-3 and LNCap cells. Wound healing assay and the transwell assay demonstrated that the cells with BPS-treated could signi�cantly promote the cell invasion capacity of prostate cells. Two key genes, COL1A1 and COL1A2, were signi�cantly upregulated with BPS-treated in the PC-3 and LNCap cells.


Introduction
In the past few decades, with urbanization and industrialization levels continue to rise, environmental pollution has become a worldwide problem.Endocrine disrupting chemicals (EDCs), which can be found in agriculture and industry, have been proved to affect hormone synthesis, metabolism and function (Ma et al.2019).Bisphenol, as a widely used EDC, can be found in plastics, receipts, food packaging, and other products.regulate hormone signal paths and various other biological functions in a wide range of products (Murata et al. 2017).Relevant studies have shown that EDs can interact with some endocrine signaling pathways, such as estrogen receptor (Legeay et al. 2017).Estrogen receptor can regulate estrogen synthesis and related enzyme activities by transcription factors.Classical estrogen receptors include two subtypes of ERα and ERβ, both of which have similar structures.ERα is mainly expressed in hormone-targeted tissues, such as breast, uterus, testes and ovaries.On the contrary, ERβ is mainly expressed in the prostate, bladder, lung and intestine (Rochester et al. 2015).
In the past 40 years, the incidence and mortality of tumors have been increasing.As is estimated by the American Cancer Society, more than 1.7 million new cancers will be diagnosed in 2019, and there will be approximately 607,000 cancer-related deaths (Giri et al. 2020).Prostate cancer is one of the most common malignant tumors and the leading cause of cancer death in men.As is reported, prostate cancer comes in the second place after lung cancer in men (Labbé et al. 2018).In addition to age, race and family history of prostate cancer, prostate cancer susceptibility genes are also a key factor in the development of prostate cancer (Sathianathen et al. 2018).Relevant studies have con rmed that some DNA repair genes, including CHEK2, PALB2, BRIP1 and NBS1, may lead to a higher risk of prostate cancer.
The risk assessment of prostate cancer usually combines the patient's age, clinical tumor stage, serum PSA concentration, and Gleason score for a comprehensive analysis (Barata et al. 2019).Among them, EDCs play a key role in the occurrence and development of prostate cancer (Hu et al. 2012).Androgen receptor (AR) is very important for prostate development and tumorigenic growth.A large amount of evidence demonstrated that human exposure to EDCs can change the function of hormones, and nally promote the occurrence and development of prostate cancer (Litwin et al. 2017).As an endocrine disruptor that interact with hormone, BPA can directly interact with both ER-α and ER-β and may enhance prostate cancer growth (Shafei et al. 2018).Besides, patients with prostate cancer showed higher levels of BPA in their urine, which suggests that urine BPA levels may have prognostic value for prostate cancer (Hess-Wilson.2005).As an alternative for BPA, BPS has been increasingly detected in consumer products, foodstuffs, water, sediment and sludge in recent years (Fossati et al. 2016).Studies have reported the occurrence of BPS in human urine, among 315 urine samples collected from the United States, BPS was found in 81% of the total samples, and the average concentration of BPS was 0.17 ng/ml, which indicated that BPS may have toxic effects on the genitourinary system (Eladak et al. 2014).
Although the epidemiological literature on the toxicity of BPS are growing rapidly, however, its quantity and quality are still limited.Many studies are essentially cross-sectional studies, which makes results di cult to verify.The Comparative Toxicogenomics Database (CTD) is a literature-based and manually curated public resource and give us advance understanding regarding the effect of environmental exposure on human health.Besides, two tumor databases are also applied for further analysis.The Cancer Genome Atlas (TCGA), as a landmark cancer genomics program, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types.GEO is a public functional genomics data repository supporting MIAME-compliant data submissions and provide arraybased and sequence-based data.The development of prostate cancer involves the damages to the normal prostate transcription network.Here, we evaluate the key genes that associated with prostate cancer.Besides, the interaction of toxins with their interaction genes helps understand the multiple pharmacology of bisphenols and their effects on biological networks.TCGA and GEO database were applied to explored the possible mechanism underlying the association between prostate cancer and bisphenols-exposure.

Datasets downloaded
The data of interactive genes that associated with bisphenols were obtained from the CTD database (http://ctdbase.org/) in July 2021.Data on 8 bisphenols (Bisphenol A, Bisphenol AF, Bisphenol B, Bisphenol C, Bisphenol E, Bisphenol F, Bisphenol S, Bisphenol Z), as well as their gene, protein, and disease interactions, were downloaded.On the basis of 2D structure of bisphenols, they have similarities in structure.Two phenol groups are connected to a carbon atom at the 4'position.(Fig. 1)

GO and KEGG pathway enrichment analysis
In order to explore the potential biological processes (BPs), cellular component (CC) and molecular functions (MFs) of the interactive genes associated with bisphenols, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotation analyses were applied.The "clusterPro ler" and "ggplot2" packages in R were used to perform GO and KEGG functional annotation pathway enrichment analysis.P < 0.05 was selected as the threshold in the analysis.

Protein-protein interaction network
In order to further explore the relationship between interactive genes and protein, the protein-protein interaction (PPI) network of interaction genes associated with bisphenol was conducted by using STRING (https://www.string-db.org/).We next counted the nodes associated with each interaction gene.Then, Cytoscape 3.7.2 was applied to analyze the PPI network.

Target prediction
The structures of bisphenols (Bisphenol A, Bisphenol AF, Bisphenol B, Bisphenol C, Bisphenol E, Bisphenol F, Bisphenol S, Bisphenol Z) were obtained from PubChem (https://pubch em.ncbi.nlm.nih.gov/).Then, we the structure of the bisphenols was transferred into SDF format.PharmMapper database (http://lilab.ecust.edu.cn/pharmmapper/) is an updated integrated pharmacophore matching platform with statistial method for potential target identi cation.The Images of the 3D bisphenols metabolites structure are uploaded to the Pharmmapper database to determine the prediction targets related to bisphenols metabolites.Finally, in order to increase the credibility of the results, all the predicted targets in bisphenols metabolites were uploaded to the UniProt database (https://www.uniprot.org/)for screening.

Identi cation of prognosis based on interactive genes
In order to explore the key genes associated with prostate cancer and bisphenol.Kaplan-Meier survival curves was used to compare the prognosis difference between two groups.We downloaded the data from the TCGA-PC datasets.P < 0.05 was de ned as a signi cant difference.

Validation of the interactive genes involved in prostate cancer and bisphenol
Receiver operating characteristic (ROC) curve was used to assess the stability of the prognosis model.The area under roc curve (AUC) value > 0.6 was considered to have good predictive capabilities.

Gene set variation analysis of key genes
Gene set variation analysis of key genes was conducted to compare the difference of biological pathways between two groups.The reference gene set was "Hallmark.7.4.symbol.gmt".
Cell Counting Kit-8 Assay CCK8 assay were performed using a CCK8 kit (Dojindo, Shanghai, China).Brie y, cells were plated into 96-well plates with 2000 cell per well and added 10µL CCK8.Cells were maintained for 1.5h in a humidi ed incubator at 37°C with 5% CO2.The absorbance at 450nm was used to express the cell proliferation ability.

RNA isolation and quantitative real-time PCR
The total RNA was isolated from cell lines using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer's protocol.The cDNA was synthesized with PrimeScript™ RT reagent Kit (Takara).The purity and concentration of RNAs were assessed by Nanodrop ND-2000 spectrophotometer (Thermo Scienti c, EUA).Quantitative real-time PCR (qRT-PCR) was conductedusing AceQ® qPCR SYBR Green Master Mix (vazyme) and according to the manufacturer's instructions with the Roche LightCycler 480 (Applied Biosystems, CA, USA).GAPDH was used as internal control and each sample was repeated four times.The fold change in the expression of genes was calculated with the formula 2-ΔΔCT.The primer sequence used was as follows: COL1A1, forward primer, 5'-GAGGGCCAAGACGAAGACATC-3', reverse primer, 5'-CAGATCACGTCATCGCACAAC-3'; COL1A2, forward primer, 5'-GTTGCTGCTTGCAGTAACCTT-3', reverse primer, 5'-AGGGCCAAGTCCAACTCCTT-3'.

Cell migration assays
Cell migration was determined by the wound healing assay and the transwell assay.(1) The woundhealing assay.4 × 10 5 cells were seeded on 12 well culture plates in DMEM supplemented with 10% FBS for 24 hours.After the cells reached to about 80-90% con uence in a monolayer, a pipette tip was used to make a straight scratch line in the cell monolayer.The cells were incubated for indicated times and treated with BPS at the concentration of 10 − 8 M. (2) Transwell assay.Transwell cell with a pore size of 8.0um was added to a 24-well plate to form upper and lower cells.Next, 600ul of medium containing 10% FBS were added in the lower chamber, 200 portions of serum-free conditioned medium containing 2*104 cells and BPS with different concentration gradients were added in the upper chamber.After 24 hours, wipe the upper chamber cells with a cotton swab, x the invading cells, and count by crystal violet staining.

Statistical analysis
All the analysis were conducted in R software v4.0.0,SPSS v13.0 and GraphPad Prism 8. P-value was two-side and < 0.05 was regarded as statistically signi cant.Student T-test were used to compare the difference between the two groups.

Enrichment analysis of interactive genes affected by bisphenols
On the basis of CTD database, we obtained interactive genes of Bisphenol A, Bisphenol AF, Bisphenol B, Bisphenol C, Bisphenol E, Bisphenol F, Bisphenol S and Bisphenol Z.After summarizing the data, a total of 24041 interactive genes were obtained.Among the genes with interactions, we selected interactive genes with more than 15 interaction counts for further analysis, includes some genes with strong interaction counts INS (3934 interaction counts), ESR1 (771 interaction counts), AR (263 interaction counts), ESR2 (229 interaction counts), MAPK1 (165 interaction counts), MAPK3 (164 interaction counts), LIF (128 interaction counts), PPARG (115 interaction counts), Then we select the interactive genes that demonstrate more than 15 interactions for GO and KEGG analysis to gure out the BPs, CC and MFs.
In terms of BPs (Fig. 2A), response to estradiol, response to drug, response to peptide hormone, regulation of lipid metabolic process and gland development were the most enrichment results.Estradiol is a key hormone which plays an important role in hypothalamic-pituitary-testicular axis regulation, reproductive function and body composition (Luine et al. 2014).Estradiol, ERα, ERβ can limit the growth and spread of prostate cancer cells and prevent prostate cancer from progressing towards castration-resistance forms (Russell et al. 2019).Lipid metabolism may promote the growth of prostate cancer and the resistance to endocrine therapy (Stoykova et al. 2019).
The GO enrichment analysis of interactive genes demonstrated that CC were most associated with vesicle lumen, membrane raft, membrane microdomain, membrane region, secretory granule lumen and cytoplasmic vesicle lumen.Membrane raft is the microdomain which is enriched in cholesterol and sphingomyelin (Sonnino et al. 2013).Prostate cancer contains high levels of cholesterol, which means that larger rafts can be formed in these cell membranes, stimulating signaling pathways to promote tumor growth and progression (Hryniewicz-Jankowska et al. 2019).(Fig. 2B) For MF, the most enriched terms included nuclear receptor activity, ligand-activated transcription factor activity, steroid binding, DNA-binding transcription factor binding, RNA polymerase II-speci c DNA-binding transcription factor binding and receptor ligand activity.As a ligand-activated transcription factor, the androgen receptor (AR) drives prostate cancer and plays an important role in the mainstay treatment for men with metastatic disease (Gabriely et al. 2017).(Fig. 2C) We nally conducted KEGG enrichment analysis by identifying 219 interactive genes to gure out the most relative diseases associated with bisphenol.After screening the rst 20 pathways related to bisphenol, the most enriched pathway was prostate cancer.Besides, bisphenol is also closely related to hepatitis B, reproductive system development and endocrine resistance.Moreover, KEGG pathway enrichment analysis demonstrated that bladder cancer may also associated with bisphenol.(Fig. 2D)

PPI network of key genes associated with prostate cancer and bisphenol
The PPI network was constructed by String (https://www.string-db.org/).Next, the cytoscape 3.8.2 was applied to analyze PPI network and further explore the roles of interactive genes involved in bisphenol.The minimum required interaction score was set as 0.9 and the degree value is correlated with node size.After comprehensive analysis of interactive genes, a total of 219 nodes and 2103 edges were identi ed.Our results revealed that the interactive genes demonstrated strong protein-protein interaction.Some interactive genes, including JUN, CREBBP, AKT1, MAPK1, STAT3, TP53 NCOA1, MAPK3, ESR1, SRC, PPARG, VEGFA, NCOA2 showed more than 30 degrees.(Fig. 3)

Target prediction for bisphenol
The 2D chemical structure of bisphenols were obtained from PubChem.With the help of Pharmmapper database, a total of 796 predicted targets of bisphenols were found.In order to gure out the potential relationship between prostate cancer and bisphenols, we used the KEGG analysis to gure out its underly pathway.We found that AR was considered as key targets of bisphenols, which also play a key role in the progression of prostate cancer.(Fig. 4)

Screening of differential expression genes DEGs
RNA-seq data from the GEO database and TCGA database were applied to compare the DEGs between high-risk group and low-risk group based on gleason score.After data download and preprocessing by using the Perl script, we obtained a total of 213 patients with low-gleason score and 209 patients with high-gleason score.The results demonstrated that 57 downregulated and 22 upregulated genes in the prostate of prostate cancer patients with high-gleason score compared with those in the prostate of prostate cancer patients with low-gleason score.(Fig. 5A) The standard was P < 0.05 and absolute log2FC > 0.5.After screening the prostate cancer-associated genes, venn diagram was conducted to explored the genes associated with prostate cancer and bisphenol.COL1A1 and COL1A2 were de ned as bisphenol-exposure genes associated with prostate cancer.(Fig. 5B)

Function and interaction analysis of interactive genes
In order to gure out the impact of bisphenol-exposure genes associated with prostate cancer on prostate cancer patients, transcriptome data and clinical data of prostate cancer patients were downloaded from TCGA-PRAD.Disease-free survival (DFS) was applied for the measurement of survival analysis.The results showed that prostate cancer patients with high expression of COL1A1 and COL1A2 were related to the poor DFS. (Fig. 6A,B) On the basis of survival analysis of the bisphenol-exposure genes associated with prostate cancer, COL1A1 and COL1A2 showed strong function of the prediction in the prognosis in prostate cancer.Time-dependent ROC curves showed AUC of COL1A1 and COL1A2 are 0.685 and 0.625 respectively.Based on the survival analysis and AUC values, COL1A1 and COLIA2 showed strong value in the prediction of prostate cancer.(Fig. 6C) Gene set variation analysis (GSVA) between low-expression and high-expression groups of COL1A1 and COLIA2 In the previous work, survival analysis and ROC curves have conducted to prove the function of COL1A1 and COL1A2 in prostate cancer and bisphenol-exposure.In order to further explore the potential function of COL1A1 and COL1A2, GSVA of hallmark gene sets were conducted by "limma","stringr" packages.The results demonstrated that collagen trimer, cell cycle and transporter activity were highly activated.(Fig. 7A,B)

Cell proliferation and invasion of prostate cells with BPStreated
The CCK8 assay was applied to gure out the cell proliferation ofboth PC-3 and LNCap cells with the BPS-treated (10 − 8 M).After treated with BPS for 72 hours, we found that BPS signi cantly increased the cell proliferation ability in both PC-3 and LNCap cells compared with the control group (p < 0.05) (Fig. 8A-B).The results of transwell assay showed that the BPS could promote the cell invasion of PC-3 and LNCap cells (Fig. 8C-D).Wound-healing assay demonstrated that BPS signi cantly promoted the migratory capacities of in both PC-3 and LNCap cells(Fig.8E-F).Besides, COL1A1 and COL1A2 in the PC-3 and LNCap cells treated with 10 − 8 M BPS were signi cantly upregulated compared with control group (Fig. 8G-H).

Discussion
As an industrial component with good strength and hardness, thermal stability and resistance to acids and oils, bisphenol has been widely used in industrial production, including materials for the production of phenolic resin, polyacrylate, polyester, epoxy resin and polycarbonate plastics (Mao et al. 2020).However, due to its toxicity and estrogen-like effects, manufacturers start to remove BPA from their products and look for alternative chemicals (Huang et al. 2020).As the two mostly used substitute for BPA, BPS and BPF has been detected in many daily products, including personal care products, paper products and food.BPS and BPF are structural analogues of BPA, therefore, their effects in the physiological system may be similar (Xie et al. 2020).The rodent studies reported that BPS or BPF exposure increased female to male sex ratio, reduced the testosterone in male rodents, damaged the sperm parameters (Thoene et al. 2020).Prostate cancer demonstrated a higher prevalence in western and immigrant population, which implies lifestyle and environmental may be the risk factors (Komura et al. 2018).As a tumor with the high mortality rate in men, it has long been known that prostate cancer is closely related to androgen receptor (Chang et al. 2014).Androgen-receptor signaling is altered in castration-resistant prostate cancer (Liu et al. 2018).EDCs with estrogenic activity is associated with the increased risk of prostate cancer.Cancer stem cells are potential factors for the occurrence and progression of cancer.EDCs with estrogenic activity may in uence the occurrence and progression of prostate cancer by reprogramming and transforming prostate stem and early progenitor cells (Zhang et al. 2018).
In recent years, due to the rapid development of network databases, many disciplines have begun to use databases for forecasting work.As an emerging research tool, analysis using the database can provide us with the correct direction for our subsequent experiments.In this work, in order to gure out the potential relationship between prostate cancer and bisphenol, bioinformatics approach was applied by silico analysis.On the basis of CTD and TCGA databases, we successfully discovered the potential tumor pathway of bisphenol.KEGG enrichment analysis demonstrated that many types of tumors, including prostate cancer, bladder cancer, colorectal cancer and thyroid cancer were highly enriched, which prostate cancer is the most enriched cancer.In order to discover the most related genes between prostate cancer and bisphenol, TCGA and GEO database were applied.Further investigations suggested that several interactive genes were strongly correlated with the progression of prostate cancer.We believe that these two interactive genes (COL1A1 and COL1A2) have predictive effects and may play an important role in promoting the development of prostate cancer induced by bisphenols.By comprehensive analyses of public databases, our work successfully evaluated the associations between bisphenols in the CTD database with the molecular markers of prostate cancer.Our study results provided novel clues for studies on the pathogenetic and risk factors of prostate cancer.
In previous studies, COL1A1 and COL1A2 was considered to be a downstream effector of cytoglobin associated with many common malignancies (Ma et al. 2019).In addition, the increase of COL1A1 and COL1A2 expression in colon cancer is signi cantly related to serous membrane in ltration, lymphatic metastasis and hematogenous metastasis.COL1A1 and COL1A2 can be used as a candidate diagnostic biomarker for colon cancer.Our study found that COL1A1 and COL1A2 can promote the progression of prostate cancer.Our data demonstrated that BPS signi cantly increase the cell proliferation ability and promote cell invasion capacity.COL1A1 and COL1A2 in the PC-3 and LNCap cells were signi cantly upregulated with BPS treated.Our research will help further understand the adverse effects of bisphenols and provide direction to future studies on bisphenols.

Conclusion
As a substitute for BPA, the application of BPS has been increasing over the decades.However, since BPS is considered to be highly toxic to mammals, BPS has attracted more and more attention.In this research, we explored the genes affected by BPS and its potential targets.Our results demonstrated that BPS was closely related to prostate cancer, Further studies proved COL1A1 and COL1A2 played an important role in predicting the prognosis and progression of prostate cancer.Besides, we veri ed that BPS can promote the proliferation and invasion of prostate cancer cell lines.

Figure 2 Circle
Figure 2

Figure 3 Construction
Figure 3