Abstract
Breast cancer is the most frequent malignant tumor in women, and the estrogen receptor (ER) plays a vital role in the vast majority of breast cancers. The purpose of the present study was to identify the significant genes regulated by ER in ER-positive breast cancer and to explore their expression pattern changes when tamoxifen or fulvestrant resistance occurs. For this purpose, the gene expression profiles GSE11324, GSE27473, and GSE5840 from the Gene Expression Omnibus database were used, which contain gene expression data from MCF7 cells treated with estrogen, MCF7 cells with silencing of ER, and tamoxifen- and fulvestrant-resistant MCF7 cells treated with estrogen (17β-estradiol), respectively. Differentially expressed genes (DEGs) between the treatment group and negative control were identified and subjected to pathway enrichment and protein–protein interaction (PPI) analyses. There were 230 DEGs in common among the three datasets, including 160 genes positively regulated by ER and 70 genes negatively regulated by ER. DEGs mainly showed enrichment for pathways in cancer, progesterone-mediated oocyte maturation, RNA transport, glycerophospholipid metabolism, oocyte meiosis, platelet activation, and so on. PPI network and modular analysis selected three significant clusters containing 19 genes. A total of 44 genes were involved in Kyoto Encyclopedia of Gene and Genome pathway results or PPI modular analysis, and 16 of them were found to correlate with relapse-free survival in patients with ER+/human epidermal growth factor receptor 2-negative breast cancer who had undergone endocrine therapies only. Some of the genes’ expression patterns were different among wild-type, tamoxifen-resistant, and fulvestrant-resistant MCF7 cells such as DDX18, ANAPC7, MAD2L1, RSL1D1, and CALCR, etc., indicating different resistance mechanisms and potential prognostic markers or therapeutic targets for fulvestrant- or tamoxifen-resistant breast cancer.
Introduction
Breast cancer is the most frequently diagnosed cancer and the leading cause of cancer-related death in women worldwide (). Based on current clinical findings and basic medical research, estrogen receptor (ER) exists in the majority of breast cancers and has a profound influence on the occurrence and development of breast cancer (). Endocrine therapies, such as tamoxifen, fulvestrant, and aromatase inhibitors, which target estrogen or ER, are useful in breast cancer treatment (; ). However, nearly half of ER-positive breast cancer patients will eventually fail one or more of these endocrine interventions; the incidence of endocrine resistance is relatively high (). Therefore, more reliable endocrine therapies with efficacy- and drug resistance-related biomarkers should be explored. These biomarkers could not only be used to convey prognostic information, but may also offer breakthroughs in overcoming endocrine resistance.
Gene chips, which have been used for more than a decade, are a quick and reliable technique for exploring differentially expressed genes (DEGs) (), and large quantities of chip data have been produced and stored in public databases. These databases represent a valuable achievement that can be retrospectively analyzed based on new perspectives and approaches. Several bioinformatics studies on malignant tumors have been published in recent years, which provide new methods to uncover the underlying mechanisms of the development and progression of different types of cancers ().
In the present study, we retrieved the gene expression profiles GSE11324, GSE27473, and GSE5840 from the Gene Expression Omnibus (GEO) as the primary research datasets. These files contain gene expression data from MCF7 cells (ER-positive breast cancer cell line) treated with estrogen, MCF7 cells with silencing of the ER, and tamoxifen- and fulvestrant-resistant MCF7 cells treated with estrogen (17β-estradiol). These datasets were utilized to identify the DEGs between the treatment group and negative control, followed by gene ontology (GO), pathway enrichment, and protein–protein interaction (PPI) analyses. Next, the interested DEGs were interrogated for prognostic information associated with the relapse-free survival (RFS) of patients with ER-positive and human epidermal growth factor receptor 2 (HER2)-negative breast cancer who had undergone endocrine therapies only, which can be the best to reflect the effect of endocrine therapies and whether patients develop endocrine resistance. Finally, we analyzed these RFS-related genes in GSE5840 gene expression profiles, which contain the differential expression gene information of tamoxifen- and fulvestrant-resistant MCF7 cells treated with estrogen. We wanted to find out significant genes regulated by ER in ER-positive breast cancer and explore their expression pattern changes when tamoxifen or fulvestrant resistance occurs. These genes may reflect the functional status of ER in endocrine-resistant breast cancer cells and promise therapeutic targets for tamoxifen or fulvestrant resistance breast cancer.
Materials and Methods
Microarray Data
The gene expression profiles used in this study were obtained from the NCBI-GEO, a free public data repository of microarray and other genomic data. Three ER-positive breast cancer cell line gene expression profiles, GSE11324, GSE27473, and GSE5840, were chosen as the primary research datasets. GSE11324 contains gene expression data from MCF7 cells that were stimulated with estrogen for 0, 3, 6, or 12 h. All experiments were performed in triplicate (). We compared the gene expression data between 0 and 3, 0 and 6, and 0 and 12 h, respectively. There are nearly consistent trends in genetic change among the three comparison groups. So we chose the 0- and 6-h time-points to compare the gene expression data and conducted the follow-up studies. GSE27473 contains gene expression data of MCF7 parental cells and MCF7 cells silencing of the ER by shRNA. These experiments were also performed in triplicate (). GSE5840 compared the gene expression patterns of 17β-estradiol-responsive genes in wild-type MCF7 cells, tamoxifen-resistant MCF7 cells, and fulvestrant-resistant MCF7 cells (the cells were treated with 17β-estradiol for 4 h). Four replicate experiments were performed with biologically independent samples (). These three microarray datasets were generated with the use of the GPL570 Affymetrix GeneChip Human Genome U133 Plus 2.0 Array. The comparison data of each dataset were uploaded as the Supplementary Tables 1–3.
Data Processing of Differentially Expressed Genes
The up-regulated genes when ER is activated by estrogen and the down-regulated genes when ER is silenced may be genes positively regulated by ER. In contrast, the down-regulated genes when ER is activated by estrogen and the up-regulated genes when ER is silenced may be genes that are negatively regulated by ER. Three gene expression datasets (GSE11324, GSE27473, and GSE5840) containing paired gene expression data were used for the comparative analyses. DEGs in the three datasets were identified with the GEO2R online tool () using a adjust P < 0.05. Values for log2FC > 0.5 in GSE11324 and GSE5840 (comparison data for wild-type MCF7 cells) and log2FC < −0.5 in GSE27473 were considered to indicate genes positively regulated by ER, whereas log2FC < −0.5 in GSE11324 and GSE5840 (comparison data for wild-type MCF7 cells) and log2FC > 0.5 in GSE27473 were considered to indicate genes negatively regulated by ER. Next, the raw data were plotted with the Venn diagram software1, and the common DEGs among the three datasets above were obtained.
GO and Pathway Enrichment Analyses
Gene ontology analysis is commonly used to define genes and their expression products to identify unique biological properties of high-throughput transcriptome or genome data (). Kyoto Encyclopedia of Gene and Genome (KEGG) is a collection of databases dealing with genomes, diseases, biological pathways, drugs, and chemical materials (). The Database for Annotation, Visualization, and Integrated Discovery (DAVID2) is an online bioinformatics tool for functional annotation of large numbers of genes or proteins (). We used DAVID to perform GO and KEGG analyses of the list of DEGs we had identified, in terms of molecular function (MF), cellular component (CC), and biological process (BP); the terms with a P < 0.05 would be taken into consideration.
PPI Network and Modular Analyses
We used the online tool, STRING (Search Tool for the Retrieval of Interacting Genes3) (), to extract PPI information from the list of identified DEGs. Cytoscape () was used to visualize these DEGs with STRING, and the PPI network was analyzed with the MCODE (Molecular Complex Detection) plug-in to identify the significant gene clusters (degree cutoff = 2, maximum depth = 100, k-core = 2, and node score cutoff = 0.2). MCODE is a graph theoretic clustering algorithm that detects densely connected regions in PPI networks that may represent molecular complexes, and larger, more dense complexes have higher MCODE scores ().
Relapse-Free Survival Analyses
We used the Kaplan–Meier plotter online tool4 to explore the RFS of patients with ER-positive/HER2-negative breast cancer who had undergone endocrine therapies only. It was a specific group of patients who were not interfered with by other therapeutic strategies, and the recurrence of breast cancer was more likely due to endocrine therapy resistance. The Kaplan–Meier plotter for breast cancer is commonly used to assess the effect of genes on survival; the background database was established using gene expression data and survival information of breast cancer patients obtained from the GEO database (). The median of each gene expression probe as a cutoff value splitting the patient into high- and low-expression groups, and when a logrank P < 0.05 was considered to have the RFS difference. In this way, the candidate genes related to endocrine therapy resistance can be screened.
Expression Patterns Analyses of the RFS-Related Genes
We compared the gene expression data of the RFS-related genes in GSE5840 dataset, which contains the gene expression information of wild-type, tamoxifen-resistant, and fulvestrant-resistant MCF7 cells treated with estrogen. Differences among groups were analyzed by the Student’s t-test, one-way analysis of variance and least significant difference, and Student–Newman–Keuls post hoc test. P < 0.05 was considered to be statistically significant.
Workflow
The diagram of the workflow and datasets is shown in Figure 1.
FIGURE 1
Results
Identification of Genes Regulated by ER in MCF7 Cells
Three datasets containing paired gene expression data, namely, GSE11324, GSE27473, and GSE5840, were analyzed in the present study. The up-regulated genes when ER is activated by estrogen and the down-regulated genes when ER is silenced may be genes positively regulated by ER. In contrast, the down-regulated genes when ER is activated by estrogen and the up-regulated genes when ER is silenced may be genes that are negatively regulated by ER. Using the GEO2R online tool, we identified 1,833, 5,209, and 2,163 genes meeting the standard to be positively regulated by ER in GSE11324, GSE27473, and GSE5840, respectively. Meanwhile, there were 993, 4,887, and 1,845 genes to be negatively regulated by ER in GSE11324, GSE27473, and GSE5840, respectively. Venn diagrams were used to obtain the common DEGs among the three datasets. There were 230 common DEGs associated with ER activation or silence, comprising 160 genes positively regulated by ER and 70 genes negatively regulated by ER (Table 1 and Figure 2).
TABLE 1
| ER regulating genes | Genes name |
| Positively regulated by ER | PNPT1, MREG, MAX, MED13L, CAMTA1, THADA, DDX18, CENPU, MAK16, TPBG, SLC9A3R1, FKBP4, LONRF2, KRT15, MAD2L1, IFITM10, SYTL5, JAK2, C6orf141, ELF1, MAST4, CA12, ZC3H14, ADCY1, SLC25A24, SLC22A5, BRI3BP, LOC441155///ZC3H11A, NKAIN1, SIAH2, CALM3///CALM2///CALM1, POLR1B, CHML, TMPRSS3, ZBTB21, RPRD1A, EIF3J, THBS1, CALCR, SLC26A2, PGR, TMEM64, ZDHHC21, ZMYM2, PIK3R1, PTRH2, SKIL, PLA2G12A, RHOBTB1, GEMIN5, AFF3, MYB, MPPED2, NUFIP2, RPP40, DGKH, KRT13, DCLRE1B, PUS1, RET, WWC1, SMG1P5///BOLA2///SMG1P2, KIF21A, NOP16, LRPPRC, IL17RB, EGR3, SLC7A2, NIFK, THRIL///BRI3BP, USP31, SNRPA1, DSCAM, AKAP1, EIF5, NBPF4, TBC1D30, HSPA9, GTPBP4, SLC1A1, PUS7, LOC101928589///TMEM164, CXCL12, PPP4R3A, MLLT10, BAZ1A, PPAT, C1orf226, LOC101060386///BOLA2///SMG1, PTP4A1, SNX24, STC2, SGO2, ADCY9, MDM4, LOC101060275///NPIPA5///NPIPB5///NPIPB11///LOC613037///NPIPB4///NPIPB3, TTC3P1///TTC3, CAND1, CHRNA5, MTRR, MSMB, ZNF75A, ABHD2, PKIB, AMMECR1, GABPB1-AS1, BARD1, CNOT7, MIR4657///PURB, TPD52L1, SLC19A1, GNA13, GLCCI1, ZFX, NUP35, HDDC2, DMXL2, PIP4K2A, ELOVL2, MED13, TRNT1, GSPT1, DGKE, AGPAT5, MRPS25, C5orf22, RSL1D1, SLC16A1, SMIM13, CYP1B1, SLC19A2, AGR3, METTL8, RFC3, GREB1, MBOAT1, NRIP1, RAPGEFL1, FAM208B, ANAPC7, NBPF1, SOX3, TPM1, TFRC, KPNA1, RLN2, RARA, RASGRP1, MGA, TET2, SRSF1, PCDH10, PA2G4, RAB30, DEPTOR, PDZK1, IGFBP4, CHAC2, ONECUT2, SCAF11 |
| Negatively regulated by ER | HIC2, TMEM42, SPTBN1, CYP1A1, AHNAK, SOX9, PLEKHM3, JRK, TPGS1, SOCS3, LOC100288911, LOC100996760, PBXIP1, AJUBA, SNORA11E///SNORA11D///MAGED4///MAGED4B, MIR4746///UBXN6, KLHL5, IL1R1, ABTB1, CTTN, RGS3, CHEK1, FZD2, PARVA, TP53INP2, KIAA1217, CAPG, PCGF2, C10orf10, CTGF, GLMP, DDIT4, LYN, BRI3, SERPINB1, ING4, PANX2, CDKN1C, HEMK1, TBC1D2, SHISA2, ENC1, LMNA, FBXO32, SH3BP4, HS3ST1, GATA6, ORAI3, APOBEC3B, PXN, NEDD4L, BCAR3, LOC102725292///DTX2, EXT1, NGEF, KRT7, ALDH1A3, UBALD2, C11orf96, SLC17A5, INO80C, TLDC1, TAGLN, IFIT1, PLEKHM1, PARP10, BCL2L1, MIR4800///MXD4, NUMA1, GPR39 |
List of the 230 common DEGs in response to ER activation or silence identified from three datasets.
FIGURE 2
GO and KEGG Pathway Analyses for ER Regulating Genes
All 230 common DEGs were analyzed with the DAVID functional annotation tool, which included the functional categories of BP, CC, MF, and KEGG pathways. For BP, the genes positively regulated by ER mainly showed enrichment for the negative regulation of apoptotic process, negative regulation of cell migration, and adenylate cyclase-activating G-protein-coupled receptor signaling pathway (P < 0.05). The genes negatively regulated by ER mainly showed enrichment for the negative regulation of transcription, tissue homeostasis, and negative regulation of intracellular signal transduction (P < 0.05). For CC, the ER positively regulating genes were mainly enriched for the nucleoplasm, nucleolus, and membrane (P < 0.05). Moreover, the ER negatively regulating genes were mainly enriched for the cytoplasm, cytosol, and cell–cell adherens junction (P < 0.05). For MF, the ER positively regulating genes were mainly enriched for RNA binding and transcription coactivator activity (P < 0.05). The ER negatively regulating genes were mainly enriched for protein binding, actin binding, and cadherin binding involved in cell–cell adhesion (P < 0.05) (Figure 3).
FIGURE 3
Kyoto Encyclopedia of Gene and Genome pathway enrichment analysis demonstrated that ER might regulate genes along pathways associated with cancer, progesterone-mediated oocyte maturation, RNA transport, glycerophospholipid and glycerolipid metabolism, oocyte meiosis, platelet activation, phosphatidylinositol signaling system, estrogen signaling pathway, and nuclear factor κB (NF-κB) signaling pathway (P < 0.05). The genes involved in each pathway are listed in Table 2.
TABLE 2
| Term | Count | P-value | Genes |
| hsa04914: progesterone-mediated oocyte maturation | 6 | 7.80E-04 | PGR, ADCY1, MAD2L1, ADCY9, ANAPC7, PIK3R1 |
| hsa00564: glycerophospholipid metabolism | 5 | 0.008 | AGPAT5, DGKE, PLA2G12A, MBOAT1, DGKH |
| hsa00561: glycerolipid metabolism | 4 | 0.013 | AGPAT5, DGKE, MBOAT1, DGKH |
| hsa04114: oocyte meiosis | 5 | 0.014 | PGR, ADCY1, MAD2L1, ADCY9, ANAPC7 |
| hsa03013: RNA transport | 6 | 0.015 | TRNT1, EIF5, EIF3J, NUP35, RPP40, GEMIN5 |
| hsa05200: pathways in cancer | 9 | 0.016 | GNA13, MAX, ADCY1, RET, ADCY9, RASGRP1, RARA, CXCL12, PIK3R1 |
| hsa04611: platelet activation | 5 | 0.024 | GNA13, ADCY1, ADCY9, RASGRP1, PIK3R1 |
| hsa04064: NF-κB signaling pathway | 3 | 0.045 | IL1R1, LYN, BCL2L1 |
| hsa04070: phosphatidylinositol signaling system | 4 | 0.049 | DGKE, DGKH, PIP4K2A, PIK3R1 |
| hsa04915: estrogen signaling pathway | 4 | 0.049 | ADCY1, ADCY9, FKBP4, PIK3R1 |
KEGG pathway analysis of the common DEGs among the three datasets.
PPI Network and Modular Analyses
The common DEGs identified among the three datasets were imported into STRING online tool to analyze the PPI network. A total of 89 of the 230 DEGs were contained into the PPI network complex, including 89 nodes and 149 edges (Figure 4A). Moreover, we applied MCODE, a plug-in of the Cytoscape software for further analysis. Three major clusters were identified according to the MCODE score. The cluster with the top MCODE score (score = 8.0) contained 8 nodes and 28 edges: PUS7, POLR1B, MAK16, RSL1D1, NIFK, PA2G4, GTPBP4, DDX18 (Figure 4B). Followed by the cluster (MCODE score = 6.0), which contained 6 nodes and 15 edges: NEDD4L, SOCS3, SIAH2, KLHL5, ANAPC7, and,FBXO32 (Figure 4C). The cluster with the third-highest MCODE score (score = 5.0) included 5 nodes and 10 edges: ADCY9, ADCY1, GPR39, CALCR, and,RLN2 (Figure 4D).
FIGURE 4
Relapse-Free Survival Analyses
Genes involved in the KEGG pathways or major PPI clusters were imported into the Kaplan–Meier plotter online tool to analyze associations between the genes’ expression and RFS from patients with ER-positive/HER2-negative breast cancer who had undergone endocrine therapies only. It was a specific group of patients who were not interfered with by other therapeutic strategies, and the recurrence of breast cancer was more likely due to endocrine therapy resistance. Among these, the expression of 16 genes was found to be significantly associated with RFS (P < 0.05), whereas the other 28 genes had no significant associations (Table 3). The 16 genes were ADCY9, ANAPC7, CALCR, CXCL12, DDX18, EIF3J, FKBP4, GEMIN5, GTPBP4, MAD2L1, MAX, NUP35, POLR1B, PUS7, RSL1D1, and SOCS3 (Figure 5). In this way, the candidate genes related to endocrine therapy resistance could be screened.
TABLE 3
| Category | Genes |
| Genes significantly associated with RFS (P < 0.05) | ADCY9, ANAPC7, CALCR, CXCL12, DDX18, EIF3J, FKBP4, GEMIN5, GTPBP4, MAD2L1, MAX, NUP35, POLR1B, PUS7, RSL1D1, SOCS3 |
| No significant associations with RFS (P > 0.05) | ADCY1, AGPAT5, BCL2L1, DGKE, DGKH, EIF5, FBXO32, GNA13, GPR39, IL1R1, KLHL5, LYN, MAK16, MBOAT1, NEDD4L, NIFK, PA2G4, PGR, PIK3R1, PIP4K2A, PLA2G12A, RARA, RASGRP1, RET, RLN2, RPP40, SIAH2, TRNT1 |
The prognostic value (RFS) of the selected genes in patients with ER-positive/HER2-negative breast cancer who had undergone endocrine therapies only.
These genes were involved in the KEGG pathways or major PPI clusters. RFS, relapse-free survival.
FIGURE 5
Expression Patterns Analyses of the Screened Genes Among Wild-Type, Tamoxifen-Resistant, and Fulvestrant-Resistant MCF7 Cells
We analyzed the expression data of the 16 RFS-related genes in GSE5840 dataset, which contains comparison data, including wild-type MCF7 cells treated with 17β-estradiol versus negative control [dimethyl sulfoxide (DMSO)], tamoxifen-resistant MCF7 cells treated with 17β-estradiol versus negative control (DMSO), and fulvestrant-resistant MCF7 cells treated with 17β-estradiol versus negative control (DMSO). First, we compared the gene expression data among wild-type, tamoxifen-resistant, and fulvestrant-resistant MCF7 cells in the negative control groups. There were different gene expression patterns among the three cell lines. For instance, the genes such as ANAPC7 and DDX18 were significantly up-regulated (log2FC > 0.5 and P < 0.05) in the fulvestrant-resistant cells, and the MAD2L1, RSL1D1, and CALCR were significantly up-regulated (log2FC > 0.5 and P < 0.05) in the tamoxifen-resistant cells; all the comparison results are shown in Figure 6. Then, we analyzed the gene expression data among the three kinds of cell lines treated with estrogen or DMSO (negative control) (Figure 7). Comparing with wild-type cells, nearly all of the genes’ expression levels changed less in fulvestrant-resistant cells when treated with estrogen. However, in tamoxifen-resistant cells, the change patterns of gene expression under the estrogen treatment were similar to that of wild-type cells except for DDX18, CXCL12, FKBP4, NUP35, MAD2L1, RSL1D1 (changed less), and CALCR (changed more).
FIGURE 6
FIGURE 7
Discussion
Estrogen and ER play essential roles in the development and progression of ER-positive breast cancer. At the same time, the phenomenon of endocrine therapy resistance occurs frequently and complicates patient management (). However, the regulatory network of the ER has not yet been fully elucidated (). The present study used bioinformatics methods to interrogate three gene expression datasets from MCF7 cells treated with estrogen, MCF7 cells subjected to the silencing of the ER, and tamoxifen- and fulvestrant-resistant MCF7 cells treated with 17β-estradiol, respectively. By using these datasets, the function of the ER could be explored from different perspectives. Not only did we perform GO and KEGG pathway analyses, but also screened RFS-related genes and explored their expression pattern changes in tamoxifen- and fulvestrant-resistant cells.
There are two classes of ER, including nuclear ER (intracellular receptor) and membrane ER (mostly G-protein-coupled receptor) (). For nuclear ER (ERα and ERβ), once activated by estrogen, the ER can translocate into the nucleus and regulate the activity of different genes (). In the present study, GO and KEGG analyses showed that ER positively regulating genes were mainly localized in the nucleoplasm, nucleolus, and membrane, involving biological processes such as negative regulation of apoptotic process and adenylate cyclase-activating G-protein-coupled receptors signaling pathway, which were corresponding to the two forms of ER as the nuclear receptor and the membrane receptor. The molecular functions such as RNA binding and transcription coactivator activity could create conditions for genetic transcription and cell proliferation. Meanwhile, the biological processes for ER negatively regulating genes were negative regulation of transcription and negative regulation of intracellular signal transduction, further reflecting the main functions of the two forms of ER.
DEGs involved in the KEGG pathways or major PPI clusters were analyzed with Kaplan–Meier plotter to determine whether there were any associations between these genes and the RFS of patients with ER-positive/HER2-negative breast cancer who had undergone endocrine therapies only. It is a specific group of patients who were not interfered with by other therapeutic strategies; such analysis may best reflect the interplay between ER-regulated genes and the development of endocrine resistance. ADCY9, ANAPC7, CALCR, CXCL12, DDX18, EIF3J, FKBP4, GEMIN5, GTPBP4, MAD2L1, MAX, NUP35, POLR1B, PUS7, RSL1D1, and SOCS3 were found to be significantly associated with RFS. We subsequently analyzed the expression patterns of these genes among wild-type, tamoxifen-resistant, and fulvestrant-resistant MCF7 cells. Interestingly enough, nearly all of the genes’ expression level changed less when treated with estrogen in fulvestrant-resistant cells compared with the wild type, indicating that the function of ER is weakened or disappeared in fulvestrant-resistant MCF7 cells. Fulvestrant is a selective ER degrader; it works by binding to and destabilizing the ERs (). If ERs ceased to function in breast cancer cells, fulvestrant resistance would be occurring. ANAPC7 and DDX18 were ER positively regulating genes in wild-type MCF7 cells. However, they were significantly up-regulated in the fulvestrant-resistant cells without estrogen treatment, speculating that these genes may play essential roles in promoting the survival of fulvestrant-resistant breast cancer cells.
When tamoxifen-resistant cells were treated with estrogen, the expression change patterns of many genes were similar to that of wild-type cells, suggesting that the function of ER exists in tamoxifen-resistant MCF7 cells. As the selective ER modulator, tamoxifen blocks the effects of estrogen by attaching to the ERs in breast cells (). The existence of functional ERs explains why fulvestrant treatment may still be effective when tamoxifen resistance occurs. Meanwhile, in the negative control groups, MAD2L1, RSL1D1, and CALCR were found to be significantly up-regulated in the tamoxifen-resistant cells compared with wild-type cells. It was speculated that these genes might play crucial roles in tamoxifen-resistant breast cancer cells.
Dead-box RNA helicase 18 (DDX18) is an essential factor in cell cycle progression in zebrafish hematopoietic cells and is mutated in some patients with acute myeloid leukemia (). identified DDX18 as a novel driver of endocrine resistance in breast cancer. They found a significant inhibition of the proliferation of endocrine-resistant cell lines based on an increased G1 phase cell population when DDX18 was silenced. By analyzing clinical samples, they discovered that the mRNA levels of DDX18 were significantly correlated with poor prognosis in breast cancer patients (). It is consistent with our results. High expression of DDX18 was correlated with significantly worse RFS in ER-positive/HER2-negative breast cancer patients who had undergone endocrine therapies only, and the expression level of DDX18 was continuously high in the fulvestrant-resistant breast cancer cells. investigated six microarray datasets from the Gene Expression Omnibus consisting of 405 triple-negative breast cancer (ER/PR/HER2-negative, TNBC) and 463 non-TNBC samples and identified 1,075 DEGs; DDX18 was 1 of 12 up-regulated genes identified as essential by a machine learning-based feature selection method (). In the present study, we found that the function of the ER is weakened or disappeared in fulvestrant-resistant MCF7 cells, suggesting that these cells might possess the characteristics of TNBC. With the DDX18 exhibits significantly higher expression in TNBC and is no longer regulated by ER in fulvestrant-resistant MCF7 cells. These phenomena suggest that DDX18 might play a vital role in endocrine therapy resistance and that its activity might confer characteristics of TNBC onto ER-positive breast cancer.
ANAPC7 gene encodes a tetratricopeptide repeat-containing component of the anaphase-promoting complex/cyclosome (APC/C), a multisubunit ubiquitin ligase that is essential for mitosis by targeting a number of cell cycle regulators such as cyclin B1 and promoting timely degradation (). Reducing the activity of APC/C delays mitotic progression, whereas APC/C activity loss is lethal (; ). In our study, high expression of ANAPC7 was correlated with significantly worse RFS in the specific cohort of breast cancer patients, and the expression level of ANAPC7 was continuously high in the fulvestrant-resistant breast cancer cells. Considering the vital function of APC/C in eukaryotic cells mitosis (), ANAPC7 might play an essential role in fulvestrant resistance, and ANAPC7 (anaphase-promoting complex subunit 7) was expected to be a novel biomarker or therapeutic target for fulvestrant-resistant breast cancer.
Other up-regulated genes such as GEMIN5 and POLR1B (values of log2FC were between 0 and 0.5, and P < 0.05) might also be associated with fulvestrant resistance. Gem nuclear organelle associated protein 5 (Gemin5) is an RNA-binding protein that was identified as a component of the survival of motor neurons (SMN) complex (). The SMN complex plays a critical role in mRNA splicing and may mediate the assembly and transport of other classes of ribonucleoproteins (). GEMIN5 codes Gemin5; alteration of GEMIN5 expression may play a role in alternative mRNA splicing and tumor cell motility (). Recently, a study reported that POLR1B is up-regulated in non-small cell lung cancer and may serve an important modulator of lung cancer cell proliferation (). Our research added to the understanding of GEMIN5 and POLR1B genes, and they might serve as biomarkers for fulvestrant resistance in ER-positive breast cancer.
In tamoxifen-resistant cells, MAD2L1, RSL1D1, and CALCR were found to be significantly up-regulated compared with the wild-type cells. MAD2L1 (mitotic arrest deficient 2 like 1) plays an essential role in supervising chromosome segregation as the component of spindle checkpoint during mitosis (). Studies showed that MAD2L1 presented overexpression in breast cancer and was significantly associated with higher clinical stage, higher histological grade, aggressive tumors, and worse disease-free survival (; ). RSL1D1 (ribosomal L1 domain containing 1) is a nucleolar protein that has been demonstrated to delay cellular senescence and serve as an independent prognostic factor in prostate cancer (). CALCR (coding calcitonin receptor) has similar gene expression patterns to MAD2L1 and RSL1D1 in tamoxifen-resistant MCF7 cells; the expression of these genes was all correlated with worse RFS in ER-positive/HER2-negative breast cancer patients who had undergone endocrine therapies only. In addition, compared with the wild-type MCF7 cells, the expression level of MAD2L1, RSL1D1, and CALCR was significantly up-regulated in tamoxifen-resistant cells even without estrogen treatment. However, there has not been sufficient research on the relationship between the proteins coded by these genes and endocrine therapy resistance in breast cancer. Therefore, MAD2L1, RSL1D1, and CALCR might be promising candidates for further research in endocrine therapy-resistant breast cancer as potential therapeutic targets or prognostic markers.
By identifying the significant genes interacting with the ER and their expression pattern changing when tamoxifen or fulvestrant resistance occurs, we not only can obtain a better understanding of this essential receptor’s regulatory network but also can speculate some possible mechanisms of endocrine resistance. The main assumptions are as follows: (1) in fulvestrant-resistant breast cancer cells, the expression change of most estrogen-regulated genes was not evident under estrogen treatment, indicating that the function of ER is weakened or disappeared in fulvestrant-resistant cells. Fulvestrant works by destabilizing the ERs, and if ERs ceased to function in breast cancer cells, fulvestrant resistance would occur. (2) DDX18 and ANAPC7 might play a vital role in fulvestrant-resistant breast cancer cells because of their crucial functions in cell cycle progression and eukaryotic cells mitosis (; ) and continuous overexpression in the fulvestrant-resistant breast cancer cells. So they were expected to be novel biomarkers or therapeutic targets for further study. (3) In tamoxifen-resistant breast cancer cells, the expression changes of many genes were similar to wild-type cells under estrogen treatment, suggesting that the ER’s function exists in tamoxifen-resistant cells, and explaining that fulvestrant treatment might still be valid when tamoxifen resistance occurs. (4) Some genes such as MAD2L1, RSL1D1, and CALCR were found to be significantly up-regulated in the tamoxifen-resistant cells, indicating that tamoxifen might not completely block the estrogen signaling in tamoxifen-resistant breast cancer cells. These highly expressed genes would be potential biomarkers for tamoxifen resistance, and suppression of them might be an important direction to overcome tamoxifen resistance in the future.
In summary, this bioinformatics analysis study identified the significant genes regulated by ER in ER-positive breast cancer cells. As endocrine resistance in ER-positive breast cancer is likely to be attributable to the abnormal regulation of ER network, the expression pattern changes of these genes were explored in tamoxifen- and fulvestrant-resistant breast cancer cells. Although these predictions need to be further validated, the present study provided useful insights regarding potential biomarkers and the pathomechanisms of ER-positive breast cancer resistant to endocrine therapy.
Statements
Data availability statement
The datasets generated for this study can be found in the GEO database, GSE11324, GSE27473, GSE5840.
Author contributions
JW, YF, and ZW conceived the idea and revised the final manuscript. RC and LQ designed the work, analyzed the data, and drafted the manuscript. XK interpreted the data and revised the manuscript. All authors read and approved the final manuscript. They agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work were appropriately investigated and resolved.
Funding
This study was supported by the National Natural Science Foundation of China (No. 81872160). The Ph.D. Innovation Funding of Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College (No. C2019- 1051-09).
Acknowledgments
We would like to acknowledge the contributors of GSE11324, GSE27473, and GSE5840. Furthermore, we would like to thank the NCBI-GEO database, Venn diagrams, DAVID, STRING, and Kaplan–Meier plotter online tools.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.538734/full#supplementary-material
SUPPLEMENTARY TABLE 1The comparison results from GSE11324 dataset (MCF7 cells were stimulated with estrogen for 6 vs. 0 h).
SUPPLEMENTARY TABLE 2The comparison results from GSE27473 dataset (MCF7 parental cells vs. MCF7 cells silencing of ER).
SUPPLEMENTARY TABLE 3The comparison results from GSE5840 dataset (wild-type MCF7 cells were treated with 17β-estradiol for 4 h vs. DMSO for 4 h).
References
1
Al SalehS.Al MullaF.LuqmaniY. A. (2011). Estrogen receptor silencing induces epithelial to mesenchymal transition in human breast cancer cells.PLoS One6:e20610. 10.1371/journal.pone.0020610
2
AshburnerM.BallC. A.BlakeJ. A.BotsteinD.ButlerH.CherryJ. M.et al (2000). Gene ontology: tool for the unification of biology. The gene ontology consortium.Nat. Genet.2525–29.
3
BaderG. D.HogueC. W. (2003). An automated method for finding molecular complexes in large protein interaction networks.BMC Bioinform.4:2. 10.1186/1471-2105-4-2
4
BhuvaD. D.CursonsJ.SmythG. K.DavisM. J. (2019). Differential co-expression-based detection of conditional relationships in transcriptional data: comparative analysis and application to breast cancer.Genome Biol.20:236.
5
CarrollJ. S.MeyerC. A.SongJ.LiW.GeistlingerT. R.EeckhouteJ.et al (2006). Genome-wide analysis of estrogen receptor binding sites.Nat. Genet.381289–1297.
6
ClarkeR.TysonJ. J.DixonJ. M. (2015). Endocrine resistance in breast cancer–an overview and update.Mol. Cell. Endocrinol.418(Pt 3), 220–234. 10.1016/j.mce.2015.09.035
7
DavisS.MeltzerP. S. (2007). GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor.Bioinformatics231846–1847. 10.1093/bioinformatics/btm254
8
DeSantisC. E.MaJ.GaudetM. M.NewmanL. A.MillerK. D.Goding SauerA.et al (2019). Breast cancer statistics, 2019.CA Cancer J. Clin.69438–451.
9
DuJ.YuanZ.MaZ.SongJ.XieX.ChenY. (2014). KEGG-PATH: Kyoto encyclopedia of genes and genomes-based pathway analysis using a path analysis model.Mol. Biosyst.102441–2447. 10.1039/c4mb00287c
10
Early Breast, Cancer Trialists, Collaborative GroupDaviesC.GodwinJ.GrayR.et al (2011). Relevance of breast cancer hormone receptors and other factors to the efficacy of adjuvant tamoxifen: patient-level meta-analysis of randomised trials.Lancet378771–784. 10.1016/s0140-6736(11)60993-8
11
FanM.YanP. S.Hartman-FreyC.ChenL.PaikH.OyerS. L.et al (2006). Diverse gene expression and DNA methylation profiles correlate with differential adaptation of breast cancer cells to the antiestrogens tamoxifen and fulvestrant.Cancer Res.6611954–11966. 10.1158/0008-5472.can-06-1666
12
FengH.GuZ. Y.LiQ.LiuQ. H.YangX. Y.ZhangJ. J. (2019). Identification of significant genes with poor prognosis in ovarian cancer via bioinformatical analysis.J. Ovarian Res.12:35.
13
GoodsellD. S. (2002). The molecular perspective: tamoxifen and the estrogen receptor.Stem Cells20267–268. 10.1634/stemcells.20-3-267
14
GuoY.ZhangX.YangM.MiaoX.ShiY.YaoJ.et al (2010). Functional evaluation of missense variations in the human MAD1L1 and MAD2L1 genes and their impact on susceptibility to lung cancer.J. Med. Genet.47616–622. 10.1136/jmg.2009.074252
15
GyorffyB.LanczkyA.EklundA. C.DenkertC.BudcziesJ.LiQ.et al (2010). An online survival analysis tool to rapidly assess the effect of 22,277 genes on breast cancer prognosis using microarray data of 1,809 patients.Breast Cancer Res. Treat.123725–731. 10.1007/s10549-009-0674-9
16
Huang, daW.ShermanB. T.LempickiR. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.Nat. Protoc.444–57. 10.1038/nprot.2008.211
17
JelovacD.MacedoL.GoloubevaO. G.HandrattaV.BrodieA. M. (2005). Additive antitumor effect of aromatase inhibitor letrozole and antiestrogen fulvestrant in a postmenopausal breast cancer model.Cancer Res.655439–5444. 10.1158/0008-5472.can-04-2782
18
JordanV. C.BrodieA. M. (2007). Development and evolution of therapies targeted to the estrogen receptor for the treatment and prevention of breast cancer.Steroids727–25. 10.1016/j.steroids.2006.10.009
19
LaiA. C.CrewsC. M. (2017). Induced protein degradation: an emerging drug discovery paradigm.Nat. Rev. Drug Discov.16101–114. 10.1038/nrd.2016.211
20
LeeJ. H.HorakC. E.KhannaC.MengZ.YuL. R.VeenstraT. D.et al (2008). Alterations in Gemin5 expression contribute to alternative mRNA splicing patterns and tumor cell motility.Cancer Res.68639–644. 10.1158/0008-5472.can-07-2632
21
LiX. P.JiaoJ. U.LuL. I.ZouQ.ZhuS.ZhangY. (2016). Overexpression of ribosomal L1 domain containing 1 is associated with an aggressive phenotype and a poor prognosis in patients with prostate cancer.Oncol. Lett.112839–2844. 10.3892/ol.2016.4294
22
MagnusonT.EpsteinC. J. (1984). Oligosyndactyly: a lethal mutation in the mouse that results in mitotic arrest very early in development.Cell38823–833. 10.1016/0092-8674(84)90277-0
23
MassenetS.PellizzoniL.PaushkinS.MattajI. W.DreyfussG. (2002). The SMN complex is associated with snRNPs throughout their cytoplasmic assembly pathway.Mol. Cell Biol.226533–6541. 10.1128/mcb.22.18.6533-6541.2002
24
NaoremL. D.MuthaiyanM.VenkatesanA. (2019). Integrated network analysis and machine learning approach for the identification of key genes of triple-negative breast cancer.J. Cell. Biochem.1206154–6167. 10.1002/jcb.27903
25
PayneE. M.BolliN.RhodesJ.Abdel-WahabO. I.LevineR.HedvatC. V.et al (2011). Ddx18 is essential for cell-cycle progression in zebrafish hematopoietic cells and is mutated in human AML.Blood118903–915. 10.1182/blood-2010-11-318022
26
PineiroD.Fernandez-ChamorroJ.Francisco-VelillaR.Martinez-SalasE. (2015). Gemin5: a multitasking RNA-binding protein involved in translation control.Biomolecules5528–544. 10.3390/biom5020528
27
PinesJ. (2011). Cubism and the cell cycle: the many faces of the APC/C.Nat. Rev. Mol. Cell Biol.12427–438. 10.1038/nrm3132
28
RaniA.StebbingJ.GiamasG.MurphyJ. (2019). Endocrine resistance in hormone receptor positive breast cancer-from mechanism to therapy.Front. Endocrinol.10:245. 10.3389/fendo.2019.00245
29
RazandiM.PedramA.GreeneG. L.LevinE. R. (1999). Cell membrane and nuclear estrogen receptors (ERs) originate from a single transcript: studies of ERalpha and ERbeta expressed in Chinese hamster ovary cells.Mol. Endocrinol.13307–319. 10.1210/me.13.2.307
30
RedmondA. M.ByrneC.BaneF. T.BrownG. D.TibbittsP.O’BrienK.et al (2015). Genomic interaction between ER and HMGB2 identifies DDX18 as a novel driver of endocrine resistance in breast cancer cells.Oncogene343871–3880. 10.1038/onc.2014.323
31
ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303
32
SzklarczykD.FranceschiniA.WyderS.ForslundK.HellerD.Huerta-CepasJ.et al (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life.Nucleic Acids Res.43D447–D452.
33
VogelsteinB.PapadopoulosN.VelculescuV. E.ZhouS.DiazL. A.Jr.et al (2013). Cancer genome landscapes.Science3391546–1558. 10.1126/science.1235122
34
WangZ.KatsarosD.ShenY.FuY.CanutoE. M.BenedettoC.et al (2015). Biological and clinical significance of MAD2L1 and BUB1, genes frequently appearing in expression signatures for breast cancer prognosis.PLoS One10:e0136246. 10.1371/journal.pone.0136246
35
WildT.BudzowskaM.HellmuthS.EibesS.KaremoreG.BarisicM.et al (2018). Deletion of APC7 or APC16 allows proliferation of human cells without the spindle assembly checkpoint.Cell Rep.252317.e8–2328.e5.
36
WirthK. G.RicciR.Gimenez-AbianJ. F.TaghybeegluS.KudoN. R.JochumW.et al (2004). Loss of the anaphase-promoting complex in quiescent cells causes unscheduled hepatocyte proliferation.Genes Dev.1888–98. 10.1101/gad.285404
37
YangF.LiuH.ZhaoJ.MaX.QiW. (2020). POLR1B is upregulated and promotes cell proliferation in non-small cell lung cancer.Oncol. Lett.19671–680.
38
YasarP.AyazG.UserS. D.GupurG.MuyanM. (2017). Molecular mechanism of estrogen-estrogen receptor signaling.Reprod. Med. Biol.164–20. 10.1002/rmb2.12006
39
ZhuX. F.YiM.HeJ.TangW.LuM. Y.LiT.et al (2017). Pathological significance of MAD2L1 in breast cancer: an immunohistochemical study and meta analysis.Int. J. Clin. Exp. Pathol.109190–9201.
Summary
Keywords
bioinformatics analysis, microarray, differentially expressed genes, ER-positive breast cancer, tamoxifen or fulvestrant resistance
Citation
Cheng R, Qi L, Kong X, Wang Z, Fang Y and Wang J (2020) Identification of the Significant Genes Regulated by Estrogen Receptor in Estrogen Receptor-Positive Breast Cancer and Their Expression Pattern Changes When Tamoxifen or Fulvestrant Resistance Occurs. Front. Genet. 11:538734. doi: 10.3389/fgene.2020.538734
Received
17 March 2020
Accepted
27 August 2020
Published
29 September 2020
Volume
11 - 2020
Edited by
Paola Lecca, Free University of Bozen-Bolzano, Italy
Reviewed by
Damu Tang, McMaster University, Canada; Chiara Damiani, University of Milano-Bicocca, Italy
Updates
Copyright
© 2020 Cheng, Qi, Kong, Wang, Fang and Wang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jing Wang, wangjing@cicams.ac.cnYi Fang, fangyi@cicams.ac.cnZhongzhao Wang, wangzhongzhao@cicams.ac.cn
†These authors have contributed equally to this work
This article was submitted to Computational Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.