Identification and Analysis of Estrogen Receptor α Promoting Tamoxifen Resistance-Related lncRNAs

70-75% breast cancer patients are estrogen receptor alpha positive (ERα+), and the antiestrogen drug tamoxifen has been used for the past three decades. However, in 20-30% of these patients, tamoxifen therapy fails due to intrinsic or acquired resistance. A previous study has showed ERα signaling still exerts significant roles in the development of tamoxifen resistance and several lncRNAs have been demonstrated important roles in tamoxifen resistance. But ERα directly regulated and tamoxifen resistance related lncRNAs remain to be discovered. We reanalyze the published ERα chromatin immunoprecipitation-seq (ChIP-seq) and RNA-seq data of tamoxifen-sensitive (MCF-7/WT) and tamoxifen-resistant (MCF-7/TamR) breast cancer cells. We demonstrate that there are differential ERα recruitment events and the differentials may alert the expression profile in MCF-7/WT and MCF-7/TamR cells. Furthermore, we make an overlap of the ERα binding lncRNAs and differentially expressed lncRNAs and get 49 ERα positively regulated lncRNAs. Among these lncRNAs, the expression levels of AC117383.1, AC144450.1, RP11-15H20.6, and ATXN1-AS1 are negatively correlated with the survival probability of breast cancer patients and ELOVL2-AS1, PCOLCE-AS1, ITGA9-AS1, and FLNB-AS1 are positively correlated. These lncRNAs may be potential diagnosis or prognosis markers of tamoxifen resistance.


Introduction
Breast cancer is the most commonly diagnosed cancer and the leading cause of cancer death among female cancers, with an estimated 2.1 million newly diagnosed female breast cancer cases in 2018 around the world [1]. About 75% of all breast cancer cases are estrogen receptor alpha positive (ERα+); thus, the selective ER modulator tamoxifen has been the standard treatment for ERα+ breast cancer patients. However, up to 30% of these cases relapse due to tamoxifen resistance [2]. Various mechanisms have been proposed to explain tamoxifen resistance, including the abnormal ERα signaling pathway, alterations in cell cycle and apoptosis molecules, and activation of TGF-β or NF-кB signaling pathway [3,4]. Noncoding RNAs (ncRNAs) also play an important role in tamoxifen resistance [5]. Long noncoding RNAs (lncRNAs) are a new type ncRNAs with more than 200 nucleotides, and they are involved in a wide variety of physiological and pathological processes [6]. What is more, lncRNAs have been demonstrated to play a significant role in drug resistance. Urothelial carcinoma associated 1 (UCA1) has been confirmed to contribute to multiple cancer drug resistance, including cisplatin resistance in bladder cancer and ovary cancer, gefitinib resistance in lung cancer, and tamoxifen resistance in breast cancer [7]. Besides, breast cancer antiestrogen resistance 4 (BCAR4), lncRNA-ROR (ROR, regulator of reprogramming), colon cancer associated transcript 2 (CCAT2), DSCAM-AS1, and LINC00894 are also reported to enhance or attenuate tamoxifen resistance [8]. But the important role and complex working mechanism of lncRNAs in tamoxifen resistance still need to be illustrated, and there is also an urgent need to discover more novel diagnosis and prognosis markers.
In this work, we reanalyze the published ERα ChIP-seq and RNA-seq data of tamoxifen-sensitive (MCF-7/WT) and tamoxifen-resistant (MCF-7/TamR) breast cancer cells [9]. We find that there are differential ERα recruitment events, and these events may lead to changes of expression profile between MCF-7/WT and MCF-7/TamR cells. And we also find that the ERα signaling is also estradiol (E2) dependent in MCF-7/TamR cells. Here, we provide insight into how ERα signaling promotes tamoxifen resistance through lncRNAs and get eight lncRNAs that may act as potential diagnosis or prognosis markers.

Materials and Methods
2.1. ChIP-seq and RNA-seq Data Analysis. The ChIP-seq and RNA-seq data were downloaded from the GEO dataset (GSE86538). Quality control of the raw reads was performed by FastQC [10]. ChIP-seq reads were aligned to the hg19 genome assembly using Bowtie 2 [11]. ChIP-seq peaks were called using MACS, and the p value 1.00e-05 was used as the cut-off [12]. RNA-seq reads were aligned and annotated using hg19 with Hisat2, and SAMtools was used to sort and reorder the SAM files [13,14]. HT-seq and DESeq2 were used to analyze the differential expression, and fold change ≥ 2 or ≤0.5 and p value ≤ 0.05 were used as the cut-off [15,16].

2.2.
Visualization of the ChIP-seq Peaks and RNA-seq Differential Expression Genes. The aligned ChIP-seq reads were annotated and visualized using ChIPseeker, and (-3 kb_3 kb) was set as the transcription start site (TSS) region [17]. The RNA-seq differential expression genes (fold change ≥ 2 or ≤0.5, p value ≤ 0.05) were visualized as a heatmap using R package of ggplot2 [18].

Functional Enrichment
Analysis. The ERα ChIP-seq binding genes in the (-3 kb_3 kb) region and upregulated genes in MCF-7/TamR cells (fold change ≥ 2, p value ≤ 0.05) were selected for the further functional enrichment analysis. We did the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the candidate genes using R package of ClusterProfiler, and p value ≤ 0.05 was used as the cut-off [19].

The Correlation Analysis and Survival Analysis.
We analyzed the correlationship between expression level of candidate lncRNAs and ERα status using the breast invasive carcinoma (BRCA) data from The Cancer Genome Atlas (TCGA). We also analyzed whether the expression levels of candidate lncRNAs correlated with patient survival time (based on p value from the univariate Cox proportional hazards model and log-rank test) and visualized their correlationship through Kaplan-Meier plot. Spearman's rank correlation analysis was used to detect significant correlations between lncRNA and mRNA with a coefficient (absolute value) cut-off of 0.3 and a p value cut-off of 0.01. The correlation analysis and survival analysis were performed through the TANRIC website [20].

BioMed Research International
These functions were consistent with the morphological changes and indicative of a metastatic phenotype that may promote tamoxifen resistance.   Figure 1: Proposed work flow of the procedure to identify the ERα promoting tamoxifen resistance-related lncRNAs. In breast cancer tamoxifen-sensitive cell MCF-7/WT (WT) and tamoxifen-resistant cell MCF-7/TamR (TAMR), (a) estrogen receptor α (ERα) binding or no binding events are determined by ChIP-seq and the differential ERα redistribution is visualized by UCSC genome browser; (b) the differential gene expression profile is investigated by RNA-seq and presented in a heatmap. (c) By analyzing the ChIP-seq and RNA-seq data from GSE86538, the overlap genes (candidate genes) may be ERα promoted and tamoxifen resistance related. And further, these genes will be verified by functional enrichment analysis by clusterProfiler package and survival analysis via the TANRIC website.  (Figures 6(a)-6(h)). To our knowledge, these eight lncRNAs have not been previously studied, and they may be promising prognosis markers of tamoxifen resistance.  3.5. Correlation Analysis between Candidate lncRNAs and mRNAs. By Spearman's rank correlation analysis between lncRNA and mRNA, we identified tens of significantly corre-lated mRNAs of the eight candidate lncRNAs (Table S3)

BioMed Research International
ITGA9-AS1, and FLNB-AS1). A previous study demonstrated that antisense lncRNA can hybridize with sense RNA and form RNA duplexes, which promoted the stability of the sense RNA [25]. In our results, we found that ATXN1-AS1, ELOVL2-AS1, and FLNB-AS1 were separately positively correlated with their sense RNA ATXN1, ELOVL2, and FLNB, which were essential genes in cellular activities. What is more, the rest five lncRNAs were also significantly correlated with key genes. AC117383.1 was positively correlated with ligand-dependent corepressor LCOR, AC144450.1, with ceramide synthase 6 LASS6, RP11-15H20.6 with zinc finger protein ZNF431, PCOLCE-AS1 with RNA-binding protein RBM5, and ITGA9-AS1 with C-terminal domain small phosphatase-like protein CTDSPL. The correlation analysis may demonstrate the important function of these eight lncRNAs and give us a clue to study their working mechanism.

Discussion
Tamoxifen therapy for five years reduced the annual breast cancer death rate by 34% and an absolute reduction in mortality of 9.2% at 15 years [2]. Despite that, there were 30% ERα+ breast cancer patients relapsed resulting from tamoxifen resistance [22]. In our work, we focused on ERα+ breast cancer, which accounted for three quarters of all breast cancer cases. We used the ERα+ breast cancer cells MCF-7/WT and its derived tamoxifen resistant cells MCF-7/TamR to investigate ERα-regulated and tamoxifen resistance-related lncRNAs. By analyzing the ERα ChIP-seq data of MCF-7/WT-E2 and MCF-7/TamR-E2 cells and RNA-seq data of MCF-7/WT and MCF-7/TamR cells, we may further understand the molecular mechanism of tamoxifen resistance and find more diagnosis and prognosis markers.    BioMed Research International We supposed whether the ERα signaling in MCF-7/TamR cells was E2-dependent, and we analyzed the ERα ChIP-seq data of MCF-7/TamR cells without E2 stimulation and with E2 stimulation. We found that the ERα directly regulated genes in MCF-7/TamR cells had two-thirds overlap of MCF-7/TamR-E2 cells, and in contrast, the ERα directly regulated genes in MCF-7/TamR-E2 cells just had a quarter overlap of MCF-7/TamR cells, which indicated that ERα signaling in MCF-7/TamR cells was also E2-dependent. We wanted to find the differential signaling pathway that ERα regulated in MCF-7/WT and MCF-7/TamR cells. We did the GO enrichment analysis of the ERα binding genes and compared their differential functional annotation. A differential enrichment of GO terms included that autophagy and I-κB kinase/NF-κB signaling pathway, and these terms were demonstrated to take an important part in tamoxifen resistance [23,24]. We proposed the upregulated lncRNAs in MCF-7/TamR cells may result from the ERα binding and the downregulated lncRNAs in MCF-7/TamR cells may due to the loss of ERα binding. So we analyzed the RNA-seq data of MCF-7/WT and MCF-7/TamR cells and made an overlap of the ERα binding genes and misregulated genes in MCF-7/TamR cells. But the overlap just made up 5.69% and 5.22% of these genes, which meant other signaling pathways except ERα signaling also played an important role in tamoxifen resistance. Despite that, we demonstrated the positive relation between the ERα status and the 143 overlap lncRNAs by analyzing the TCGA data [20]. These above steps and further survival analysis helped us get eight candidate lncRNAs. Furthermore, we did Spearman's rank correlation analysis of these eight lncRNAs and may demonstrate their important function in cellular activities. This bioinformatics process facilitated the selection of most possible tamoxifen resistance-related lncRNAs.
To our knowledge, these eight lncRNAs have not been previously studied, and they may be promising tamoxifen resistance diagnosis or prognosis markers. Further understanding the function of the eight lncRNAs will help the clinician to early diagnose whether the patients are tamoxifen resistant and bring some clinical indications in the development of novel prognostic factors in breast cancer.

Conclusion
In this study, we show the differential ER distribution events and differential expression profile in tamoxifen-sensitive and tamoxifen-resistant breast cancer cells and further demonstrate that the ER signaling pathway in tamoxifen-resistant breast cancer cells is also E2-dependent. Besides, we find eight candidate lncRNAs that may serve as diagnosis or prognosis markers.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments
The study was supported by the Central Government Guides the Development of Local Science and Technology Special Funds of China (Z135050009017).