Introduction
Long non-coding RNAs (lncRNAs) are encoded similarly to coding genes but do not contain a protein-coding sequence. The transcripts of lncRNA are >200 nucleotides long and expressed at lower levels than protein-coding transcripts. Some lncRNAs directly modulate gene expression in a cis manner; these lncRNAs may be expressed from the promoter (1), intron (2) and enhancer regions (3) of certain genes and then regulate neighboring protein-coding genes on the same chromosome. Furthermore, some lncRNAs distally regulate gene expression across multiple chromosomes in a trans manner. These lncRNAs can facilitate enhancer function via long-range DNA looping interactions (4). Also, some lncRNAs appear to modulate the DNA-binding of certain transcription factors (TFs) and non-DNA-binding cofactors at several target sites, thus affecting gene expression in trans (5).
In recent years, several lncRNAs have been implicated in cancer development and progression (6–8). Chemoresistance is an important feature of cancer progression, since chemoresistant cells are insensitive to the apoptotic signals delivered by cytotoxic chemotherapeutic agents. They also show a strong ability to proliferate. A large degree of chemoresistance is acquired during the response of cancer cells to unfavorable niches, and their rapid response depends on the effective and heritable epigenetic regulation of gene expression, including DNA methylation and microRNA regulation (9). For example, methylation of microRNA-200c inhibits microRNA expression, so the targets of the microRNA, including modulators of epithelial mesenchymal transition pathways that essentially promote cell proliferation and diminish apoptosis, are activated (10–14). In the current decade, emerging evidence indicates that 70–90% of the mammalian genome produces lncRNAs, which are another important component of epigenetic regulation (15,16). However, although >10,000 mammalian lncRNAs have been catalogued (5), functional studies are as yet limited in number, and how the network of lncRNAs is involved in cancer chemoresistance still requires elucidation.
In the present study, we set out to analyze lncRNA expression in adriamycin-resistant breast cancer cells using microarrays, and compare their lncRNA expression profile with that of parental chemosensitive cells in order to identify and characterize dysregulated lncRNAs that may be involved in breast cancer chemoresistance (Fig. 1A).
Materials and methods
Cell culture
MCF-7/WT human breast cancer cells were obtained from the American Type Culture Collection (ATCC; Manassas, VA, USA). Adriamycin (ADM)-resistant cells (MCF-7/ADM) and paclitaxel (PTX)-resistant cells (MCF-7/PTX) were derived by treating MCF-7 cells with stepwise increasing concentrations of ADM or PTX over 8 months (17). The cells were cultured in RPMI-1640 supplemented with 10% FBS, 100 μg/ml penicillin and 100 U/ml streptomycin.
Microarray
The microarray profiling was conducted in the laboratory of OE Biotechnology Co. (Shanghai, China). RNA from MCF-7/ADM and MCF-7/WT cells was separately extracted using the acid-phenol and chloroform method. Cyanine-3-CTP-labeled cRNA was obtained using a Quick Amp Labeling kit (Agilent Technologies, Santa Clara, CA, USA) and then purified with an RNeasy Mini kit (Qiagen, Valencia, CA, USA). The labeled cRNAs were then hybridized onto Agilent-062918 OE Human lncRNA Microarray V4.0 028004 (Agilent Technologies), which is a Custom Gene Expression Array for OE Biotechnology Co. and detects 46,506 lncRNAs. After washing, the arrays were scanned with an Agilent scanner (G2505C).
Quality control of the microarray data
Total RNA was quantified by the NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and the RNA integrity was assessed using Agilent Bioanalyzer 2100 (Agilent Technologies). RNA with RNA integrity number (RIN) value >7 and 28S/18S >0.7 was used for microarray analysis. After the raw data extraction, the median CV (%) for each probe set was reported as the array reproducibility. To analyze the biological repeatability of the microarrays between MCF-7/ADM (n=3) or MCF-7/WT (n=3) cells, two-dimensional principal component analysis was performed.
Data deposition
The microarray data have been submitted to GEO (GSE81971).
Data analysis
Raw data of microarray was generated using Agilent Feature Extraction software (Agilent Technologies) and then normalized using GeneSpring’s quantile normalization (version 12.5; Agilent Technologies). Differentially-expressed lncRNAs were identified with a fold change ≥2.0 and a P-value <0.05. WebGestalt (http://bioinfo.vanderbilt.edu/webgestalt/) was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.
Real-time PCR
LncRNA expression was analyzed using qRT-PCR. Briefly, total RNA was extracted from MCF-7/WT and /ADM cells with TRIzol (Invitrogen, Carlsbad, CA, USA). cDNA was synthesized from total RNA (3 μg) using the SuperScript First Strand Synthesis system (Invitrogen) with Oligo (dT) primers. Primers used for real-time PCR were as previously described for UCA1 (18,19), HOTAIR (6), GAS5 (20), HULC (21), MEG3 (22), SPRY4-IT1 (23) and CRNDE (24). Primer for NONHSAT028712 was: forward, 5′-AAATACCTCACCCTCATCTATACCAAC-3′ and reverse, 5′-TTTCCCGTTGCCATTGAT-3′; for NONHSAT057282 was forward, 5′-AGCCGGAGGTGAGGAAGTT-3′ and reverse, 5′-AAGATTTTATTAGATTTTGGAACCTGAG-3′; for NONHSAG023333 was forward, 5′-GTTGGGAAATCAAG CATCGT-3′ and reverse, 5′-TTTAGCAAAAATGCAACTA CATCC-3′.
The RT-PCR values were normalized to GAPDH and calculated using the 2−ΔΔCT method.
LncRNA inhibition and functional studies
LncRNA Smart Silencer was synthesized by Guangzhou RiboBio Co., Ltd., (Guangzhou, China) and used to inhibit lncRNA. The inhibitor was then transfected into MCF-7 cells using Lipofectamine (Invitrogen).
MCF-7/ADM or MCF-7/PTX cells transfected with lncRNA inhibitor were seeded onto 96-well plates and exposed to different concentrations of ADM or PTX for 48 h. Then cell viability and IC50 were assessed as previously described (25,26).
In experiments analyzing the cell cycle, MCF-7/ADM cells transfected with lncRNA inhibitor were fixed and stained with 100 μg/ml propidium iodide (PI; Sigma Life Science) containing RNase. PI fluorescence was detected using a FACSCalibur flow cytometer on the PE-Texas Red channel for DNA content.
Results
Expression of lncRNAs in chemoresistant breast cancer cells
To gain insights into the role of lncRNAs in chemoresistance, we used microarray-based profiling to analyze the lncRNAs and mRNAs in adriamycin-resistant MCF-7/ADM cells and their chemosensitive parental control MCF-7/WT cells (17). In MCF-7/ADM cells, 4030 lncRNAs and unannoted transcripts were upregulated and 3708 were downregulated (Fig. 1B; Submitted online as GSE81971), while 3423 mRNAs were upregulated and 2950 were downregulated (Fig. 1C; fold-change ≥2, P<0.05), suggesting that lncRNAs may be dysregulated and participate in the development of chemoresistance.
We then validated several cancer-related lncRNAs with RT-PCR (6,7,18–24). Among these, SPRY4-IT1 was upregulated in chemoresistant MCF-7/ADM cells by both microarray analysis and RT-PCR, suggesting that it may play a role in chemoresistance; other lncRNAs were not significantly changed either in microarray analysis or RT-PCR (Fig. 1D). These data, thus, confirm the accuracy of microarrays, while prompted us to search for new lncRNAs that may mediate chemoresistance.
In order to identify new candidate lncRNAs in chemoresistance, we then chose the top 200 most significantly changed lncRNAs for further analysis.
LncRNAs correlate with mRNA expression
LncRNAs influence gene expression by regulating chromatin remodeling, transcription and post-transcriptional processing (5). To identify potential lncRNA targets, we calculated the Pearson correlation of each significantly changed lncRNA with each significantly changed mRNA. An lncRNA and an mRNA were considered to be correlated when the coefficient was >0.7 (P<0.05), so such mRNAs might be regulated by their correlated lncRNAs. The most correlated mRNAs for top 10 changed lncRNAs are exemplified in Table I.
The 500 lncRNA-mRNAs with the highest Pearson correlation coefficient values were chosen for functional analysis in GO and KEGG using the method described by Guttman et al (27). The GO and KEGG functions of each lncRNA-correlated mRNA were analyzed, then a hypergeometric cumulative distribution function was applied to calculate the enrichment of functional terms in the annotation of these mRNAs. The most enriched GO processes and KEGG pathways are shown in Fig. 2. Both analyses showed that pathways directly associated with apoptosis and cell proliferation were frequently regulated by lncRNAs, including release of cytochrome c from mitochondria (GO:0001836), cell cycle (GO:0007049), cell cycle arrest (GO:0007050), apoptosis (KEGG 04210), negative regulation of cell proliferation (GO:0008285), MAPK signaling pathway (KEGG 04010) and p53 signaling pathway (KEGG 04115).
Cis-regulation of lnRNAs
Cis-regulation by lncRNAs of their correlated mRNAs was then analyzed in chemoresistant MCF-7/ADM cells vs. MCF-7/WT cells. Several lncRNAs were located 100K windows upstream or downstream of the given mRNA, and the mRNA expression was significantly correlated with the lncRNA. Possible cis-regulation of their correlated mRNAs is exemplified in Table II. Among these, we found that NONHSAT028712 was significantly increased in MCF-7/ADM cells (Table II and Fig. 3A). We also analyzed the expression of NONHSAT028712 in another chemoresistant breast cancer cell line MCF-7/PTX, which is paclitaxel-resistant (9). This lncRNA also significantly increased in MCF-7/PTX cells (Fig. 3A). These data suggest a possible role of NONHSAT028712 in mediating chemoresistance.
The cis-regulation of NONHSAT028712 is shown in Fig. 3B. Genes of several significantly-changed mRNAs was found to locate near the coding sequence of NONHSAT028712. Among these mRNAs, the cell cycle kinase CDK2 showed a high mRNA level in MCF-7/ADM and MCF-7/PTX cells (Fig. 3A), and this has been associated with cancer progression and chemoresistance (28,29). We therefore inhibited the expression of NONHSAT028712 with a synthesized inhibitor (Fig. 3C), then analyzed the chemoresistance of the treated MCF-7/ADM cells. The IC50 significantly decreased in MCF-7/ADM cells when NONHSAT028712 was inhibited (Fig. 3D), along with a lower CDK2 mRNA level (Fig. 3E) and a higher rate of cell cycle arrest in G1 (Fig. 3F). These data strongly suggest that NONHSAT028712 regulates chemoresistance via a CDK2-related pathway, while the mechanism requires further exploration.
Interaction of transcription factors with lncRNAs
To identify the possible role of lncRNA-TF interactions in regulating gene expression, we first predicted the TFs of lncRNA-correlated mRNAs using data from Gerstein et al (30) that showed the genomic binding information of different TFs. Then the intersections of lncRNA-mRNA and mRNA-TF were calculated with a hypergeometric cumulative distribution. Each lncRNA was significantly associated with several TFs (data not shown). For instance, NONHSAT057282 and NONHSAG023333 were significantly correlated with ELF1 and E2F1, and enriched the most mRNAs of all lncRNA-TF interactions. NONHSAT057282-ELF1 co-regulated 241 genes and NONHSAG023333-E2F1 co-regulated 308 genes. We then drew a ternary relation graph of the two interactions with Cytoscape 3.01 software (Agilent) (Fig. 4A). Furthermore, NONHSAT05728 was upregulated in chemoresistant MCF-7/ADM and MCF-7/PTX cells vs. chemosensitive MCF-7/WT cells, while NONHSAG023333 was upregulated in MCF-7/WT vs. MCF-7/ADM and MCF-7/PTX cells (Fig. 4B), suggesting that NONHSAT05728 may enhance chemoresistance, but NONHSAG023333 may negatively regulate chemoresistance. Indeed, when we knocked down NONHSAT05728, chemoresistance decreased in both MCF-7/ADM and MCF-7/PTX cells. On the other hand, knockdown of NONHSAG023333 increased the chemoresistance in MCF-7/WT cells (Fig. 4C). Therefore, these results suggest that both NONHSAT05728 and NONHSAG023333 are involved in chemoresistance, and their activity may be facilitated by TFs.
The top 15 TFs with the highest enrichment of lncRNAs are summarized in Fig. 4D, indicating the potential involvement of certain TFs in regulating chemoresistance via lncRNAs. In order to visualize these most significantly-related lncRNA-TF interactions, the top 100 with the lowest Q-values were then used to draw a two-element relation graph (Fig. 4E). ELF1 was still most frequently associated with several lncRNAs, and PBX3 and ZEB1 were also intensively associated with lncRNAs.
Discussion
In the present study, for the first time we assessed the genome-wide lncRNA expression patterns in adriamycin resistant breast cancer cells using microarrays and explored their possible functions by analyzing their cis-regulated mRNAs, as well as TF-regulated mRNAs.
We first identified dysregulated lncRNAs in adriamycin-resistant MCF-7/ADM cells; these lncRNAs correlated with a list of dysregulated mRNAs. Because most of the lncRNAs in current databases have not yet been functionally annotated, we predicted their functions based on their correlated mRNAs. Chemoresistance is an important feature of cancer progression, so the lncRNAs dysregulated in chemoresistant MCF-7/ADM cells also showed functions associated with hallmarks of cancer progression (31). For instance, proteoglycans in cancer (KEGG 05205) are responsible for increased cancerous angiogenesis and provide a favorable microenvironment for cancer cells (32); and toll-like receptor signaling pathways (GO:0002224 and 0002756) provide cancer cells with sustained proliferative signals (33). Furthermore, the key feature of chemoresistant cancer cells is insensitivity to the cytotoxicity of chemotherapeutic agents. Such insensitivity could be achieved by a low efficacy of cellular drug transport, which may be associated with actin filament bundle assembly (GO:0051017) (34) and endocytosis (KEGG 04144) (35). Importantly, apoptosis inactivation and cell proliferation enhancement, whose pathways were frequently enriched in both the GO and KEGG pathways, not only support cancer growth and metastasis but also contribute to the insensitivity of cancer cells to chemotherapeutic agents. Therefore, being the ‘mission critical’ of cancer progression and chemoresistance regulated by various genetic and epigenetic mechanisms, we suggest apoptosis and cell proliferation may be still the main targets of regulation by lncRNAs.
Cis-action on target genes located at or near the same locus is one of the main mechanism by which lncRNAs regulate gene expression (36). We therefore identified genes whose expression was correlated with that of nearby lncRNAs. This analytical method greatly facilitated the identification of lncRNAs critical for chemoresistance. Based on the roles of lncRNAs in apoptosis and cell proliferation during chemoresistance, we then explored the possible mechanism of action of NONHSAT028712 in MCF-7/ADM and MCF-7/PTX cells because it may cis-regulate CDK2. The preliminary results strongly suggest an interaction between NONHSAT028712, CDK2, the cell cycle, and chemoresistance, and further studies are needed to clarify how NONHSAT028712 modulates expression of the nearby CDK2 gene; it may directly interact with the gene, facilitate the 3D folding of chromatin, or interact with other genetic (e.g. TFs) and epigenetic (e.g. microRNAs) regulators of the CDK2 gene (36).
LncRNAs frequently physically interact with TFs to regulate gene expression. We found that NONHSAT057282 and NONHSAG023333 were involved in chemoresistance. Then the lncRNA-TF interaction analysis suggests that these two lncRNAs may interact with ELF1 and E2F1 respectively, and subsequently modulated a group of chemoresistance-related genes such as GSTP1 (37,38), BTG3 (39), SOCS3 (40) and BRAC2 (41). Furthermore, among the identified lncRNA-TF interactions, ELF1 showed the highest enrichment frequency; 50 lncRNAs were significantly associated with this TF. ELF1 belong to the ETS transcription factor family, which is important for cancer progression (42) and breast cancer chemoresistance, as we demonstrated previously (25). ELF1 is associated with tumor angiogenesis (43), but its role in cancer chemoresistance has not been identified. Our results, thus, suggest ELF1 as a new participant in chemoresistance by potentially interacting with different lncRNAs. In future studies, mass spectrometry could be applied to confirm the lncRNA-TF interactions (44), and chromatin immunoprecipitation-based sequencing might be needed to verify the ELF1-related target genes (6). The TF ZEB1 also frequently interacted with lncRNAs. ZEB1 modulates the epithelial-mesenchymal transition pathway, which is essential for chemoresistance (45). Previous studies have shown that the epigenetic regulation of ZEB1 by microRNAs and DNA methylation effectively generates chemoresistance (10,46), but few studies have explored the ZEB1-lncRNA interaction in chemoresistance. Therefore, considering the importance of ZEB1, it would be worthwhile to investigate the molecular mechanism by which lncRNAs mediate chemoresistance via ZEB1; here, we have provided several target lncRNAs that are likely to interact with ZEB1. Furthermore, other significantly enriched TFs, such as PBX3 (47) and E2F1 (48–50) are also involved in breast cancer progression; their interactions with lncRNAs may be important mechanisms by which they control gene expression and enhance chemoresistance.
In summary, we provide an overview of lncRNA regulation at the combined levels of lncRNA and gene expression in breast cancer chemoresistance, and systematically identify novel dysregulated targets in chemoresistant breast cancer cells. Experimental validation of specific interactions between lncRNAs and genes, and lncRNAs and TFs allowed us to identify the key players in chemoresistance, and decipher the underlining molecular mechanism of action of lncRNAs in cancer progression. Based on this analytical approach, we have shown the relationship and possible mechanism of action of several lncRNAs in the development of chemoresistance, suggesting that the analysis is precise and valuable for the support of future studies.