HOTAIR Facilitates Endocrine Resistance in Breast Cancer Through ESR1/miR-130b-3p Axis: Comprehensive Analysis of mRNA-miRNA-lncRNA Network

Introduction

Breast cancer is the most common malignancy and the main cause of cancer death among women worldwide.¹ Seventy-five percent of breast cancer tissues express estrogen receptor-alpha (ERα). Endocrine therapy (ET) is one of the most effective adjuvant therapies for these ER-positive patients,^2–4 including selective estrogen receptor modulators (SERMs), selective estrogen receptor degraders (SERDs) and inhibitors of the enzyme aromatase (AI) converting androgens to estrogens.^4,5 Clinically, more than half of the patients benefit from ET initially, however nearly 40% experience de novo or acquired ET resistance (ETR).⁶ Therefore, it is imperative to explore the potential mechanism of endocrine-resistant breast cancer and to identify novel therapeutic targets.

Previous studies have supported that noncoding RNA including miRNA, long noncoding RNA (lncRNA) and circular RNA (circRNA) are involved in the occurrence and development of multiple tumors.^7,8 A competing endogenous RNA (ceRNA) hypothesis point out that mRNA, miRNA and lncRNA could achieve cross-talk between each other through a regulatory network.⁹ By the use of miRNA response elements, lncRNAs act as ‘sponges’ for miRNAs and result in the miRNAs-regulated mRNA levels alteration. The aberrant regulation of mRNA-miRNA-lncRNA network play key role in a variety of biological processes and molecular mechanism of tumors, including regulating gene transcription and post-transcriptional translation, epithelial-to-mesenchymal transition, signaling pathways, prognostic evaluation and therapeutic targets.^10–13 The deregulation of protein‐coding and noncoding genes in endocrine resistance in luminal breast cancer has reported in several studies.^14–16 However, understanding of the core mRNA-miRNA-lncRNA networks associated with endocrine-resistant breast cancer is still limited.

In this study, bioinformatic analysis were conducted to explore differentially expressed genes (DEGs), differentially expressed miRNAs (DEMs), and differentially expressed lncRNAs (DELs) in long-term estrogen-deprived (LTED) ER-positive breast cancer cells, to establish a mRNA-miRNA-lncRNA regulatory network associated with endocrine resistance, and to mine potential therapy targets for overcoming endocrine resistance.

Materials and Methods

Data Resources and Differentially Expression Analysis

To identify DEGs, DEMs, and DELs in endocrine resistant breast cancer cells, gene expression dataset (GSE120929) and miRNA expression dataset (GSE120930) were obtained by searching the keywords of “breast cancer” and “endocrine resistance” in GEO database (http://www.ncbi.nlm.nih.gov/geo). |log2 fold change (FC)|>2 and p<0.05 was the cutoff criteria. Volcano plots were created to visualize the expression of all genes in the datasets.

Function Enrichment Analysis

Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using Database for annotation, visualization, and integrated discovery (DAVID, https://david.ncifcrf.gov/). p < 0.05 and enrichment score> 1.0 were significant.

Comprehensive Analysis of PPI Network, Identification and Validation of Hub Genes

The protein–protein interaction (PPI) network was constructed using The Search Tool for the Retrieval of Interacting Genes (STRING) database and visualized by Cytoscape 3.6.1. The key modules in PPI network was selected using the cutoff criteria (MCODE score > 3) with the default parameters (degree cutoff = 2, node score cutoff = 0.2, K‐core = 2, and Max depth = 100). Cytoscape CentiScape was used to screen hub genes in the network, according to the degree of centrality. The expression levels of hub genes were further validated between breast cancer and normal samples in GEPIA database. The prognostic values of hub genes including overall survival (OS) and relapse-free survival (RFS) was evaluated in breast cancer patients and breast cancer patients with following systemic treatment (endocrine therapy) in GEPIA and Kaplan–Meier plotter database.

Identification of miRNA Target Genes and lncRNAs

The target genes interacted with miRNAs were predicted based on the miRTarbase database (http://miRTarBase.cuhk.edu.cn/). The upstream lncRNAs of miRNAs were predicted by miRNet database. “Organism-H.sapies”, “Breast cancerous tissues” and “target type-lncRNAs” were considered as selection criteria. FunRich 3.1 software was used to summarize the overlapping genes between the predicted miRNA targets and DEGs, the predicted lncRNAs and DELs.

Construction and Validation of mRNA-miRNA-lncRNA Regulatory Network

The interaction among mRNA, miRNA and lncRNA related to endocrine resistance was constructed according to the lncRNA targets and miRNA targets prediction using the Cytoscape software. The expression levels of miRNAs and hub genes from the network were further validated between breast cancer and normal samples in GEPIA and UALCAN database. The prognostic values of hub genes including OS and RFS was evaluated in breast cancer patients, ER-positive breast cancer patients and breast cancer patients with following systemic treatment (endocrine therapy) in GEPIA and Kaplan–Meier plotter database.

Results

Identification of DEGs, DEMs, and DELs

Dataset GSE120929 and GSE120930 were downloaded to select DEGs, DEMs and DELs in LTED ER-positive breast cancer cells (LTED MCF7) (modeling relapse on endocrine therapy) and MCF7 cells in the presence of estrogen (E2) (modeling a patient at primary diagnosis) (Figure 1A). A total of 769 DEGs including 443 up-regulated and 326 down-regulated genes were identified between LTED cells and MCF-7 cells with E2 (Figure 1B). A total of 33 DEMs, 18 upregulated and 15 downregulated miRNAs, were determined in GSE120930 (Figure 1C). As for DELs, 10 DELs (7 up-regulated and 3 down-regulated lncRNAs) were found in total (Figure 1D).

Figure 1 Identification of significant differentially expressed genes (DEGs), differentially expressed miRNAs(DEMs), and differentially expressed lncRNAs(DELs) in long-term estrogen-deprived (LTED) ER-positive breast cancer cells (LTED MCF7) (modeling relapse on endocrine therapy) and MCF7 cells in the presence of estrogen (E2) (modeling a patient at primary diagnosis). (A) Workflow to identify molecular signature markers associated with breast cancer endocrine resistance from the GEO database. Volcano plot showing the DEGs (B), DEMs (C) and DELs (D) selected from GSE120929 and GSE120930. The blue dots indicating downregulated genes, the red dots indicating upregulated genes.

Functional Analysis for the Significant DEGs

For upregulated significant DEGs, the enriched GO functions are cell adhesion, cell–cell signaling and skeletal system development in the biological process (BP) category, calcium ion binding, calmodulin binding and heme binding in the molecular function (MF) category, and integral component of membrane, plasma membrane and integral component of plasma membrane in the cellular component (CC) category (Figure S1A–C). Besides, Figure S2A indicated that the upregulated significant DEGs were enriched in calcium signaling pathway, protein digestion and absorption and gap junction.

The downregulated significant DEGs were enriched in type I interferon signaling pathway, defense response to virus and response to virus in the BP category, cytokine activity, signal transducer activity and double-stranded RNA binding in the MF category, and plasma membrane, cytosol and extracellular space in the CC category (Figure S1D–F). Additionally, KEGG pathway enrichment analysis indicated that the downregulated DEGs were particularly involved in pentose and glucuronate interconversions, hepatitis C and Jak-STAT signaling pathway (Figure S2B).

Identification and Validation of Hub Genes

313 nodes and 768 edges of upregulated significant DEGs and 257 nodes and 290 edges of downregulated significant DEGs were mapped in the PPI network. The interaction among the top 30 upregulated and downregulated hub genes were visualized by Cytoscape software (Figure S3A–B). Two significant modules in the PPI network were identified using MCODE (Figure S3C and D). The top 10 hub genes and the corresponding node degrees are shown in Table 1 and the top 10 upregulated hub genes were CXCR4, ADCY1, CD44, ESR1, GNAI1, COL1A1, GCGR, CCR1, SOCS3 and PTGER3, and the top 10 downregulated hub genes were STAT1, IRF7, HERC6, IRF9, USP18, ISG15, IFIT1, HERC5, OAS1 and IFIT3.

Table 1 The Top 10 Hub Genes in PPI Networks

The expression levels and prognostic values of those key mRNAs in breast cancer were validated in the GEPIA database. Only 1 hub genes (ESR1) was significantly overexpressed and associated with worse OS of patients with breast cancer (Figure 2A and B). The prognostic values of those hub genes were also evaluated using Kaplan–Meier plotter database. As a result, only one of these hub genes (ESR1) showed both worse OS and RFS in breast cancer patients with following systemic treatment (endocrine therapy) (Figure 2C and D).

Figure 2 The expression and prognostic values of the hub genes from GEPIA database and Kaplan-Meier plotter databases. (A) The expression level of ESR1 validated in GEPIA database (*P=7.79E-13). (B) The overexpression of ESR1 was associated with worse OS of patients with breast cancer (P=0.03). ESR1 shows both worse OS (C) (P=0.012) and RFS (D) (P=0.00066) in breast cancer patients with following systemic treatment (endocrine therapy).

Integration of mRNA-miRNA-lncRNA Regulatory Network

To explore the role of miRNAs and lncRNAs in endocrine resistant cells, the mRNA-miRNA-lncRNA interaction network was constructed based on 60 mRNA nodes, 6 miRNA nodes and 3 lncRNA nodes (Figure 3A). Six overlapping DEMs (hsa-miR-130b-3p, hsa-miR-196a-5p, hsa-miR-23b-3p, hsa-miR-10a-5p, hsa-miR-195-5p and hsa-miR-370-3p) were identified in dataset GSE120930. The target genes of hsa-miR-130b-3p, hsa-miR-196a-5p, hsa-miR-23b-3p, hsa-miR-10a-5p, hsa-miR-195-5p and hsa-miR-370-3p overlapping with DEGs are listed in Table 2. The connections between DELs and lncRNAs predicted by these six miRNAs in miRNet database are shown in Table 3. The interaction of the 60 overlapping mRNA nodes were presented in the PPI network (Figure 3B). A significant module containing 3 nodes and 3 edges was identified. The hub genes were ABCG2, ESR1 and GJA1 (Figure 3C). In the network, we found ESR1, the target genes of hsa-miR-130b-3p belonging to the top 10 hub genes, revealing the potential role of ESR1/miR-130b-3p/HOTAIR in regulating endocrine resistance in breast cancer (Figure 3B and C).

Table 2 Overlapping Genes Between the miRNA Targets and Overlapping DEGs

Table 3 Overlapping lncRNAs Between the miRNA Upstream Targets and Overlapping DELs

Figure 3 The construction of mRNA-miRNA-lncRNA regulatory network. (A) The network consists of 60 mRNA nodes, 6 miRNA nodes and 3 lncRNA nodes. Ellipses: mRNA, round rectangles: miRNA, triangles: lncRNAs, the red nodes: up-regulation, the blue nodes: down-regulation. (B) PPI network containing the 60 overlapping mRNA nodes. (C) Network of three hub genes extracted from the PPI network.

GO analysis revealed that the 60 DEGs in network were enriched in 46 GO terms. The top enriched GO terms are response to hypoxia in BP terms, protein binding in MF terms, integral component of membrane in CC terms (Figure S4A–C). The top KEGG pathways were Pathways in cancer, Proteoglycans in cancer, HIF-1 signaling pathway (Figure S4D).

The Validation of Genes in mRNA-miRNA-lncRNA Network

The expression of miRNAs in the network and the prognostic values were assessed in breast cancer patients and breast cancer patients with following systemic treatment (endocrine therapy) in the UALCAN database and Kaplan–Meier plotter database, respectively. Three of these overlapping miRNAs (hsa-miR-130b-3p, hsa-miR-10a-5p and hsa-miR-195-5p) were aberrantly expressed in breast cancer and the expression was associated with tworse OS of breast cancer patients with following systemic treatment (endocrine therapy) (Figure 4). As for lncRNAs, HOTAIR is the only one significantly expressed genes from the network in breast cancer (Figure 5A). The high expression of HOTAIR was significantly related to worse OS in breast cancer patients and worse RFS in breast cancer patients with following systemic treatment (endocrine therapy), lightly related to worse OS in ER-positive patients (Figure 5B–D). Based on the hub genes selected from 60 overlapping mRNA nodes in the network (ABCG2, ESR1 and GJA1), ESR1 and ABCG2 showed aberrant expression related to OS in breast cancer patients (Figures 2 and S5A and B). High expression of ESR1 and low expression of GJA1 showed both worse OS and RFS in breast cancer patients with following systemic treatment (endocrine therapy) (Figures 2 and S5C and D).

Figure 4 The expression and prognostic values of miR-130b (A and B), miR-10a (C and D) and miR-195 (E and F) from mRNA-miRNA-lncRNA regulatory network in UALCAN and Kaplan-Meier plotter databases.

Figure 5 The expression and prognostic values of lncRNAs from mRNA-miRNA-lncRNA regulatory network. (A) The expression level of HOTAIR validated in GEPIA database (P=5.93E-33). The overexpression of HOTAIR was related to worse OS in breast cancer patients (B) (P=0.038), lightly related to worse OS in ER-positive patients (C) (P=0.052) and related to worse RFS in breast cancer patients with following systemic treatment (endocrine therapy) (D) (P=0.041) in Kaplan-Meier plotter databases.

Discussion

ET remains the mainstream adjuvant treatment on hormone receptor-positive breast cancer, however, de novo or acquired ETR is a major limitation of treatment. It is essential to explore effective biomarkers predicting ETR and find alternative treatment on endocrine resistant tumors.

In our study, 769 DEGs, 33 DEMs, and 10 DELs were identified in the endocrine‐resistant cell by mining the datasets from the GEO database in order to set up ETR related mRNA-miRNA-lncRNA network including 60 mRNA s, 6 miRNAs and 3 lncRNAs. The role of noncoding RNAs in breast cancer and endocrine resistant breast cancers have been recently reviewed.^15,17,18 Overexpression of lncRNA breast cancer anti-estrogen resistance 4 (BCAR4) in tamoxifen sensitive cells partly arrogates the anti-proliferative effects of tamoxifen.¹⁹ BCAR4 is also considered as a biomarker for increased cancer invasiveness and tamoxifen resistance.²⁰ Long non-coding RNA UCA1 promoted tamoxifen resistance in breast cancer cells by regulating HIF1α via sponging miR-18a.¹⁴ lncRNA HOTAIR upregulation promotes ligand independent ER activities and contributes to tamoxifen resistance in breast cancer.¹⁵

Additionally, KEGG pathway enrichment analysis indicated that the enriched pathways primarily involved pathways in cancer, proteoglycans in cancer, HIF-1 signaling pathway. Besides, GO analysis and pathway analysis revealed that the significant DEGs were mostly enriched in response to hypoxia. It has been documented that hypoxia-inducible factors, HIF-1 play key roles in doxorubicin resistance of breast cancer.¹⁶

Subsequently, the target of hsa-miR-130b-3p were found in the top 10 hub genes in 769 DEGs, revealing the vital role of ESR1/miR-130b-3p/HOTAIR in regulating endocrine resistance. The prognostic values of the genes from the mRNA-miRNA-lncRNA network were validated in multiple databases. Eventually, ESR1/miR-130b-3p/HOTAIR were not only differentially expressed in patients with breast cancer but also significantly correlated with the prognosis of breast cancer patients with following systemic treatment (endocrine therapy).

miR-130b acts as an onco-miRNAs and promotes the progression of cancer.^21,22 Meanwhile, the overexpression of miR-130b is a risk factor for poor prognosis in prostate cancer patients.²³ On the contrary, few studies also identified that downregulation of miR-130b plays a role in stimulating progression and angiogenesis of prostate cancer.²⁴ In addition, miR-130b is a target miRNA of CASC15 and CASC15/miR-130b axis might be a novel therapy for bladder cancer and non-small cell lung cancer patients.²⁵ It was reported that miR-130b was downregulated in multidrug resistant ovarian cancer cells and restoration of miR-130b expression could sensitize cells to anticancer drugs.²⁶

HOX transcript antisense RNA (HOTAIR) has been identified as the first lncRNA correlated to poor prognosis of breast cancer, taking part in EMT and maintenance of breast cancer stem cells (bCSCs).^27,28 Based on the Estrogen Responsive Elements (EREs) elements in HOTAIR promoter, estradiol regulate HOTAIR expression in ER-positive breast cancer cells.²⁹ HOTAIR is reported to be significantly upregulated in tamoxifen-resistant breast cancer cells, contributing to tamoxifen resistance by interacting with ER.³⁰

Among the overlapping DEGs, ABCG2 and GJA1 were also screened as hub genes. The expression of ABCG2 was reported to be associated with triple-negative tumors.³¹ The ABCG2 signaling was involved in autophagy-promoting drug resistance in breast cancer stem cells.³² A recent study reported GJA1 is dysregulated in breast cancer and luminal tumors with high levels of GJA1 mRNA were associated with a better prognostic.^33,34

According to the above analysis results, mRNA-miRNA-lncRNA regulatory network was built up. This is a novel mRNA-miRNA-lncRNA regulatory network reported in endocrine resistance. This study provides some powerful evidence for molecular mechanism of endocrine resistance. However, there are a few limitations in our study: the whole set of DEGs was extracted from cell lines and not from tumor tissues; the regulatory network is related to endocrine therapy resistance without specifying which category or which specific medicine. Further study will validate these findings by experiment on cell lines and human tissues.