An innovative systematic approach introduced the involved lncRNA-miR-mRNA network in cell cycle and proliferation after conventional treatments in breast cancer patients

ABSTRACT

The present study aimed to explore the involved lncRNA-miRNA-mRNA network in the cell cycle and proliferation after conventional treatments in Luminal A breast cancer patients.The candidate miRNAs (miRs), lncRNAs, and mRNAs were first taken from the Gene Expression Omnibus and TCGA databases. The lncRNA–miR–mRNA network was then constructed using the high-throughput sequencing data. The expression levels of selected targets were measured in the breast cancer and healthy samples by the Real-Time PCR technique and compared with the clinical outcomes by the Kaplan-Meier method.Our analysis revealed a group of differentially expressed 3 lncRNAs, 9 miRs, and 14 mRNAs in breast cancer patients. A significant expression decrease of the selected tumor suppressor lncRNAs, miRs, and genes and a substantial expression increase of the selected onco-lncRNAs, oncomiRs, and oncogenes were obtained in the patients compared to the healthy group. The plasma levels of the lncRNAs, miRs, and mRNAs were more significant after the operation, chemotherapy, and radiotherapy than the pre-treatment. The Kaplan-Meier analysis indicated that the patients with a high expression of miR-21, miR-20b, IGF1R, and E2F2 and a low expression of miR-125a, PDCD4, and PTEN had exhibited a shorter overall survival rate.Our results suggested that the underlying mechanisms of the lncRNA, miRs, and mRNAs and relevant signaling pathways may be considered predictive and therapeutic targets for breast cancer.

KEYWORDS:

• Breast tumor
• microRNA
• lncRNA
• mRNA
• chemotherapy
• radiotherapy
• operation

Introduction

Non-coding RNAs have central roles in regulating specific genes’ expression and may involve breast cancer’s primary biological processes. They can play vital roles in regulating various biological processes, including cell growth, apoptosis, and stem cell differentiation [1]. Thus, they can be considered potential prognostic and diagnostic biomarkers, therapeutic targets, and treatment-response predictors for human cancers. Recent studies have suggested that lncRNAs interact with microRNAs (miRs) and modulate their availability for endogenous mRNA targets. Besides, they can show a particularly tissue-specific or cell-type-specific expression pattern. Thus, elucidating the function of lncRNAs in normal and cancerous cells will be one of the following milestones for understanding the molecular mechanisms of cancer. Interestingly, it has been reported that lncRNAs can also compete with miRs by binding directly to mRNAs [2,3]. This issue suggests a complex interaction between lncRNAs and miRs that can eventually determine the stability and translation of coding protein mRNAs.

Likewise, the aberrant expression of miRs in cell proliferative and cell cycle regulatory pathways is also attributed to resistance to breast cancer treatments to suppress cell growth and proliferation. It is hypothesized that interplay between lncRNAs and miRs might comparatively elucidate the cell-type-specific outcomes. Therefore, given this evidence, it has appeared that the abundance of lncRNAs may ultimately determine the effects of miRs on the expression of protein-coding genes. Generally, the interplay between lncRNAs and miRs at the intracellular and intercellular levels may provide a framework for understanding tissue-specific phenomena in cancer. Subsequently, the interaction between lncRNAs and miRs may affect the signaling cascades by regulating the abundance of proteins in these pathways in a cell-specific manner. This putative mechanism, acting on the intracellular level, can also be extended to the communication between different cells. Hence, cell proliferation and cell cycle mechanisms may be partly explained by the interaction of miRs and lncRNAs. Although there is no direct experimental evidence for this hypothesis, theoretical evidence demonstrates how this interaction could contribute to cell proliferation in specific tissues.

Overall, monitoring the miRs, lncRNA, and target genes during the treatment can be an appropriate diagnostic tool for estimating the identification and efficacy of therapeutic approaches. However, no significant progress has been made toward screening and monitoring cancer patients [4]. Several studies have shown that interactions among lncRNAs, miRs, and target genes can be used as sensitive and specific indicators for identifying different diseases, including cancer [5]. A better understanding of underlying mechanisms of interactions among them, relevant signaling pathways, and ultimately their impacts on cell function in future studies will better diagnose and identify potential sites for more selective tumor therapy in human cancers. Thus, it is worthwhile to thoroughly characterize the mechanisms underlying breast cancer’s aggressive characteristics for developing novel targets such as non-coding RNAs for the diagnosis and treatment. Here, we aimed to explore the involved lncRNA-miR-mRNA network in the cell cycle and proliferation after surgery and chemo-radiotherapy in breast cancer. Using bioinformatics and experimental methods, we evaluated their prognostic values regarding patients’ clinical outcomes, such as the 5-year overall survival rate.

Materials and methods

The data collection of breast cancer datasets

The platforms used for miRs and mRNAs, including the miRs’ expression profile (GSE81000 and GSE44124) and the mRNAs’ expression profile (GSE42568 and GSE65194) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov /geo). GSEs’ data were downloaded for use with the GEOquery R package (https://bioconductor.org/ packages/GEOquery) [6]. We analyzed the selected lncRNAs, miRs, and genes with P-value < 0.05 and |LogFC| > 1 in the datasets as differentially expressed genes (DEGs), differentially expressed miRNAs (DEMs), and differentially expressed lncRNAs (DELs). Moreover, the GEPIA2 (http://gepia2.cancer-pku.cn), the cBioPortal (https://www.cbioportal.org), and the Broad Institute’s FireBrowse (http://firebrowse.org) are websites for analyzing the differential expression genes from the TCGA and Genotype-Tissue Expression projects [7]. The used platforms for miRs included the OncomiR (http:// www.oncomir.umn.edu/omcd/), miRGator 3.0 (https://tools4mirs.org), and miRCancerdb (http://mircancer.ecu.edu) databases of TCGA dataset [8]. The used databases for lncRNAs included LncRNADisease (http://www.rnanut.net/ lncrnadisease) and Lnc2Cancer v3.0 (http://bio-bigdata.hrbmu.edu.cn/lnc2cancer) of TCGA dataset [9,10]. Furthermore, we used the Tumor Immune Estimation Resource (TIMER) and Cancer Cell Lines Encyclopedia (CCLE) to determine DEGs between breast cancer and normal tissues. The DiffExp module in the TIMER database (https://cistrome.Shinyapps.io/ timer) determined differential expression patterns between normal and tumor tissues for DEGs in all TCGA tumors [11]. The CCLE database (https://portals.broadinstitute.org/ccle) explored the expression of DEGs in breast cancer cell lines [11]. Figure 1 shows a flowchart diagram for used bioinformatics analysis.

Figure 1. A flowchart diagram for used bioinformatics analysis in the present study.

The prediction of target genes and lncRNAs of candidate miRs

The target genes of miRs were identified by using the online predictive program, such as the miRmap (https://mirmap.ezlab.org/app/) [12], miRWalk2 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) [13], and Targetscan Release 7.0 (http://www.targetscan.org) [14]. Furthermore, the lncRNAs which regulate DEMs were collected by the LncRNA2target and TANRIC datasets.

The analysis of GO term pathways by the FunRich dataset

The pathway enrichment analyses of the GO database were executed through the FunRich dataset. Likewise, the same genes were explored for pathway enrichment using the g: Profiler tool (http://biit.cs.ut.ee/gprofiler) [15]. The lncRNA-miR-mRNA network was built using the Cytoscape software.

Construction and analysis of the lncRNA–miR–mRNA network

The lncRNA–miR–mRNA network was constructed and visualized using Cytoscape software based on the ceRNA theory [16]. The nodes and edges were used to represent extensive biological data. Each node represents a biological molecule, and the edges stand for the interactions between nodes and their degrees. The number of edges was calculated to exploit the hub nodes that possess essential biological functions [17]. To explore the structure and feature of the lncRNA–miR–mRNA competing triplets, a network analysis was performed using a Cytoscape plug-in called “Network Analyzer” [18].

Correlation analysis between lncRNAs, miRs, and mRNAs in breast cancer

The correlation coefficient was calculated among lncRNAs, miRs, and mRNAs. The absolute value of the correlation coefficient equal to or more than 0.5 represented a significant correlation. The interacting miRs were matched according to the miR-code database (http://www.mircode.org/) using the differentially expressed lncRNA and miRs from the above analysis. The databases of the TargetScan (http://www.targetscan. org /vert_72/), miRDB (http://mirdb.org/), and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/ index .php) were used to predict the miRs’ target genes.

Experimental sample collection

The lncRNAs, miRs, and mRNAs were evaluated in 10 healthy plasma samples and 120 breast cancer patient samples (Luminal A; ER⁺, PR⁺, Her2⁻, and ki-67 < 15%) before and after the conventional treatments. Non-inclusion criteria were: evidence of metastatic breast cancer, planning to receive any additional neoadjuvant chemotherapy, and pregnant or breastfeeding. The patients’ characteristics included age, T classification, nodal status, Her2-, hormone-receptor-positive status, and tumor subtypes. Moreover, the second, third, and fourth samples were collected one week before chemo-radiotherapy and one month after the last radiotherapy session [19]. Written informed consent was obtained from each participant before the sample collection.

Real-time PCR analysis

Total RNA was extracted from the plasma samples (250 μl) and added to 750 μl TRIzol (Beijing Tiangen Biotech Co., Ltd.). RNA extraction was carried out according to the manufacturer’s instructions. The absorbance ratio (A260/280) of total RNA, between 1.8 and 2.2, was determined using an ultraviolet (UV) spectrophotometer. The miRcute miRNA cDNA First-Strand Synthesis kit (Beijing Tiangen Biotech Co., Ltd.) for miRNAs quantification and the cDNA Synthesis Kit Manual (TAKARA BIO INC. Cat. 6 30 v.0708) for genes and lncRNAs quantification were used according to the manufacturer’s recommendations. Then, cDNA was then used in each of the real-time PCR assays with the miRcute miRNA Fluorescence Quantitative Detection kit (Tiangen Biotech Co., Ltd.) for miRs. The cycling conditions were as follows: pre-denaturation at 94°C for 2 min, followed by 40 cycles of 94°C for 20s and 60°C for 34s. The SYBR Green method (AccuPower Green Star qPCR Master Mix; Bioneer, Korea) was used for genes and lncRNAs. PCR cycling was performed as follows: one cycle at 95°C for 10 min, 40 cycles at 95°C for 20 seconds, and 60°C for 45 seconds. The melting curve analysis was run from 60°C to 95°C to confirm specific amplification [20,21]. The expression of U6 and B-actin was used to normalize miRs, lncRNAs, and genes as the Internal Reference Gene. The list of primers has shown in Table 1. The qRT-PCR reactions were performed using an ABI StepOne plus System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The expression level of the selected lncRNAs, miRs, and mRNAs were calculated using the – ΔCT method. The ΔCT was calculated by subtracting the CT values of U6 and B-actin from the CT values of the targets [22,23].

Table 1. The list of primers for Real-Time PCR.

Genes/miRNAs/LncRNAs	Forward primer	Reverse primer
MYC	TTTCTGGAGTAGCAGCTCCTAA	ATGGCCCATTACAAAGCCG
CCND1	CATGGAGGGCGGATTGGAA	CAATGACCCCGCACGATTTC
WNT1	CCAGCCTTCACTTGCTGAG	TTCAGACACGAGAGATGGAACT
PIK3CA	TGGCGGATAGACATACATTGC	ACCAGTAGGCAACCGTGAAG
AKT1	ACCGCACATCATCTCGTACA	GTGCCACCATGAAGACCTTT
E2F2	AGTCCTGAGCCTGTTTTGTG	GGGAGGGCTTTCTTTGTGTA
ERBB2	CAGGGGTGGTATTGTTCAGC	GGGAAACCTGGAACTCACCT
FGFR1	AGATGGGACCACACTTTCCATA	GGTGGCTGAAAAACGGGAAG
IGF1R	GGTGCCAGGTTATGATGATG	GCAAAGGGGACATAAACACC
APC	GACCAGAAAAGTAGCTGCCG	GCCCGGATGTGCACTAAAAT
SOCS1	TAAGGGCGAAAAAGCAGTTCC	CACGCACTTCCGCACATTC
PTEN	GGTGGGTTATGGTCTTCAAAAGG	TGGATTCGACTTAGACTTGACCT
PDCD4	CTAGCCTGCACACAATCTACAGTT	CCTGAATTAGCACTGGATACTCCT
FOXO3	CACGGCTTGCTTACTGAAGG	TCACGCACCAATTCTAACGC
B-actin	CACCATTGGCAATGAGCGGTTC	AGGTCTTTGCGGATGTCCACGT
miR-21	CGCCATGTAAAGTGCTTATAGTGC	CGATTCATTTGTTAGCGAGCGG
miR-10b	TTGGAGTTACCCTGTAGAACCG	TAAGCACGAGACTTACGGAGGA
miR-155	CGCCATGTTTAATGCTAATCGTGA	TTCCAGAAACCGATCAGAGTGT
miR-20a	CGCCATGTAAAGTGCTTATAGTGC	CGATTCATTTGTTAGCGAGCGG
miR-20b	GCCCTAAATGCCCCTTCTGGCA	ACACTGCACAGTCCCCACCATCT
miR-224	CGTTTGCCAAGTCACTAGTGGT	TTGTAAGCACGCTACATCCTGA
miR-145	GTAGGAGGTCCAGTTTTCCCAG	TGAACTTCGCAACTACCGTTTG
miR-196b	GGTCTCCCCAGTGTTCAGACTA	AGAAGGCGATTGATACGAGTCA
miR-125a	GTTGATTCTCCCTGAGACCCTTTA	GTCCTCACAACGATTCCACAAG
U6	ATGCAGTCGAGTTTCCCACAT	CCATGATCACGAAGGTGGTTT
HOTAIR	GCTTCTAAATCCGTT	CTCCACGGTAAATCCGGCAG
SRA1	CCG ACC TTA TAG ACA CGG GAC GC	TCC CAT GTT CCT GGA CTG ACG A
ZFAS1	AACCAGGCTTTGATTGAACC	ATTCCATCGCCAGTTTCT

Clinical outcomes

We evaluated the association of the expression levels of candidate lncRNAs, miRs, and mRNAs with the clinicopathological feature of the patients. The 5-year overall survival (OS) was measured by the Kaplan-Meier method.

Data analysis

The sample size was calculated based on a study by Chen et al. (2016) and according to the differences in expression of miR-21 in healthy subjects and patients with breast cancer [20]. With considering the alpha error of less than 0.05 (α) and the research power of 95% (1-β), the sample size was obtained approximately 30 in each group. G Power software version 3 (pmid:17,695,343) was used to estimate the sample size. The data analyses were performed by GraphPad Prism 7.0. Kolmogorov-Smirnov was used to evaluate the natural distribution of data. We used the t-test and the Mann-Whitney to analyze the parametric and non-parametric data in two groups. A repaid measured ANOVA to calculate the intra-group data. To evaluate the value of circulating lncRNAs, miRs, and mRNAs in predicting relapse, receiver operating characteristics (ROC) curves were constructed, and the area under the curve (AUC) was calculated. The five-year overall survival rate of candidate lncRNAs, miRs, and mRNAs was evaluated using the Kaplan-Meier method. All data were presented as mean ± SD. A p-value < 0.05 was considered to be statistically significant.

Results

Identification of DEGs

GEO and TCGA datasets were analyzed to identify differentially expressed lncRNAs, miRs, and genes in breast cancer and healthy samples. FunRich_3.1.3 software made a Venn diagram and extracted the common selected datasets (Figure 2). 156 miRs (36 up-regulation and 120 down-regulation) (Figure 2a), 2060 mRNAs (891 up-regulation and 1169 down-regulation) (Figure 2b), and 169 lncRNAs (102 up-regulation and 67 down-regulation) (Figure 2c) were obtained from the selected datasets. The top up-regulated miRs were miR-21, miR-20a, and miR-20b, and the top down-regulated miRs were miR-145, miR-196b, and miR-125a (Table 2). The target genes of selected miRs have been represented in (Tables 3, 4, 5). Likewise, the lncRNAs of selected miRs were obtained using the LncRNADisease, Lnc2Cancer v3.0, LncRNA2target, and TANRIC datasets (Table 6).

Figure 2. Venn diagram of the differently expressed lncRNAs, miRs, and mRNAbetween GEO and TCGA datasetsAllocation of (a) the 156 differently expressed miRs (36 up-regulation and 120 down-regulation), (b) 2060 differently expressed genes (891 up-regulation and 1169 down-regulation), and (c) the 169 differently expressed lncRNAs (102 up-regulation and 67 down-regulation) found between the used datasets in the present study.

Table 2. The candidate miRs involved in the cell cycle and proliferation pathways in breast cancer patients.

miRs	adj.P.Val	logFC
UP-regulated
miR-21	5.51E-06	2.84086
miR-10b	2.16e-05	2.173494
miR-155	8.78e-04	3.20482
miR-20a	4.20e-04	2.04178
miR-20b	3.60e-05	2.02111
Down-regulated
miR-145	0.00035118	−3.771778
miR-196b	0.0046717	−2.017851
miR-224	0.009643	−2.032727
miR-125a	0.00046119	−2.3917

Table 3. The predicted genes involved in the cell cycle and proliferation pathways in breast cancer patients.

Up-regulated

AKT2, BCL2,CCND1, EGFR, IGFIR, MYC,PI3KR1, IGFB1, IGFB2, IGFBR2, BCL2, CCND1, FZD9, MAPK1, MAPK9, MYC, SMAD4, TGFBR2, FGFR1, FZD9, FZD4, MAPK1, MAPK9, SMAD4, FZD3, MYC, FGF1, FGF18, PIK3R1, FGF2, WNT16, BRAF, WNT4, WNT2B, WNT9A, AKT3, FGF23, IGF1, WNT10B, IGF1R, WNT3, WNT5A, PIK3CB, CDK6, NCOA3, MAPK1, AKT2, FGF1, WNT4, WNT9A, AKT3, WNT4, PIK3R3, WNT9A, AKT3, RAF1, PIK3CB, FGF2, PIK3R1, FGF1,CDK6, CDK4, FZD10,IGF1R, FZD2, WNT3, WNT9B, FZD3,FGF8, WNT5B, KRAS, FGFR1, EGFR, WNT2, FZD3, FGF23, WNT8B, WNT1, IGF1, MAP3K2, AKT2, MAPK1, WNT7B, RAF1, WNT5A, PIK3CB, PIK3CA, PIK3R1, FGF1, CDK6, WNT16, BRAF, MYC, AKT3, WNT5A, FGFR1, KRAS, WNT10B, FGF9, IGF1R, WNT7B, FGFR1, PIK3R3, FZD3, POLK, EGFR, BRAF, AKT3, RAF1, PIK3CB, MAPK1, FGF5, FGF2, PIK3R1, POLK, CDK6, WNT2, WNT16, BRAF, KRAS, IGF1, FGF7, IGF1R, WNT9B, NCOA3, MAPK1, WNT7B, LRP6, CDK4, IGF1, RPS6KB1, RPS6KB2, PIK3R1, WNT8A, PIK3CD, FGF1, BRAF, NCOA1, WNT5A, TFAP2A, APPL1, JAG1, SHC3, JUN, SHC1, SHC2, CTNNB1, DVL3, DVL1, LEF1, ERBB2, HES5, FLT4, CDKN1A, DLL4, GADD45B, SOCS7, SOCS5, SESN3, BCL2L2, GABRB1, NUFIP2

Down-regulated

APC, PTEN, FOXO3, PDCD4, TP53, TCF7L2, E2F2, TCF7, POLK, BAK1, TP53, BRCA1, APC2, E2F1, GSK3B, BAX, SP1, RB1, GSK3B, HOXD10,SP1, AP1, TCF7, FRAT2, FRAT1, ESR1, ESR2, BAK1, DDB2,E2F3, ESRRG, CYCS, MSH2, MSH6, MSH3, ACVR1B, KIT, HEYL, GADD45A, CSNK1A1, CSNK1A1L, TNFSF11, SCH3, AXIN2, ZDHHC21, PTEN, CCND1, FOXO1, SOCS2

Table 4. The candidate genes in breast cancer patients.

Genes	adj.P-value
E2F2	1.92e-118
AKT1	2.23e-67
MYC	1.54e-47
CCND1	3.67e-44
PIK3CA	6.65e-70
IGF1R	8.67e-168
FGFR1	2.30e-77
ERBB2	5.52e-61
PDCD4	7.19e-49
WNT1	2.51e-120
APC	9.32e-154
SOCS1	9.82e-113
FOXO3	8.34e-206
PTEN	3.92e-63

Table 5. The list of candidate genes involved in cell cycle and proliferation pathways in breast cancer patient.

Genes	Target miRs
Cell cycle
APC	miR-21, miR-125a
MYC	miR-21, miR-20a, miR-125a, miR-145
CCND1	miR-21, miR-20a, miR-20b
WNT1	miR-125a, miR-145, miR-155, miR-196b, miR-224, miR-10b
SOCS1	miR-224
Proliferation
PIK3CA	miR-125a, miR-145, miR-196b, miR-224, miR-10b
AKT1	miR-21, miR-125a, miR-145, miR-196b, miR-224
PDCD4	miR-20a, miR-20b, miR-155, miR-21, miR-196b, miR-145
E2F2	miR-125b, miR-196b, miR-224
FGFR1	miR-20a, miR-20b, miR-145, miR-196b
PTEN	miR-21, miR-125a, miR-145
ERBB2 (Her2)	miR-125a,
IGF1R	miR-21, miR-155, miR-10b, miR-20a, miR-20b, miR-224, miR-145, miR-196b
FOXO3	miR-196b, miR-224, miR-155, miR-125a, miR-10b, miR-21

Table 6. Interaction analysis between selected miRs and lncRNAs involved in cell cycle and proliferation pathways in breast cancer patients.

LncRNAs	Target miRNAs
ZFAS1	miR-196b, miR-20a
HOTAIR	miR-20b, miR-20a
SRA1	miR-20b, miR-20a

We also used TIMER and CCLE databases to reveal the transcriptional expressions of the selected DEGs (Figure 3 and Figure 4). We determined the distributions of gene expression levels using the TIMER database displayed in a box plot (Figure 3). In addition, we studied the expression levels of the selected DEGs in breast cancer cell lines using the CCLE database (Figure 4).

Figure 3. The expression of the candidate genes in breast cancer(a) A box plot of the candidate genes’ transcripts in normal and breast cancer tissues in the TIMER database. Statistical significance was calculated using the Wilcoxon test, *** p < 0.00.

Figure 4. A heatmap of the genes’ expression in breast cancer cell lines (CCLE). We used mRNA expression values from the CCLE database and then ranked them. In CCLE, purple represents overexpression (top column), and blue indicates under-expression (bottom column).

Enrichment analysis of DEGs

To examine the biological function of the DEGs, GO analysis was performed in the FunRich software. The DEGs were mainly enriched in the receptor binding, protein serine/threonine kinase activity, growth factor activity, transcription factor activity, kinase regulator activity, and DNA repair protein (Table 7). Furthermore, up-regulated DEGs were exclusively enriched in the pathways, including the EGF receptor signaling and the Glypican pathway and the FGFR signaling pathway (Figure 5a). The most significantly enriched pathways associated with cancer were VEGF, MAPK, ERBB, RAP1, and RAS signaling pathways (Figure 5a). However, down-regulated DEGs were mainly enriched in the pathways, including the mTOR, P53, IGF1, and ATM (Figure 5b). Likewise, the functional enrichment of DEG genes was analyzed using the g: Profiler software (Figure 6).

Figure 5. The enrichment analysis of the GO pathway of the DEGs (a) Top 7 functional network/pathways associated with the up-regulated DEGs through GO analysis with a p-value less than 0.05. (b) Top 7 functional network/pathways related to the down-regulated DEGs through GO analysis with a p-value less than 0.05.

Figure 6. Functional enrichment by using the g: Profiler software (a) the X-axis shows the functional terms grouped and color-coded by the data source. (B, C, and F) the position of terms in the plots fixed and terms from the same branch (GO). (d) p-values in the table outputs are color-coded from yellow (insignificant) to blue (highly significant). (e) in a multi-query case, the same term is highlighted on other plots. (g) a click allows for pinning the circles to the plot together with a numeric ID, creating a more detailed result in the table below the image.

Table 7. The GO enrichment analysis for DEGs.

Source	term_name	term_id	adj_p_value	−log10(padj)
GO: CC
	transferase complex, transferring phosphorus-containing groups	GO:0061695	0.004582707	2.338877905
	extrinsic component of membrane	GO:0019898	0.008857802	2.052674039
	protein-containing complex	GO:0032991	0.025512924	1.593239769
	phosphatidylinositol 3-kinase complex	GO:0005942	0.027452467	1.561418614
GO: MF
	protein kinase activity	GO:0004672	0.000568652	3.245154
	protein kinase binding	GO:0019901	0.001324113	2.878075
	phosphotransferase activity, alcohol group as acceptor	GO:0016773	0.001477544	2.83046
	protein kinase regulator activity	GO:0019887	0.002395176	2.620663
	kinase binding	GO:0019900	0.00251691	2.599132
	kinase activity	GO:0016301	0.003284015	2.483595
	kinase regulator activity	GO:0019207	0.004281197	2.368435
	insulin receptor substrate binding	GO:0043560	0.007291041	2.13721
	transferase activity, transferring phosphorus-containing groups	GO:0016772	0.009011408	2.045207
	transmembrane receptor protein tyrosine kinase activity	GO:0004714	0.026357557	1.579095
	growth factor binding	GO:0019838	0.036785922	1.434318
	protein tyrosine kinase activity	GO:0004713	0.039149829	1.40727
	transmembrane receptor protein kinase activity	GO:0019199	0.039959246	1.398383
GO: BP
	regulation of cell population proliferation	GO:0042127	3.20372E-06	5.494345
	cell population proliferation	GO:0008283	1.56308E-05	4.806017
	cell death	GO:0008219	5.28086E-05	4.277295
	positive regulation of cell population proliferation	GO:0008284	0.000251265	3.599869
	regulation of epithelial cell proliferation	GO:0050678	0.000442664	3.353926
	positive regulation of epithelial cell proliferation	GO:0050679	0.000595434	3.225166
	epithelial cell proliferation	GO:0050673	0.000990806	3.004011
	G1/S transition of mitotic cell cycle	GO:0000082	0.001595936	2.796984
	cell cycle G1/S phase transition	GO:0044843	0.002595649	2.585754
	regulation of cell cycle	GO:0051726	0.011959955	1.92227
	reentry into mitotic cell cycle	GO:0000320	0.01257232	1.900585
	regulation of smooth muscle cell proliferation	GO:0048660	0.012635894	1.898394
	regulation of G1/S transition of mitotic cell cycle	GO:2000045	0.013206732	1.879205
	smooth muscle cell proliferation	GO:0048659	0.013499198	1.869692
	regulation of cell cycle G1/S phase transition	GO:1902806	0.022064212	1.656312
	muscle cell proliferation	GO:0033002	0.048090747	1.317938

GO: gene ontology, CC: cellular component, MF: molecular function, and BP: biological process

PPI network analysis of DEGs

PPI analysis of DEGs was performed in the FunRich software. Two significant modules with a score ≥ 7 were screened out via the FunRich software. BRAF1, PIK3R1, TGFBR2, EGFR, IGF1R, RAF1, VEGF, TWIST, and NOTCH were hub nodes with higher node degrees in module 1 (Figure 7a), and ESR1, BRCA1, HOXD11, RB1, GSK3B, SP1, TP53, and SOS1 were hub nodes in module 2 (Figure 7b). Furthermore, VEGF, TGFBR2, TWIST, HMGA1, NOTCH, and HOXD11 hub genes were selected for further analysis owing to the high degree of connectivity (Figure 7c). As described in the method section, we identified 3 co-expression lncRNAs, 9 co-expression miRs, and 14 co-expression mRNAs selected to construct the lncRNA–miR–mRNA network (Figure 8).

Figure 7. Protein-protein interaction (PPI) network constructionPPI network was constructed with the DEGs from the four datasets of GEO and TCGA. (a-b) The significant module was identified from the PPI network using the FunRich method with ≥ 7. Panel A shows module 1 of up-regulated genes. Panel B shows module 2 of down-regulated genes. Panel C shows the interaction between modules 1 and 2.

Figure 8. The lncRNA-miRNA-mRNA network The network includes 27 nodes and 126 edges. The blue rhombus, red hexagonal, and yellow oval represent the lncRNAs, miRs, and genes.

Construction of ceRNA network of lncRNA-miR-mRNA in breast cancer

To further understand the role of differentially expressed lncRNAs in breast cancer, we performed correlation analysis with the expression variability of lncRNAs, miRs, and mRNAs. Accordingly, there was a strong correlation between some network components, including miR-196b and HOTAIR, miR-20a and IGF1R gene, and miR-196b and PIK3CA gene. Moreover, we created a heat map of expression of the selected lncRNAs, miRs, and mRNAs using CIMminer (https://discover.nci.nih.gov/cimminer/home.do) (Figure 9). The results showed that three lncRNAs were involved in regulating nine miRs. After a subset of sections was obtained from nine miRNA-regulated target genes and differentially expressed mRNAs, 14 target genes regulated by miRNAs were identified. Thus, 3 lncRNAs, 9 miRs, and 14 mRNAs are directly related to lncRNA-miRNA-mRNA in breast cancer (Figure 10).

Figure 9. A plot heatmap to show the gene expression profile of DEGs in both bioinformatics (a) and experiment data (-∆CT) (B and C) in breast cancer patients The green color indicates down-regulated genes, and the red indicates up-regulated genes between tumor and normal samples.

Figure 10. ceRNA regulatory network of lncRNAs, miRs, and mRNAs in rectal cancerNote: Red lines indicate a negative correlation, and green lines indicate a positive correlation. SRA: SRA1, HOT: HOTAIR, ZFA: ZFAS1, mR-21: miR-21, mR-10: miR-10b, mR-15: miR-155, mR-20: miR-20a, mR-20b: miR-20b, mR-14: miR-145, mR-19:miR-196b, mR-1:miR-125a, mR-24: miR-224, IGF: IGF1R, E2F: E2F2, AKT: AKT1, PI3: PIK3CA, FGF: FGFR1, ERB: ERBB2, MYC: MYC, WNT: WNT1, CCN: CCND1, APC:APC, SOC: SOCS1, FOX: FOXO3, PDC: PDCD4, PTE: PTEN.

Pre-treatment analysis of the candidate lncRNAs, miRs, and genes in breast cancer patients

The expression levels of the miRs, including miR-21, miR-10b, miR-125a, miR-155, miR-20a, miR-20b, miR-224, miR-145, and miR-196b, and also the genes, including APC, SOCS1, MYC, WNT1, CCND1, IGF1R, E2F2, AKT1, PIK3CA, FGFR1, ERBB2, FOXO3, PDCD4, and PTEN, and also lncRNAs, including SRA1, HOTAIR, and ZFAS1, demonstrated significant differences in the patients compromised with the healthy subjects (P < 0.05).

Furthermore, The results showed a similar expression pattern of the selected miRs, genes, and lncRNAs between the nodal status (−)/(+), the menopause status, and the stages of patients.

The expression of tumor suppressor- and onco-lncRNAs after the conventional treatments

The expression levels of the onco-lncRNAs, including SAR1 and HOTAIR, were significantly increased in the patients compromised with the healthy subjects (P < 0.0001) (Figure 11A-b).

Figure 11. The relative expression of the selected lncRNAs in the breast cancer patientsThe relative expression levels of the lncRNAs were normalized by using a reference RNA. The oncogenic lncRNAs included: A) SRA1 and B) HOTAIR. Tumor suppressor lncRNAsincluded: C) ZFAS1. The expression level of the lncRNAs was calculated using the – ΔCT method.

The expression level of tumor suppressor lncRNA, including ZFAS1, showed a significant decrease in the patients compromised with the healthy subjects (P < 0.0001) (Figure 11C). However, the expression levels of the tumor suppressor lncRNAs and onco-lncRNAs showed substantial changes after conventional treatments compromised with pre-treatment (P < 0.0001).

Validation of selected oncomiRs and tumor suppressor miRs after the conventional treatments

The expression levels of the oncomiRs, including miR-21, miR-20a, miR-20b, miR-155, and miR-10b, were significantly increased in the patients compromised with the healthy subjects (P < 0.0001) (Figures 12A-12e). The expression levels of the oncomiRs showed a significant decrease after conventional treatments were compromised with pre-treatment. The expression levels of tumor suppressor miRs, including miR-145, miR-224, miR-125a, and miR-196b, have shown a significant decrease in breast cancer patients compared to the healthy subjects (P < 0.0001) (Figure 12F-i). However, the expression levels of the tumor suppressor miRs showed a significant increase after conventional treatments compared to the pre-treatment (P < 0.0001).

Figure 12. The relative expression of the candidate miRs in the breast cancer patientsThe relative expression levels of the miRs were normalized by using a reference RNA.The oncomiRs included: A) miR-21, B) miR-10b, C) miR-155, D) miR-20a, E) miR-20b.Tumor suppressor miRs included: F) miR-145, G) miR-196b, H) miR-125a, I) miR-224.The expression level of the miRs was calculated using the – ΔCT method.

Figure 12. (Continued).

Validation of the selected onco- and tumor suppressor genes

The expression levels of the oncogenes, including IGF1R, E2F2, AKT1, PIK3CA, FGFR1, ERBB2, MYC, and WNT1, were significantly increased in the patients compromised with the healthy subjects (P < 0.0001) (Figure 13A-h). The expression levels of the oncogenes showed a significant decrease after conventional treatments were compromised with pre-treatment. The expression levels of tumor suppressor genes, including CCND1, APC, SOCS1, FOXO3, PDCD4, and PTEN, have shown a significant decrease in breast patients compared to the healthy subjects (P < 0.0001) (Figure 13I-n). However, the expression levels of the tumor suppressor genes showed a significant increase after conventional treatments compromised with pre-treatment (P < 0.0001).

Figure 13. The relative expression of the selected genes in the breast cancer patients The relative expression level of the genes was normalized by using a reference gene. The oncogenes included: A) IGF1R, B) E2F2, C) AKT1, D) PIK3CA, E) FGFR1, F) ERBB2, J) MYC, H) WNT1. Tumor suppressor genes included: I) CCND1, G) APC, K) SOCS1, L) FOXO3, M) PDCD4 and N) PTEN. The expression level of the mRNAs was calculated using the – ΔCT method.

Figure 13. (Continued)

Figure 14. A Kaplan-Meier analysis of 5-year overall survival of candidate miRs in the breast cancer patientsThe patients with high expression of miR-21, miR-20b, and low expression of miR-125a had exhibited a shorter overall survival rate. Cutoff point: mean value of the candidate expression.A: miR-21, B: miR-20b, and C: miR-125a.

Figure 15. A Kaplan-Meier analysis of 5-year overall survival of candidate mRNAs in the breast cancer patients The patients with a high expression of IGF1R and E2F2 and a low expression of PDCD4 and PTEN had exhibited a shorter overall survival rate. Cutoff point: mean value of the selected candidate expression. A: IGF1R, B: PDCD4, C: PTEN, and D: E2F2.

The candidate miRs and mRNAs as predictive targets

The mean follow-up duration of the patients was 4.7 ± 1.3 years. During this follow-up period, six patients were deceased due to the disease. Besides, two patients had recurrences and underwent a new course of treatment. We performed the log-rank survival analysis according to the expression levels of the candidate lncRNAs, miRs, and mRNAs to investigate the predictive targets by the Kaplan-Meier analysis. Except for candidate lncRNAs, the patients with a high expression of miR-21, miR-20b, IGF1R, and E2F2 and a low expression of miR-125a, PDCD4, and PTEN exhibited a shorter overall survival rate (Figures 14–15).

Discussion

Our analysis revealed a group of differentially expressed 3 lncRNAs, 9 miRs, and 14 mRNAs in breast cancer patients. To better understand the interactions of target genes with different signaling pathways, the GO and KEGG pathway analyses were performed. The pathway analyses showed that the selected miRs could interact with more significant genes involved in signaling pathways. Our results also demonstrated a substantial difference in the expression levels of the lncRNAs, miRs, and target genes after surgery and chemo-radiotherapy compared to pre-treatment. Collectively, these results suggested that they may be a potential prognostic biomarker for patients with breast cancer.

Our results are in line with previous studies conducted by Anastasov et al. (2012) and Wang et al. (2015) [24,25]. Their findings suggest that tumor resection may alter some of the inflammatory, physiological, and even pathological features of cells in affected tissues [26]. The surgical procedure involves the removal of the tumors localized in the tissues and nearby margins. It has been demonstrated that a breast tumor’s pathological characteristics may have significant impacts on pathways associated with homeostasis, and these remarkable effects persist even more than two weeks after complete tumor resection. Indeed, alterations in cell behavior and physiological responses to tumor removal suggest that surgical treatments can play a critical role in the consequences of cell behaviors associated with pathological processes in cancer [26].

Moreover, chemotherapy drugs are systemically or locally given. Although the goal of surgery is to remove all visible tumor lesions, the cancer cells can be left after tumor resection by chemotherapy. Thus, chemotherapy can be used as adjuvant therapy to destroy the remaining cancer cells after surgery [27]. Another option available for a particular type of cancer, including breast cancer, is radiation therapy. It can directly damage the chemical structure of the genetic material of cancer cells and, therefore, be used to prevent cancer cells’ proliferation and induce their death [28,29]. In cancer research, the efficacy of therapeutic approaches is a substantial factor in continuing treatment. In cases with no detectable tumor, an increase in the oncomiRs or a decrease in tumor suppressor-miRs may indicate treatment failure. It is worth noting that the operation, chemotherapy, and radiotherapy can shrink the tumor size, leading to a more reduction in the oncomiRs and an increase in the tumor suppressor-miRs.

Many studies have shown that miRs can play a significant role in breast tumors’ cell cycle and proliferation pathways through functional interaction with the cyclin and the kinase proteins. In this setting, the up-regulation of protein kinases and down-regulation of their inhibitors by miRs can increase breast tumors’ cell viability and result in the aberration of the cell cycle alterations [30]. Likewise, the dysregulation of miRs in the cell proliferative and the cell-cycle regulatory network is attributed to resistance to cancer treatments to suppress cell growth and proliferation [29]. Liu et al. (2018) indicated that miR-26a/E2F7 feedback loop could contribute to the regulation of tamoxifen resistance in ER-positive breast cancer. Notably, it has been demonstrated that the overexpression of miR-26a could directly inhibit E2F7 expression and indirectly repress MYC expression. On the other hand, E2F7 overexpression could decrease miR-26a expression through MYC transcription inhibition. They also showed that miR-26a overexpression and E2F7 silencing could desensitize breast tumor cells to tamoxifen, reducing cell viability and inducing G1 cell cycle arrest [31].

LncRNAs and miRs and the crosstalk between them are eminent players in the pathophysiological processes of cancers. LncRNAs have been shown to play pivotal roles in regulating multiple cellular processes in normal cells or tumor formation and development of cancerous cells through various mechanisms [32]. However, few studies have described the regulatory mechanisms of lncRNA-miR-mRNA in cancer. In this setting, HOX Antisense Intergenic RNA (HOTAIR) is a lncRNA located at the HoxC gene cluster on chromosome 12. It has been demonstrated that HOTAIR has a significant role in silencing HoxD gene expression in the tumor through interactions with chromatin-modifying enzymes and regulates the expression of the PTEN gene [33–35]. In this setting, growth arrest-specific transcript 5 (GAS5) can compete with PTEN to bind miR-21; thus, it can positively regulate the expression of PTEN and negatively regulate miR-21 expression. In other words, a decreased or increased expression of GAS5 can lead to overexpression or reduced expression of miR-21 and subsequently cause the down-regulation or up-regulation of PTEN in tumor cells [36,37]. GAS5 can also induce PTEN expression by inhibiting miR-103, miR-196b, miR-205, and miR-32 [38,39]. Furthermore, the regulation of miR-181a on the PTEN gene was mediated by a lncRNA X–inactive specific transcript (XIST). The inhibition of XIST can lead to a significant increase in cell proliferation and liver cancer invasion, leading to a decreased PTEN expression [40]. These observations collectively led us to propose the novel hypothesis that such interactions may increase cell proliferation and hinder apoptosis through inhibition or direct activation of target genes involved in these pathways.

Overall, our results demonstrated that the changes in the candidate lncRNA-miR-mRNA network were associated with poor OS in breast cancer patients. In this context, the prognosis is an important indicator of disease treatment. Therefore, we hypothesized that this lncRNA-miR-mRNA network might significantly influence breast cancer’s pathogenesis and prognosis. In the present study, the patients with a high expression of miR-21, miR-20b, IGF1R, and E2F2, and a low expression of miR-125a, PDCD4, and PTEN exhibited a shorter overall survival rate. In addition, the candidate lncRNAs were not conducive to the prognosis of patients, and the regulation between miRs and lncRNAs may be more complicated. In this setting, Gao et al. (2021) observed that, unlike high expression miRs, low expression lncRNAs/miRs were significantly associated with a better overall survival rate. They found that 10 lncRNA and 6 miRs can predict the overall survival of patients with breast cancer [41]. Like Tuersong et al. (2019) [42], Jing-Jing Wang et al. (2019) found that six lncRNAs, two miRs, and nine mRNAs were negatively/positively associated with OS; thus, these RNAs may serve as potential prognostic biomarkers of OS in patients with breast cancer [43]. These instances can show the high versatile potential of the ceRNA network in the clinic. Xiao et al. (2015) demonstrated that MALAT1 might promote the expression of ZEB2 by expanding miR-200s [44]. Chen et al. (2017) revealed that the top 10 aberrantly expressed lncRNAs have essential roles through a lncRNA-miR-mRNA network [45]. Li et al. (2018) described a novel ceRNA network including 11 lncRNAs, 9 miRs, and 41 mRNAs in cancer [46]. Accordingly, an investigation into the ceRNA network may provide insight into specific disease processes and aid the development of novel therapeutic strategies [43]. Our results highlighted the potential of the lncRNA-miR-mRNA network in understanding the pathogenesis of breast cancer. Simultaneously, we demonstrated that this panel could dynamically change during the treatment phases. They may propose innovative concepts and an experimental basis for identifying predictive biomarkers and therapeutic targets.

Conclusion

Overall, monitoring the lncRNA, miRs, and target genes during the treatment can be an appropriate diagnostic tool for estimating the identification and efficacy of therapeutic approaches. Understanding their underlying interactions will better diagnose and identify potential sites for more selective tumor therapy. Thus, it is worthwhile to thoroughly characterize the mechanisms underlying breast cancer’s aggressive characteristics for developing novel targets such as non-coding RNAs for the diagnosis and treatment.

Author contributions

MM and SK: sample processing, data analysis, and manuscript preparation. AMA: study conception and design and manuscript revision. VK, MGZ, MZ, and RO: sample collection and data analysis. EP: article revision. HR and ZM: sample processing and data analysis. ZB: article revision.

Compliance with ethical standards

The data that support the findings of this study are available from the corresponding author, AMA, upon reasonable request.

Availability of data and materials

Please contact the authors for data requests.

Ethical approval

All methods were performed under the relevant guidelines and regulations. All procedures performed in studies involving human participants were under the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. The experimental procedures and care protocols were approved by a review board committee of Dzful University of Medical Sciences (NO: IR.DUMS.REC. 1396.30) and registered by the Iranian Randomized Control Trial (IRCT) ethical board (NO: IRCT2016080818745N11). Written informed consent was obtained from each participant before the sample collection.

Acknowledgments

This study was co-funded by the Dzful University of Medical Sciences (Grant Number: 933) and the Tehran University of Medical Sciences (Grant Number: 38472). Research reported in this publication was also supported by Elite Researcher Grant Committee under award number [Grant Number: 971158] from the National Institute for Medical Research Development (NIMAD), Tehran, Iran.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

None of the funding sources had any role in the study design, the collection, analysis, and interpretation of data, or the decision to submit the article for publication;the Tehran University of Medical Sciences [38472];