Genomic Understanding Reveals the Important Role of FGFR2 as Paeoniflorin Target for Circumventing Breast Cancer Resistance to Tamoxifen.

OBJECTIVES: Paeoniflorin (PF), a compound found in Paeonia lactiflora and Paeonia suffruticosa, has anticancer potential, particularly in inhibiting migration and invasion, the resistant cancer cells hallmarks. To date, the mechanism of overcoming tamoxifen resistance in breast cancer is not yet elucidated. This research aims to explore the potential target of PF as a co-treatment for circumventing breast cancer resistance to tamoxifen with a genomic understanding-bioinformatics.
METHODS: Microarray data originating from GSE67916 and GSE85871 in the NCBI GEO database was analyzed to obtain differentially expressed genes (DEGs). Further analyses were performed on DEGs using the DAVID v6.8, STRING-DB v11.0, the Cytoscape, and cBioportal. Gene expression analysis validation in breast cancer cells and tamoxifen-resistant breast cancer cells was accomplished using GEPIA and ONCOMINE databases. Survival rate analysis of selected genes was conducted using Kaplan-Meier.
RESULTS: We obtained 175 DEGs from the two samples (tamoxifen-resistant and paeoniflorin-treated). DEG involves in 70 biological processes, 26 cellular components, and 18 molecular functions, and three pathways relevant to breast cancer. The PPI network analysis and hub genes selection obtained 10 genes with the highest degree scores. Genetic changes for selected genes, including IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 were found from 0.5% to 7% of the case population per patient case. Additional analysis using cBioportal revealed FGFR signaling pathway through Ras is important for the PF mechanism in circumventing breast cancer resistance to tamoxifen. ONCOMINE and GEPIA analysis emphasized the importance of selected genes in the tamoxifen-resistance mechanism.
CONCLUSION: PF has potential to be used as a co-treatment for circumventing breast cancer resistance to tamoxifen by targeting FGFR2 signaling, but further validation is needed.

Introduction

Breast cancer is responsible for 2.1 million cancers, the fifth leading cause of cancer related death worldwide in 2018 (Dolatkhah et al., 2020). Approximately 80% of breast cancers are estrogen receptor (ER)-positive type. Thus, endocrine therapy is considered a complement to surgery in most patients with breast cancer (Lumachi et al., 2015). Tamoxifen (Figure 1A) is a central component of ER-positive breast cancer treatment (Chang, 2012). As a selective ER modulator and a triphenylethylene derivative, it acts as an agonist in the uterine but an antagonist in the breast tissue (Hultch, 2018). However, more than half of advanced ER-positive breast cancers are intrinsically resistant to tamoxifen, and about 40% of them develop resistance during treatment (Hultch et al., 2018). Crosstalk between the ER and a set of growth factors receptor, including HER2, EGFR, FGFR, IGFIR, and impaired PI3K/PTEN and NFκB activation are among factors that responsible for tamoxifen resistance (Ali et al., 2016). Therefore, it is necessary to develop compounds that overcome the resistance to tamoxifen.

Figure 1

(A) Chemical Structure of Tamoxifen, (B) Chemical Structure of Paeoniflorin, (C) A Venn diagram of the DEGs between tamoxifen-resistant and Paeoniflorin (PF)-treated MCF-7 breast cancer cells

Paeoniflorin (PF, Figure 1B) is a promising adjunct therapy to curb tamoxifen resistance in breast cancer. PF is a pinan monoterpene glycoside with various bioactivities (Xiang et al., 2020). The anticancer activity of PF has been comprehensively studied in many types of cancer, including liver, gastric, breast, lung, pancreatic, colorectal, glioma, bladder, and leukemia (Xiang et al., 2020). The PF remarkably impeded MCF-7 breast cancer cell proliferation in a time- and concentration-dependent manner (Zhang et al., 2017). Besides, PF has been proven to prevent hypoxia-induced migration and invasion in MCF-7 cells (Zhou et al., 2016), two hallmarks of resistant cells (Fares et al., 2020). Therefore, PF has the potential to overcome tamoxifen resistance in breast cancer. Collectedly, it is necessary to conduct further studies on how PF works more broadly in circumventing breast cancer resistance to tamoxifen.

Bioinformatics analysis accelerates the drug targets identification and the drug candidates screening and refinement. This analysis also facilitates the side effects characterization and predicts drug resistance mechanisms (Xia, 2017). This study aims to explore the potential target of PF as a co-treatment for circumventing breast cancer resistance to tamoxifen with a genomic understanding-bioinformatics using several public databases, including GEO2R, DAVID v6.8, STRING-DB v11.0, Cytoscape 3.8.2, CytoHubba, and cBioPortal. The approach would facilitate and fasten the identifying potential target genes process and the PF mechanism in circumventing tamoxifen resistance in breast cancer.

Materials and Methods

Data mining

The microarray data were obtained from the NCBI Gene Expression Omnibus (GEO) database with ID GSE67916 for Tamoxifen-resistant ER+ breast cancer cells and GSE85871 for MCF7 cells treated with paeoniflorin (Lv et al., 2017). The gene expression in both gene series expressions were profiled using the microarray technology with Affymetrix Human Genome U133A 2.0 (Santa Clara, CA). The microarray data of the two GSE codes were further analyzed using the GEO2R to identify differentially expressed genes (DEGs) (Zhang et al., 2020). The adjusted p-value < 0.05 and Iog fold change > 1 were used to select the significant (Hermawan and Putri, 2020a). The two samples derived from the GSE code were then analyzed using Venny 2.1.0 to obtain a Venn diagram (Liu et al., 2019). From the Venn diagram, it is obtained that DEGs are shared by both samples.

Gene ontology and KEGG pathway enrichment analysis

The DEGs from the previous step was analyzed for gene ontology (GO) enrichment analysis consisting of biological processes, molecular functions, and cellular component; and KEGG pathway enrichment analysis using the database for annotation, visualization, and integrated discovery (DAVID) v6.8) (Wulandari et al., 2020). The value of p < 0.05 was selected as the significant value (Hermawan et al., 2021a).

Proten-protein interaction (PPI) network and hub gene selection

The PPI network analysis of the DEGs shared by the two GSE codes was executed using STRING-DB (Karimizadeh et al., 2019) with medium confidence (score = 0.400) and then visualized using Cytoscape software (Hermawan and Putri, 2020b). Comprehensive relationship analysis between nodes in the PPI was carried out using the Maximal Clique Centrality (MCC) function, the most effective method for finding hub genes in the PPI network (Li et al., 2020).

Analysis of genetic alterations of selected target genes

According to the previous study, the genetic alterations of selected target genes were analyzed using cBioPortal (Hermawan et al., 2021b). The breast cancer study with the highest genetic alterations was chosen for further connectivity analysis with the significant value at p < 0.05.

Gene expression analysis validation in non-resistant and tamoxifen-resistant breast cancer cells

Validation analysis of IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 was performed using the GEPIA database (http://gepia.cancer-pku.cn). The GEPIA database analysis was carried out to determine the difference in expression between cancer cells and normal cells based on TCGA and GTEx data (Tang et al., 2017). The default settings were used except for the BRCA with logFC > 1 and p-value < 0.01 dataset.

Validation of gene expression analysis was conducted according to Hermawan et al., (2020). Analysis was carried out on IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 genes in the ONCOMINE database (https://oncomine.org) (Rhodes et al., 2004). The settings used in the ONCOMINE database were cancer type analysis vs. cancer analysis, breast cancer, the clinical outcome of recurrence and metastatic status, and molecular subtype of the estrogen receptor. The analysis was performed to determine the expression level in the recurrence group and the nonrecurrence group after the tamoxifen treatment.

Survival rate analysis

Survival rate analysis of IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 was conducted using Kaplan–Meier survival rate plots in the Kaplan–Meier plotter database (http://kmplot.com) (Gyorffy, 2021). Kaplan–Meier plotter provides a correlation between survival and gene expression, that can be used for clinical cohorts’ evaluation (Gyorffy, 2021). The survival rate analysis was executed using breast cancer mRNA expression datasets. The default settings were used in the analysis with the exception of ER-positive only samples and restricted endocrine therapy only with tamoxifen. The samples number in the analysis is shown in each figure, and the p-value used was < 0.05.

Results

Data collection and processing

The objective of this study was to identify potential target genes and molecular mechanisms of paeoniflorin compounds in circumventing tamoxifen resistance in breast cancer. There was a good GSE67916 and GSE86871 data distribution (Supplementary Figure 1). A total of 3388 DEGs were obtained from the GSE67916 sample of gene expression analysis in Tamoxifen-resistant ER+ breast cancer cell lines, with details: 976 DEGs in the form of genes with increased expression and 2412 DEGs in the form of genes with decreased expression (Supplementary Table 1). Then 1210 DEGs were obtained from the GSE85871 sample of gene expression profiles of MCF-7 cells treated with traditional Chinese medicine components, with details: 857 DEGs of genes with increased expression and 353 DEGs of genes with decreased expression (Supplementary Table 2). Then, 175 DEGs owned by the two samples (tamoxifen-resistant and paeoniflorin-treated) were obtained from the Ven diagram by Venny 2.1.0 (Figure 1C, Supplementary Table 3). The DEGs of these two samples were analyzed further (Figure 1C).

Gene ontology and KEGG pathway enrichment analysis

Enrichment analysis of the Ontology Gene (GO) pathway and KEGG pathway was carried out to study DEGs’ biological function and molecular mechanism. The GO and KEGG pathway analyses were performed with DAVID v6.8 using a cut-off p-value <0.05 based on three criteria: biological processes, cellular components, and molecular functions. The results showed that DEG involves in 70 biological processes, 26 cellular components, and 18 molecular functions (Figure 2, Supplementary Table 4). Three KEGG pathway analyses with DAVID revealed three pathways relevant to breast cancer (Table 1).

Figure 2

GO Enrichment Analysis on DEGs with DAVID, Consist of Biological Process, Molecular Function, and Cellular Component

Table 1

KEGG Pathway of DEGs Relevant to Breast Cancer

Pathway	P-value	DEGs
hsa04630:Jak-STAT signaling pathway	0.012149	IFNB1, STAT2, PTPN11, PIK3CB, IL6ST, PRLR, IRF9
hsa05200:Pathways in cancer	0.014094	KITLG, CDK6, COL4A3, TPR, BCL2, CTNNA1, GNAS, LAMB1, PIK3CB, SKP2, FGFR2, PML
hsa04151:PI3K-Akt signaling pathway	0.037994	RBL2, KITLG, CDK6, IFNB1, COL4A3, BCL2, LAMB1, PIK3CB, PRLR, FGFR

PPI network and hub gene selection

Protein-protein Interaction (PPI) has a central role in running the organisms’ cellular systems. Problems with PPIs can be a cause or indication of disease (Kuzmanov & Emili, 2013). The PPI network analysis was conducted to study the interactions between DEGs. PPI network analysis showed that the PPI was made up from 188 proteins with a confidence level of 0.4 (medium confidence), consisting of 173 nodes, 364 edges, PPI enrichment p-value: < 1.0e-16 and avg. local clustering coefficient: 0.402 (Figure 3A).

Figure 3

(A). PPI network of the DEGs, as constructed by STRING-DB v11.0 (B). PPI network between hub genes and surrounding DEGs, as analyzed by CytoHubba plugin on Cytoscape 3.8.2

Further analysis with the CytoHubba plugs in to the Cytoscpape 3.8.2 application to select the hub gene, a gene that has a high level of intra-module connectivity (Zu et al., 2019). Present study defined the hub gene as the ten genes with the highest connectivity in the obtained PPI network. A comprehensive analysis of the relationship between nodes in the PPI was carried out using the MCC function, the most effective method for finding hub genes in the PPI network (Li et al., 2020). The top 10 genes with the highest degree scores (Fig. 3B, Table 2), in a particular order, are OAS1, OASL, IRF9, MX1, OAS3, IFI27, IFIT1, IFI6, IFI44, and STAT2. Different colors represent different degree scores of the hub genes. The red color indicated a higher degree score (Xiao et al., 2018).

Table 2

Selection of hub Genes, as Analyzed by CytoHubba plugin of CytoScape Based on MCC

Rank	Name	Score
1	OAS1 (2'-5'-oligoadenylate synthetase 1)	1.71E+08
2	OASL (2'-5'-Oligoadenylate Synthetase Like)	1.71E+08
3	IRF9 (Interferon regulatory factor 9)	1.71E+08
4	MX1 (Myxovirus resistance 1)	1.71E+08
5	OAS3 (2'-5'-Oligoadenylate Synthetase 3)	1.71E+08
6	IFI27 (Interferon Alpha Inducible Protein 27)	1.71E+08
7	IFIT1 (Interferon Induced Protein With Tetratricopeptide Repeats 1)	1.71E+08
8	IFI6 (Interferon alpha-inducible protein 6)	1.71E+08
9	IFI44 (Interferon Induced Protein 44)	1.27E+08
10	STAT2 (Signal Transducer and Activator of Transcription 2)	1.23E+08

Analysis of genetic alterations of selected genes

Among the DEGs, IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 were selected to explore their genome changes across breast cancer studies using cBioPortal. FGRF2, BCL2, and CDK6 were preferred among the DEGs of the KEGG pathway enrichment (cancer signaling pathway and PI3K-Akt). IFNB1 was chosen based on Jak-STAT and PI3k-Akt pathways. STAT2 was selected from the hub gene and DEGs from the Jak-STAT pathway, while OAS1 was the hub gene in the first order. A study, BRCA (INSERM) 2016, showed the second-highest genetic alteration among breast cancer studies and had more variation of DEGs-associated genetic changes than the study that ranked first (Figure 4A). Genetic changes for each target gene were found from 0.5% to 7% of the case population per patient case, including STAT2 (0.5%), OAS1(1.4%), BCL2 (2.8%), FGFR2 (3 %), IFNB1 (3%), and CDK6 (7%) (Figure 4B). In addition, most of the gene changes include amplification and deep deletion. Additional analysis using cBioportal revealed FGFR signaling pathway through Ras is important for the PF mechanism in circumventing breast cancer resistance to tamoxifen (Figure 4C).

Figure 4

(A) Genetic Alterations of DEGs (IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2) presented on a chart of 14 breast cancer studies (B) Oncoprint cBioportal in the 2016 BRCA study (INSERM) : Genetic alteration of selected DEGs (IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2) in the patient population. (C). Results of DEGs analysis by cBioportal in the pathway section related to gene alteration. The results showed that DEGs are involved in the FGFR2 signaling pathway through Ras

Validation of gene expression analysis in non-resistant and tamoxifen-resistant breast cancer cells

Validation analysis with GEPIA showed that the mRNA expression in non-resistant breast cancer cells and normal cells were not significantly different except that OAS1 was expressed significantly higher in tumor cells than in normal cells (p-value < 0.01) (Figure 5A). ONCOMINE analysis showed that CDK6 and STAT were overexpressed in the tamoxifen-resistant group or recurrent group (Figure 5B). However, the expression of IFNB1 did not differ between recurrent and nonrecurrent groups. Furthermore, FGFR2 and OAS1 were expressed slightly lower than nonrecurrent, but it is not apparent. In this analysis using ONCOMINE, the BCL2 was excluded as it is not available in the dataset used.

Figure 5

(A). Validation analysis of gene expression of IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 in non-resistant breast cancer cells and normal cells using GEPIA. (B). Gene expression analysis of IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 in recurrent and non-recurrent group after treated with tamoxifen using ONCOMINE

Survival rate analysis

Kaplan–Meier curves showed that the survival rate of low mRNA expression of CDK6 was significantly lower than that of the high mRNA expression group (p = 0.026) (Figure 6). The survival rate of STAT2 between the low expression and high expression groups was not significantly different, but both had a fairly high survival rate, above 80%. Overall survival rates of the high-expressed mRNA of IFNB1 (p = 0.72), FGFR2 (p = 0.17), OAS1 (p = 0.14), and BCL2 (p = 0.74) were worse than those of the lower mRNA expression, but the results were not significant.

Figure 6

Survival Rate Analysis of IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2 in ER+ Breast Cancer and Tamoxifen-Treated Only, Using Kaplan–Meier Curve in KMPlotter

Discussion

The results of this study using genomic understanding-based bioinformatics analysis highlighted the potential target genes, namely: IFNB1, CDK6, FGFR2, OAS1, BCL2, and STAT2. Furthermore, by using cBioportal, it was found that the PF mechanism in circumventing breast cancer resistance to tamoxifen is via the FGFR signaling pathway through Ras. PF circumvent breast cancer resistance to tamoxifen by targeting FGFR2 via inhibition of its downstream signaling or inhibition of its transcriptional activity through Notch signaling (Figure 7). The FGFR pathway is initiated when the FGFR receptor or Fibroblast Growth Factor Receptor binds to its ligand, Fibroblast Growth Factor, changing the FGFR conformation, resulting in dimerization and activation of FGFR. Activated FGFR phosphorylates FRS2a and FRS2a, which subsequently bind to the adapter containing the SH2 domain, Grb2. Grb2 further bind to SOS, GAB1, and Cbl via its SH3 domain to activate Ras/Raf/MAPK, including ERK MAPK, p38 MAPK, and JNK MAPK (Xie et al., 2020) and its downstream protein, including cyclin-dependent kinases (CDKs) and BCL2 (Chambard et al., 2007; Modi et al., 2012; Wang et al., 2015).

Figure 7

Proposed Molecular Mechanisms of PF in Circumventing Breast Cancer Resistance to Tamoxifen. PF inhibits the FGFR2 pathway through its downstream signaling. PF also inhibits the transcriptional activity of FGFR2 through inhibition of Notch signaling

The RAS protein serves as the center of a primary intracellular signaling pathway that regulates cell growth, motility, angiogenesis, and immune evasion (Sever and Brugge, 2015). Oncogenic RAS mutations constitutively increase the interaction of breast cells with the basement membrane, induce a phenotype independent of anchorage, invasion, tumorigenic potential, TGF-β and IGF-1 secretion, EGFR activation, mitogen-activated protein kinase (MAPK), and estrogen insensitivity. (Galiè, 2019). Studies have also shown that the MAPK pathway facilitates tamoxifen resistance (Eanes and Patel, 2016). ERK phosphorylation and increased activity are associated with endocrine resistance and decreased survival in breast cancer patients. One possible mechanism is that ERK phosphorylates Ser 118 in the ER, leading to ligand-independent ER activation (Haagenson and Wu, 2014).

In a study conducted by Guo et al., (2017), PF treatment hindered VSMC proliferation by decreasing the expression of four proteins involved in the G0/G1 and G1/S transitions: cyclin D1, cyclin E, CDK4 and CDK2. Furthermore, PF increases the expression of p21, which is known to form a heterotrimetric complex with the cyclin-CDK complex, thereby inhibiting its activity (Guo et al., 2017). PF also exhibited antitumor activity in gastric cancer cells by inhibiting cell proliferation and stimulating apoptosis by reducing the expression of PI3K, p-Akt and Akt (Zheng et al., 2015). Besides, the effect of PF in inhibiting the Akt pathway is also known to occur in primary human liver sinusoidal endothelial cells (HHSECs) (Gong et al., 2015). Additionally, PF is known to have effects on proteins associated with apoptosis. Research conducted by Hu et al. (2018) on mouse chondrocytes shows that PF reduces the expression of mRNA and protein Bcl-2, a protein involved in the apoptosis. Future studies on the FGFR2 and its downstream, including PI3K/Akt, Ras, and MAPK upon PF treatment in tamoxifen-resistant breast cancer cells are needed to the present study’s results.

The second proposed mechanism of PF in circumventing breast cancer resistance to tamoxifen is via inhibition of FGFR2 transcriptional activity through Notch signaling. In invasive breast cancer, increased expression of the Notch signaling pathway components has been reported, including Jagged1-2, Dll1, Dll3, and Dll4, Notch receptors, Hes and Hey target genes, and accumulation of NICD (Wang et al., 2010). Notch signaling is known to regulate many cellular processes, including proliferation, apoptosis, angiogenesis, hypoxia, cancer stem cell activity, EMT, and metastasis (Acar et al., 2016). Crosstalk between estradiol and Notch signaling has a significant role in human breast carcinogenesis and angiogenesis (Kontomanolis et al., 2018). NOTCH-1 activates ERα-dependent transcription in the presence or absence of estradiol. Therefore, Notch signaling regulates estrogen signaling in breast cancer cells, and the combination of antiestrogens and Notch inhibitors may be more effective in ERα(+) breast cancer (Kontomanolis et al., 2018). The Notch-1 target gene, SOX9, is expressed in multi-potent cells in the liver. Overexpression of SOX9 is frequently detected in patients with liver cancer (Huang et al., 2019). SOX9 is a transcription factor (Lefebvre and Ginzberg, 2018) and its overexpression is frequent in breast cancer subtypes (Jana et al., 2020).

SOX9 is elevated in breast cancer patients after endocrine therapy failure (Domeneci et al., 2018). Recently, mouse genetic studies have revealed a feed-forward loop on early pancreatic FGF10/FGFR2b/SOX9 progenitors to sustain their proliferation. SOX9 acts as a link between FGF10 and Notch in progenitor proliferation. (Li et al., 2015). The SOX9 ablation resulted in significant downregulation of FGFR2 and FGFR4, which indicates that SOX9 is required to maintain the expression of FGFR2 and FGFR4 (Seymor et al., 2012). A previous study by Zhang et al., (2018) showed that using different concentrations of PF significantly decreased NOTCH-1 and HES-1 mRNA and protein expression in MCF-7 breast cancer cells. Overall, the axis between Notch signaling-FGFR2 and PF treatment in circumventing breast cancer to tamoxifen are warranted for future studies.

The advantage of using the genomic understanding-based bioinformatics in this study is that it can accelerate the identifying the target of a compound for a drug candidate of a specific disease. In addition, in terms of duration, time needed is relatively shorter than other techniques such as screening in wet laboratory experiments. However, the results of this study, in which targets and molecular mechanisms of PF in circumventing breast cancer resistance to tamoxifen still need to be validated within vitro, in vivo, toxicity test, and clinical trials.

In Conclusion, PF circumvent breast cancer resistance to tamoxifen by targeting FGFR2 via inhibition of its downstream signaling or inhibition of its transcriptional activity through Notch signaling. The results of this study need further confirmation for the development of PF as drugs for circumventing breast cancer resistance to tamoxifen.

Author Contribution Statement

Data collection, analysis and interpretation of data, and writing draft: Moh. Arikus Wibowo and Eri Prasetyo Nugroho. Study concept and design, supervision, funding acquisition, editing, and final approval of manuscript: Adam Hermawan.

Funding statement

This study was supported by Rekognisi Tugas Akhir (RTA) year 2021 Universitas Gadjah Mada, SK Rektor UGM No.1185/UN1.P.III/SK/HUKOR/2021.

Availability of data

The data that support the findings of this study are available in the supplementary files and available in http://ugm.id/PFRTA

Conflict of interest

The authors declare no conflict of interest.