Expression profiles and functional prediction of long non-coding RNAs LINC01133, ZEB1-AS1 and ABHD11-AS1 in the luminal subtype of breast cancer

Background Luminal breast cancer (BC) is the most frequent subtype accounting for more than 70% of BC. LncRNAs, a class of non-coding RNAs with more than 200 nucleotides, are involved in a variety of cellular processes and biological functions. Abberant expression is related to the development of various cancers, such as breast cancer. LINC01133, ZEB1-AS1, and ABHD11-AS1 were reported to be dysregulated in different cancers. However, their expression level in luminal BC remains poorly known. The aim of the present study was to evaluate the potential roles of these lncRNAs in BC, especially in luminal subtypes. Methods A comprehensive analysis was performed using the Lnc2Cancer database to identify novel cancer-associated lncRNA candidates. After conducting a literature review, three novel lncRNAs named LINC01133, ZEB1-AS1, and ABHD11-AS1 were chosen as target genes of the present study. Quantitative real‐time polymerase chain reaction (qRT-PCR) was used to evaluate the expression level of the mentioned lncRNAs in both luminal BC tissues and cell lines. Then, the correlation of the three mentioned lncRNAs expression with clinicopathological characteristics of the patients was studied. Moreover, several datasets were used to discover the potential roles and functions of LINC01133, ZEB1-AS1 and ABHD11-AS1 in luminal subtype of BC. Results According to the qRT-PCR assay, the expression levels of LINC01133 and ZEB1-AS1 were decreased in luminal BC tissues and cell lines. On the other hand, ABHD11-AS1 was upregulated in the above-mentioned samples. The expression levels of LINC01133, ZEB1-AS1, and ABHD11-AS1 were not associated with any of the clinical features. Also, the results obtained from the bioinformatics analyses were consistent with qRT-PCR data. Functional annotation of the co-expressed genes with the target lncRNAs, protein–protein interactions and significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways across luminal BC were also obtained using bioinformatics analysis. Conclusions Taken together, our findings disclosed the dysregulation of LINC01133, ZEB1-AS1, and ABHD11-AS1 in luminal BC. It was revealed that LINC01133 and ZEB1-AS1 expression was significantly downregulated in luminal BC tissues and cell lines, while ABHD11-AS1 was upregulated considerably in the mentioned tissues and cell lines. Also, bioinformatics and systems biology analyses have helped to identify the possible role of these lncRNAs in luminal BC. However, further analysis is needed to confirm the current findings. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-021-03026-7.


Introduction
Breast cancer (BC) is the most common cancer among women worldwide, accounting for 30% of all female cancers alone [1,2]. It is also the second leading cause of cancer death with 11.6% of the total cancer deaths [3]. The luminal subtype of BC, expressing both estrogen and progesterone receptors (ER + /PR + ), accounts for a large percentage of this cancer (more than 70%) [4,5]. Luminal type of BC is classified into two groups, luminal A and luminal B, according to human epidermal growth factor receptor 2 (HER2) status and levels of ki-67 [5]. The leading choice for luminal BC treatment is endocrine therapy. However, resistance to endocrine therapy is a major challenge for clinicians [6]. Also, similar treatment strategies can provide a variety of responses in luminal subtype patients [7]. In general, despite the better prognosis, due to the heterogeneity and resistance to hormone therapy of luminal subtype, the treatment effectiveness is still limited [7]. Therefore, identifying novel biomarkers and molecular mechanisms of BC is required to perform early diagnosis and personalized therapeutic strategies. In addition, the limited number of studies on the dysregulation of long non-coding RNAs (lncRNAs) in luminal BC, as well as the high frequency of this subtype, increase the importance of further research on this subgroup.
The GENCODE project results have indicated that only about 2% of the human genome encodes proteins, whereas the vast majority of it is transcribed into noncoding RNAs [8]. LncRNAs, a class of non-coding RNAs with more than 200 nucleotides, have been shown to play essential roles in many different cellular processes, such as genetic imprinting, transcriptional responses, development, etc. [9,10]. Therefore, the aberrant expression of lncRNAs can lead to the pathogenesis of various diseases including cancers [11]. They have been identified as crucial molecules involved in cancer proliferation, invasion, and resistance to therapy [12,13]. Some lncRNAs are expressed highly specific in different cancer types, and some of them are associated with the survival of patients. These two features make them ideal prognostic and diagnostic tools [14]. Due to the development of lncRNAbased therapeutic approaches, determining the role of different lncRNAs in the tumorigenic process is of great importance. Hence, in the present study, we aimed to figure out the expression level and the underlying molecular mechanisms of LINC01133, ZEB1-AS1, and ABHD11-AS1 in luminal BC.
LINC01133, a novel lncRNA, is located on chromosome 1q23.2. It was reported that LINC01133 exerted a tumor suppressive role in gastric cancer [15], ovarian cancer [16], and esophageal squamous cell carcinoma [17] progression. It was also found that LINC01133 played oncogenic role in the development of some other cancers [18][19][20]. These dual roles of LINC01133 could attribute to the tissue-specific expression of this lncRNA. ZEB1-AS1 is associated with the progression and development of several cancers. Its aberrant expression was detected in many tumors [21], such as gastric cancer [22], hepatocellular carcinoma [23], and glioma [24]. ABHD11-AS1, which is located at 7q11.23, is a newly discovered lncRNA. Previous studies have reported the upregulation of this lncRNA in different tumors, including papillary thyroid carcinoma [25] and pancreatic cancer [26]. To the best of our knowledge, the present research is the first study that determines the dysregulation of ABHD11-AS1 in BC so far. Despite all the studies that have been done on these lncRNAs, their biological processes and potential mechanisms, specifically in luminal subtype of BC, are still unclear.
In the present study, we examined the expression level of three lncRNAs in luminal BC by qRT-PCR assay. Furthermore, we performed a correlation analysis between the expression of these lncRNAs and clinicopathological parameters of patients. Finally, to investigate the potential role and mechanism of the lncRNAs in luminal BC, different bioinformatics and systems biology analyses were accomplished using various databases.

Identification of novel lncRNAs and performing literature mining
Lnc2Cancer database was used to discover novel lncR-NAs that are associated with different types of cancers.
Lnc2Cancer is a comprehensive database which provides lncRNA-cancer correlations between 165 different cancer types and 1614 lncRNAs [27,28]. Then, to find lncR-NAs that have been less or not studied in breast cancer, a comprehensive PubMed data mining was performed on the candidate lncRNAs. The following keywords were searched: "long non-coding RNA" or "lncRNA", "breast cancer" or "breast carcinoma", and "the candidate lncR-NAs". LncRNAs whose expression and role in luminal BC were not studied were selected. Also, among the selected lncRNAs, those that played an important role in several cancer-associated signaling pathways were included in the study.

In-vitro analysis Breast cancer tissue samples
To perform the present case-control study, luminal A and B BC tissue samples and paired adjacent non-cancerous tissues were obtained from 79 patients at Breast Cancer Research Biobank (BCRC-BB) [29]. BCRC-BB follows every ethical guideline about storing and utilizing human biological samples for biobanks. Tissue specimens were snap-frozen in liquid nitrogen and stored at − 80 °C until performing tests. The inclusion criterion was that patients should not receive any local or systemic treatment before surgery. Clinicopathological parameters including age at diagnosis, tumor size, tumor subtype, grade, stage, pathology of tumors, and ER/PR/HER2 status were retrieved from the hospital records. Written consent was signed by all of the patients, and all the samples were approved by relevant pathologists. Also, the Ethics committee of Tehran University of Medical Sciences (TUMS) approved the present research (Code: IR.TUMS.MEDICINE.REC.1398.792).

Isolation of total RNA and qRT-PCR
Total RNA isolation from breast tissues and cultured cell lines was performed by RiboEx ™ (GeneAll), applying the manufacturer's instructions. Then, complementary DNA (cDNA) synthesis was carried out by performing reverse transcription process of 1 µg of isolated RNA samples using 5× All-In-One RT MasterMix kit (Applied Biological Materials). Next, q-RT PCR was carried out by using AMPLIQON Real Q Plus 2× Master Mix Green Low ROX on LightCycler96 Roche system. Moreover, using serially diluted samples, the efficiencies for each pair of primers of the genes of interest were calculated. The expression levels of the genes of interest were normalized to the expression level of Beta-2-Microglobulin (β2M) as an internal control. The primer sequences are listed in Table 1.

Bioinformatics analysis
Several web servers and databases were used to select the genes of the current study and their role in luminal BC ( Table 2).

Gene correlation analysis
While lncRNAs have constituted a rapidly expanding field of research, only a small number of them has been identified functionally thus far. Analyzing the mechanisms of co-expressed genes with lncRNAs can be very useful in predicting the potential function of novel RNAs. For this purpose, the UCSC Xena Browser [30] was used to select a list of luminal samples. The following parameters were included: "Luminal A", "Luminal B", and "The TCGA Sample ID". Then, the TCGA ID of the specimens   was used to search the FPKM file of the selected samples in the Genomic Data Commons (GDC) database [31]. Finally, a correlation analysis was performed to obtain the co-expressed genes with LINC01133, ZEB1-AS1, and ABHD11-AS1 across the selected samples. Pearson correlation coefficient for the co-expressed genes has been computed by a standard method (R ≥ 0.4).

Expression analysis among molecular subtypes of BRCA
The expression pattern of the lncRNAs LINC01133, ZEB1-AS1 and ABHD11-AS1 was further evaluated using TANRIC database [32] across TCGA Breast invasive carcinoma (BRCA) dataset. The analysis was performed among different BC subtypes of the PAM50 classification (Normal-like, Luminal A, Luminal B, Basal and HER2 + ). TANRIC database has characterized the expression pattern of several lncRNAs using TCGA and different independent datasets (> 8000 samples overall). In addition, the expression level of the mentioned lncR-NAs was evaluated among luminal and also non-luminal cell lines using GENEVESTIGATOR database [33]. GEN-EVESTIGATOR integrates thousands of RNA-Seq and microarrays data in order to evaluate gene expression among several cancers.

Oncogenomic analysis
The International Cancer Genome Consortium (ICGC) portal [34] was used to identify genomic alterations of the selected lncRNAs across luminal BC. The mutational data was obtained from BRCA-EU project including ER + HER2 − tumor subtypes (luminal A and B groups). Also, figures of the lncRNAs were illustrated by trackViewer R Package [42] according to the location of the mutations.

Functional enrichment analysis
To investigate the associated pathways and functions of LIN01133, ZEB1-AS1 and ABHD11-AS1, functional enrichment analysis of the genes co-expressed with these three lncRNAs was performed using various databases. First, the Database for Annotation, Visualization and Integrated Discovery (DAVID) [35,36] was used for Gene Ontology (GO) term enrichment analysis. To this end, the list of co-expressed genes with LINC01133, ZEB1-AS1 and ABHD11-AS1 was used which was obtained in the previous step. Then, to remove redundant GO terms based on similar measures, the list of GO terms obtained from DAVID was submitted to the REVIGO [37]. At this level, Cytoscape software [38] was utilized for performing a detailed visualization of the networks specified by REViGO. The enrichr database [41] was also used to evaluate pathways in which three lncRNAs may be involved. Calculating and visualizing protein-protein interactions (PPI) of the co-expressed genes with LINC01133, ZEB1-AS1 and ABHD11-AS1 were done using the STRING database [40]. Also, the hub genes of the lncRNAs coexpressed genes were explored by the degree method available in CytoHubba [39], a Cytoscape application.

Statistical analysis
IBM SPSS Statistics (Statistical Package for the Social Sciences, Version 24.0) was used for analyzing the data obtained from the qRT-PCR assay. The tumoral and adjacent non-tumoral tissues were extracted from paired sources. The differences between ∆Ct values of LINC01133, ZEB1-AS1, and ABHD11-AS1 expression levels in tumor and adjacent non-tumor tissue specimens were analyzed using paired sample t-test. For clinicopathological correlation analysis, the median expression of the lncRNAs was used for dividing 79 patients into two categories: relatively high and relatively low expression of the lncRNAs. Also, correlation analysis of LINC01133, ZEB1-AS1, ABHD11-AS1 and clinicopathological features was carried out by the chi-square and independent t-test. Moreover, Kaplan-Meier method and log-rank test were used to calculate the overall survival rates and compare differences between survival curves, respectively. The data is represented as mean ± standard deviation (SD). P-values of less than 0.05 were considered significant. The diagnostic value of each gene was calculated via the receiver operating characteristic (ROC) curve analysis. ROC curves were illustrated by IBM SPSS Statistics version 24.0.

Evaluation of diagnostic value of LINC01133, ZEB1-AS1 and ABHD11-AS1 in luminal BC
The outcome of ROC analysis was defined based on the type of tissue (tumor or adjacent) and the predictor was characterized as gene expression level (upper or below the median expression). According to the ROC analysis of ZEB1-AS1(T5,6) shown in Fig. 3, the area under ROC curve is 0.729 across luminal A and B subtypes of BC (p ZEB1-AS1(T5,6) < 0.001). On the other hand, LINC01133, ZEB1-AS1(T1-4), and ABHD11-AS1 were not suggested as good biomarkers for diagnosis of luminal subtypes of BC (p

Expression pattern of lncRNAs LINC01133, ZEB1-AS1 and ABHD11-AS1 among PAM50-breast cancer subtype classifications
According to the differential expression analysis, lower expression of LINC01133 was observed in luminal A and B BC subtypes. The expression of ZEB1-AS1 in luminal B and HER2 subtypes was lower than other groups. Moreover, significant differences of the expression of these lncRNAs were observed among different PAM50 subtypes (p LINC01133 = 1.13e−32) (p ZEB1-AS1 = 9.82e−5).
On the other hand, luminal A and B had higher expression of ABHD11-AS1 followed by Normal-like, HER2, and Basal subtypes. The expression of ABHD11-AS1 was also significantly different among molecular subtypes of BC (p ABHD11-AS1 = 2.79e−8) (Additional file 2). Besides, results obtained from the GENEVESTIGA-TOR database across luminal BC cell lines were consistent with the experimental data of the present study which indicated that LINC01133 and ZEB1-AS1 were downregulated in luminal cell lines. Also, expression analysis of these lncRNAs in several non-luminal BC cell lines showed that their expression was also downregulated across these cell lines (Additional files 3, 4). However, as ABHD11-AS1 is a novel lncRNA, no information was available about this lncRNA across GENEVESTIGATOR database. According to the results obtained from the ICGC data portal, LINC01133 mutations in ER + HER2 -BC take place most frequently in intronic, downstream, and upstream regions, respectively. Also, most LINC01133 mutations have been identified to be of the substitution type (Fig. 4A). It was also found that ZEB1-AS1 mutations occur more frequently on introns and exons and are of the substitution type (Fig. 4B).
Most of ABHD11-AS1 mutations are substitutions which mostly take place in upstream and downstream regions (Fig. 4C).

Co-expression gene network and functional annotation analysis
A list containing 420 luminal samples (230 luminal A and 190 luminal B) was obtained using UCSC Xena Browser and the Genomic Data Commons (GDC) databases (Additional file 5). Results obtained from the co-expression analysis of LINC01133, ZEB1-AS1, and ABHD11-AS1 across the selected luminal A and B BC samples suggested that these lncRNAs are significantly co-expressed with 256, 149 and 126 genes, respectively (R ≥ 0.4) (Additional file 6).
Data retrieved from GO term enrichment analysis by DAVID included a list for categories of biological process (BP), cellular component (CC), and molecular function (MF) (Additional file 7). Also, results of summarizing and excluding redundant GO terms of the target lncRNAs of the present study, which were performed by REVIGO, were concordant with the data obtained from DAVID GO term enrichment analysis. The top 5 BP, CC, and MF terms related to each lncRNA are shown in Fig. 5. According to it, the genes co-expressed with LINC01133 are significantly involved in biological processes such as ion transport, response to external stimulus, ion transmembrane transport, positive regulation of phosphorylation, etc. It was also determined that molecular function of these genes considerably consists of ion/cation/metal binding, cation transmembrane transporter activity, etc. Each of the above activities takes place inside and/or in the vicinity of a cell. According to CC terms, it was predicted that the co-expressed genes with LINC01133 are considerably found in extracellular region/exosome, membranebounded vesicle and so on (p < 0.05) (Fig. 5A).
It has been found that the co-expressed genes with ZEB1-AS1 are considerably involved in the regulation of RNA metabolic process, regulation of nucleobasecontaining compound metabolic process, nucleic acidtemplated transcription, etc. As shown in Fig. 5B, the most significant molecular functions of these genes are signal transducer activity, kinase activity, phosphotransferase activity, etc. These genes are mostly enriched in the synapse and Golgi lumen (p < 0.05).
Moreover, the genes co-expressed with ABHD11-AS1 are considerably involved in biological processes such as transmembrane transport, reproductive process, ion transport, neurological system process, etc. They are also considerably involved in various molecular functions including signal transducer activity, transmembrane signaling receptor activity, molecular transducer activity, etc. Additionally, the GO cellular component terms suggested that the co-expressed genes with ABHD11-AS1 are considerably enriched in the intrinsic component of nuclear membrane part, apicolateral plasma membrane, etc. (p < 0.05) (Fig. 5C).
The results of gene set enrichment analysis by Enrichr are demonstrated in Fig. 6. According to KEGG, Reactome and WikiPathways analysis, the co-expressed genes with LINC01133 are mainly involved in sodium/ proton exchanges, mucin-type O-glycan biosynthesis, and nuclear factor erythroid 2-related factor 2 (NRF2) pathway (Fig. 6A). Moreover, ZEB1-AS1 co-expressed genes are significantly enriched in tight junction interactions, estrogen biosynthesis, glycerophospholipid catabolism, and Wnt signaling pathway across luminal BC (Fig. 6B). The co-expressed genes with ABHD11-AS1 are mainly involved in the processes of Wnt ligand biogenesis, FGFR1 mutant receptor activation, mesodermal commitment pathway, and Wnt signaling (Fig. 6C). Finally, the association analysis of the genes coexpressed with LINC01133 demonstrated that 42 out of 256 genes had a strong PPI (interaction score > 0.4) with each other (Fig. 7A). Also, according to PPI network analysis of the co-expressed genes with ZEB1-AS1 and ABHD11-AS1, 20 and 14 genes have been recognized with strong interactions with each other, respectively (Fig. 7B, C).

Discussion
The luminal subtype of BC, which has shown better prognosis among the other subtypes, accounts for approximately two-thirds of this cancer [43,44]. However, it is reported that the luminal B subtype has demonstrated poorer recurrence-free survival in adjuvant treatment categories in comparison to luminal A [44]. Also, it was proved that a large percentage of ER + patients with lymph node positive benefit less from adjuvant chemotherapy [43]. Due to the heterogeneity of this cancer, the same therapies may demonstrate diverse outcomes in patients. So, distinguishing the underlying mechanisms and specific biomarkers can help the individualized treatment of especially highrisked luminal patients.
In this study, LINC01133 and ZEB1-AS1 were significantly downregulated in luminal A and B BC tissues compared to their adjacent non-tumoral tissues. Also, they were downregulated in luminal cell lines, namely T47D and MCF7, compared to a normal breast cell line (MCF10A). Moreover, ABHD11-AS1 was upregulated considerably in the mentioned tissues and cell lines.
LncRNAs have shown greater expression variability among different individuals and a greater degree of tissue-specificity, compared with coding-genes [45]. The downregulation of LINC01133, ZEB1-AS1(T1-4), and ZEB1-AS1(T5,6) was detected in 64%, 71% and 80% of patients, respectively. Also, the upregulation of ABHD11-AS1 was found in a large proportion of patients (70% of cases). Song et al. reported that the downregulation of LINC01133 is considerably associated with the poor prognosis of BC patients [46]. On the other hand, Luo et al. suggested that ZEB1-AS1 promotes triple-negative breast cancer progression [47]. So, this lncRNA may also play an important role in luminal BC patients, although further investigation is needed. Also, ZEB1-AS1(T5,6) may be used as a biomarker for luminal BC progression. However, more invitro tests are needed for more accurate conclusions.
Limited studies have been performed on the dysregulation of these lncRNAs and their roles in cancers so far. So, the combination of various datasets as well as the analysis of their co-expressed genes could provide better identification of functional roles of the mentioned lncR-NAs in luminal BC. According to the result of TCGA data analysis available in TANRIC database, LINC01133, ZEB1-AS1, and ABHD11-AS1 expression were significantly different among PAM50 subtype classification. Besides, the expression pattern of LINC01133 and ZEB1-AS1 across luminal and non-luminal BC cell lines was consistent with the data from TANRIC database and also experimental data of the present study. According to the data obtained from the ICGC data portal, most of the mutations of LINC01133, ZEB1-AS1, and ABHD11-AS1 in ER + HER2patients are from substitution category. The effect of mutations and variations of lncRNAs have not been extensively studied as they constitute a novel class of RNAs. Ponjavic et al. demonstrated that functional lncRNAs show reduced single substitutions, deletions and insertions in their sequences. Although nucleotide substitutions occur mostly in protein-coding sequences compared with noncoding sequences, their effect on lncRNAs are so significant and should not be overlooked [48]. Substitution mutations in lncRNAs may contribute to the pathogenesis of various cancers by disrupting the secondary structure of RNAs [18,49]. The RNA secondary structure plays important roles in multiple cellular processes such as gene regulation and localization, splicing, stability and also translation [50]. Also, mutations in the upstream regions of lncRNAs, which constitute a considerable percentage of LINC01133 and ABHD11-AS1 mutations, can have significant effects on them. For example, the presence of single-nucleotide polymorphisms in the promoter of a lncRNA may have an effect on its expression pattern [49]. So, further studies are recommended to determine the effect of these mutations on the structure of the mentioned lncRNAs.
Remarkably, the protein-coding genes are significantly better annotated than lncRNAs [51]. So, in this study, an attempt was made to determine the function of the target lncRNAs with the help of the known functional information of some coding and non-coding genes, using "guiltby-association" principle. According to this principle, Fig. 4 Genetic alterations within 5 kb upstream and downstream of the candidate lncRNAs across ER + HER2 -BC, studied in ICGC. A LINC01133 mutations occur more frequently in intronic, downstream, and upstream regions, respectively. They are mostly of the substitution type. B ZEB1-AS1 mutations occur more frequently in intronic and exonic regions and are mostly of the substitution type. C ABHD11-AS1 mutations occur more frequently in upstream and downstream regions, respectively. ABHD11-AS1 mutations are also mostly of the substitution type. All of the data were obtained from ICGC and were illustrated by trackViewer R Package  19:364 genes that are involved in some related and/or similar biological pathways may show similar expression patterns across different experimental conditions [52]. According to the "guilt-by-association" principle and the GO term enrichment analysis of the co-expressed genes with LINC01133, this lncRNA might be involved in ion transport, response to external stimulus, and positive regulation of phosphorylation, all of which are cancerassociated biological processes. Also, ion/cation/metal binding and cation transmembrane transporter activity are known as the most crucial molecular functions in which LINC01133 is involved. LINC01133 might be mostly found in extracellular region/exosome and membrane-bounded vesicles.
According to the results obtained from the GO term enrichment analysis of ZEB1-AS1 co-expressed genes, this lncRNA is involved in the regulation of RNA metabolic process, regulation of nucleobase-containing compound metabolic process, nucleic acid-templated transcription, etc. Deregulation of all of these biological processes can lead to the development and progression of various cancers, including BC. In addition, ZEB1-AS1 might be plausibly involved in molecular functions including signal transducer activity, kinase activity and phosphotransferase activity, all of which can be linked to cancer progression. ZEB1-AS1 might be mostly enriched in synapse and Golgi lumen.
ABHD11-AS1 might be involved in transmembrane transport, reproductive process, ion transport, neurological system process, etc. ABHD11-AS1 also might be considerably involved in different molecular functions, like signal transducer activity, transmembrane signaling receptor activity, molecular transducer activity, etc. Dysregulation of all of the mentioned functions can lead to the development of various cancers. Moreover, ABHD11-AS1 is likely to be found in intrinsic component of nuclear membrane part, apicolateral plasma membrane, etc.
Significantly enriched KEGG, Reactome and WikiPathways of the co-expressed genes with LINC01133 include sodium/proton exchanges, mucin-type O-glycan biosynthesis, and NRF2 pathway. Abnormal glycoprotein structure of tumor cells can affect the survival, growth and metastasis of cancer cells [53]. Nuclear receptors (NR) are known as regulators of physiological processes and play oncogenic or anti-oncogenic roles in cancerous cells [54]. Growing evidence support the involvement of several NRs in the regulation of various pathways related to the initiation and development of BC [55]. Overexpression of Nrf2 has been shown to increase the expression of glucose-6-phosphate dehydrogenase (G6PD) and Hypoxia-inducing factor 1α (HIF-1α) in BC cell lines, including MCF7. Also, overexpression of Nrf2 increases Notch1 expression via the G6PD/HIF-1α pathway. Notch signaling pathway regulates BC cell proliferation and migration by affecting the downstream gene, HES-1, and the epithelial-to-mesenchymal transition (EMT) pathway, respectively [56]. All of these evidences suggest that LINC01133 could play an essential role in BC.
Moreover, according to the pathway enrichment analysis, ZEB1-AS1 might be involved in tight junction interactions, estrogen biosynthesis, glycerophospholipid catabolism, and Wnt signaling pathway. All of these pathways can be associated with the progression and development of BC. Steroid hormones increase cell proliferation in BC [57]. Also, alterations in tight junction complexes could facilitate BC initiation and progression by impairing their control over crucial cellular processes such as cell motility and polarity [58]. Moreover, the Wnt pathway is significantly activated in breast tumors [59].
ABHD11-AS1 might also be mainly involved in Wnt ligand biogenesis, FGFR1 mutant receptor activation, mesodermal commitment pathway, and Wnt signaling pathway. The involvement of fibroblast growth factors (FGF) has been found in various cellular processes, including proliferation, drug resistance, anti-apoptosis and angiogenesis. Also, amplification of FGFR1 has been found in ER + BC [60].
Data from STRING database showed that 42, 20, and 14 co-expressed genes with LINC01133, ZEB1-AS1 and ABHD11-AS1 have strong interactions with each other. These data confirmed the involvement of these lncR-NAs and their co-expressed genes in similar biological pathways.
Taken together, the dysregulation of the three potential lncRNAs, LINC01133, ZEB1-AS1, and ABHD11-AS1 (See figure on next page.) Fig. 5 Top 5 gene ontology (GO) terms related to the co-expressed genes with the candidate lncRNAs for categories of biological process (BP), cellular component (CC), and molecular function (MF). A Bar chart of the GO term enrichment analysis of the co-expressed genes with LINC01133. B Bar chart of the GO term enrichment analysis of the co-expressed genes with ZEB1-AS1. C Bar chart of the GO term enrichment analysis of the co-expressed genes with ABHD11-AS1. The top GO terms were selected based on the amount of gene count and the chart was illustrated according to the data obtained from enrichment analysis. The columns were sorted from top to bottom according to significancy of the terms (p < 0.05). All of the data was obtained from DAVID database across luminal A and B subtypes of BC was reported in the present study. The bioinformatics analyses performed in this study helped us better identify the possible role of these lncRNAs in luminal BC. However, more experimental studies are needed to confirm these findings and to verify the exact roles of the mentioned lncRNAs.

Conclusions
All our findings have demonstrated the downregulation of LINC01133 and ZEB1-AS1 across luminal BC tissue specimens and cell lines compared to adjacent non-tumoral tissues and cell line. On the other hand, ABHD11-AS1 was upregulated in the mentioned BC tissues and cell lines. The pathway analysis of the hub genes showed involvement of the co-expressed genes with the candidate lncRNAs in cancer-related pathways. A Bar chart of the KEGG, Reactome and WikiPathways enrichment analysis of the co-expressed genes with LINC01133. B The Reactome pathway enrichment analysis of the co-expressed genes with ZEB1-AS1. C The Reactome and WikiPathways enrichment analysis of the co-expressed genes with ABHD11-AS1. The columns were sorted from top to bottom according to significancy of the terms (p < 0.05). All of the data was obtained from Enrichr web server Fig. 7 The analyzed protein-protein networks of the co-expressed genes with the target lncRNAs confirmed strong interactions between a number of genes. A 42 genes co-expressed with LINC01133 have strong interactions with each other. B 20 genes co-expressed with ZEB1-AS1 have strong interactions with each other. C 14 genes co-expressed with ABHD11-AS1 have strong interactions with each other. The thicker edge and the stronger node color are representative of the higher STRING combined-score and the higher degree, respectively. Data was obtained from the STRING database and was illustrated by Cytoscape (interaction score > 0.4)