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

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 (DeSantis et al., 2019). 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 (Early Breast Cancer Trialists' Collaborative Group et al., 2011). Endocrine therapies, such as tamoxifen, fulvestrant, and aromatase inhibitors, which target estrogen or ER, are useful in breast cancer treatment (Jelovac et al., 2005;Jordan and Brodie, 2007). 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 (Clarke et al., 2015). 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) (Vogelstein et al., 2013), 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 (Feng et al., 2019).
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 (ERpositive breast cancer cell line) treated with estrogen, MCF7 cells with silencing of the ER, and tamoxifen-and fulvestrantresistant 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.

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 (Carroll et al., 2006). 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 (Al Saleh et al., 2011). 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 (Fan et al., 2006). 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 (Davis and Meltzer, 2007) using a adjust P < 0.05. Values for log 2 FC > 0.5 in GSE11324 and GSE5840 (comparison data for wild-type MCF7 cells) and log 2 FC < −0.5 in GSE27473 were considered to indicate genes positively regulated by ER, whereas log 2 FC < −0.5 in GSE11324 and GSE5840 (comparison data for wild-type MCF7 cells) and log 2 FC > 0.5 in GSE27473 were considered to indicate genes negatively regulated by ER. Next, the raw data were plotted with the Venn diagram software 1 , 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 (Ashburner et al., 2000). Kyoto Encyclopedia of Gene and Genome (KEGG) is a collection of databases dealing with genomes, diseases, biological pathways, drugs, and chemical materials (Du et al., 2014). The Database for Annotation, Visualization, and Integrated Discovery (DAVID 2 ) is an online bioinformatics tool for functional annotation of large numbers of genes or proteins (Huang da et al., 2009). 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 Genes 3 ) (Szklarczyk et al., 2015), to extract PPI information from the list of identified DEGs. Cytoscape (Shannon et al., 2003) 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 FIGURE 3 | Gene ontology (GO) analysis of the common DEGs among the three breast cancer datasets. (A) DEGs positively regulated by ER were analyzed by GO enrichment, the top 3 in P-value ranking (sorting from small to large) in biological process (BP), cellular component (CC), and molecular function (MF). (B) DEGs negatively regulated by ER were analyzed by GO enrichment, the top 3 in P-value ranking (sorting from small to large) in BP, CC, and MF (P < 0.05). 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 (Bader and Hogue, 2003).

Relapse-Free Survival Analyses
We used the Kaplan-Meier plotter online tool 4 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 (Gyorffy et al., 2010). 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 RFSrelated genes in GSE5840 dataset, which contains the gene expression information of wild-type, tamoxifenresistant, 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.

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).

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-proteincoupled 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 cellcell 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). 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.

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.

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 (log 2 FC > 0.5 and P < 0.05) in the fulvestrant-resistant cells, and the MAD2L1, RSL1D1, and CALCR were significantly up-regulated (log 2 FC > 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).

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 (Rani et al., 2019). However, the regulatory network of the ER has not yet been fully elucidated (Bhuva et al., 2019). 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 fulvestrantresistant cells.
There are two classes of ER, including nuclear ER (intracellular receptor) and membrane ER (mostly G-proteincoupled receptor) (Razandi et al., 1999). For nuclear ER (ERα and ERβ), once activated by estrogen, the ER can translocate into the nucleus and regulate the activity of different genes (Yasar et al., 2017). 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 cyclaseactivating 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, tamoxifenresistant, and fulvestrant-resistant MCF7 cells. Interestingly enough, nearly all of the genes' expression level changed FIGURE 5 | Kaplan-Meier plotter was used to analyze associations between the expressions of the genes involved in the KEGG pathways or major PPI clusters and relapse-free survival (RFS) in patients with ER-positive/HER2-negative breast cancer who had undergone endocrine therapies only. This analysis identified 16 of 44 genes as being associated with a significantly different RFS (P < 0.05).
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 (Lai and Crews, 2017). 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, FIGURE 6 | The comparisons of the screened RFS-related genes' expression among the wild-type, tamoxifen-resistant, and fulvestrant-resistant MCF7 cells (treated cells with DMSO) in GSE5840 dataset (*P < 0.05, **P < 0.01, ***P < 0.001).
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 fulvestrantresistant 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 (Goodsell, 2002). 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 tamoxifenresistant 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 (Payne et al., 2011). Redmond et al. (2015) 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 (Redmond et al., 2015). 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. Naorem et al. (2019) 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 (Naorem et al., 2019). 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 (Wild et al., 2018). Reducing the activity of APC/C delays mitotic progression, whereas APC/C activity loss is lethal (Magnuson and Epstein, 1984;Wirth et al., 2004). 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 (Pines, 2011), 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 log 2 FC 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 RNAbinding protein that was identified as a component of the survival of motor neurons (SMN) complex (Pineiro et al., 2015). The SMN complex plays a critical role in mRNA splicing and may mediate the assembly and transport of other classes of ribonucleoproteins (Massenet et al., 2002). GEMIN5 codes Gemin5; alteration of GEMIN5 expression may play a role in alternative mRNA splicing and tumor cell motility (Lee et al., 2008). 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 (Yang et al., 2020). Our research added to the understanding of GEMIN5 and POLR1B genes, and they might serve as biomarkers for fulvestrant resistance in ERpositive 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 (Guo et al., 2010). 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 (Wang et al., 2015;Zhu et al., 2017). 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 (Li et al., 2016). CALCR (coding calcitonin receptor) has similar gene expression patterns to MAD2L1 and RSL1D1 in tamoxifenresistant 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 upregulated 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 (Pines, 2011;Redmond et al., 2015) and continuous overexpression in the fulvestrantresistant 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 tamoxifenresistant 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 fulvestrantresistant 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.

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.