Potential Diagnostic and Prognostic Utility of miR-141, miR-181b1, and miR-23b in Breast Cancer

miRNAs, a group of short noncoding RNAs, are key regulators of fundamental cellular processes and signaling pathways. Dysregulation of miRNA expression with known oncogenic or tumor suppressor functions has been associated with neoplastic transformation. Numerous studies have reported dysregulation of miRNA-141, miR-181b1, and miR-23b in a wide range of malignancies, including breast cancer. To the best of our knowledge, no previous study had demonstrated the expression of miR-141-3p, miR-181b1-5p, and miR-23b-3p in different histological grades and molecular subtypes of breast cancer. Here, we identified differential expression of these three miRNAs in breast cancer tissues compared with benign breast fibroadenomas. In addition, high expression levels of miR-141-3p and miR-181b1-5p are strongly associated with aggressive breast carcinomas. We also confirmed the clinical potential of using the three miRNAs individually or combined as diagnostic and prognostic markers in breast cancer. Using bioinformatics analyses, we identified 23 hub genes of these three miRNAs which are involved in key signaling pathways in breast cancer. Furthermore, the KM plotter online database analysis demonstrates the association between elevated expression of miR-141 and miR-181b and shorter overall survival of breast cancer patients. Together, our data suggest an oncogenic role of the studied miRNAs and highlight their molecular roles and potential clinical applications in breast cancer.


Introduction
Breast cancers, which are heterogeneous and complex in nature, have been classified into different categories with unique genetic, epigenetic, and molecular characteristics [1]. Despite the recent advances in diagnosis and therapeutic approaches of breast cancer, some subsets remain deadly due to a high incidence of relapse and metastases. Benign breast diseases are non-cancerous disorders that encompass distinct histologic entities-non-proliferative, proliferative without atypia, and atypical hyperplasia. Although benign breast diseases have been recognizable as distinct from breast cancers, the association of these benign conditions to the development of breast cancer has been a controversial topic [2]. Therefore, identifying new biomarkers with diagnostic and prognostic values in breast cancer is urgently needed. Recently, many reports have demonstrated specific patterns of micro-RNA (miRNA) expression as potentially useful therapeutic targets, as well as diagnostic, and prognostic biomarkers in breast cancer [3][4][5]. 2 of 17 miRNAs are a class of small, noncoding RNAs that vary in size from 19 to 25 nucleotides. They posttranscriptionally control gene expression by regulating mRNA through binding to complementary sequences with subsequent mRNA degradation and/or translational repression [6]. Dysregulation of miRNAs is often associated with cancer-related signaling pathways by functioning as oncogenes (oncomiRs) or tumor suppressor miRNAs [7][8][9]. To date, a growing number of studies have shown the association of aberrant miRNA expression with various types of cancer such as leukemia, prostate, lung, and breast cancer [10][11][12]. Among miRNAs that are implicated with breast tumorigenesis, miR-141, miR-23b, and miR-181b are critical regulators [13]. Functionally, overexpression of these three miRNAs significantly induces tumor growth, proliferation, and metastasis in breast cancer cell lines [14][15][16]. Despite many studies on mechanisms of miR-141, miR-181b, and miR-23b in breast cancer pathogenesis were carried out, their accurate expression in tissues of breast cancers has not yet been elucidated. To the best of our knowledge, no previous study has demonstrated the expression of miR-141, miR-23b, and miR-181b1 in different histological grades and molecular subtypes of breast cancer.
The aim of the current study was to assess the expression pattern and clinical significance of these three miRNAs in breast cancer tissues compared to benign breast fibroadenomas. Here, our results show that miR-141-3p, miR-181b1-5p, and miR-23b-3p are differentially expressed in breast cancer tissues compared with benign breast fibroadenomas. We demonstrated that high expression levels of miR-141-3p and miR-181b1-5p are strongly associated with highly aggressive breast carcinomas: grade III and triple-negative molecular subtype. Furthermore, we confirmed that using the three miRNAs individually or combined can serve as potential diagnostic and prognostic markers in breast cancer. Our data highlight the functional and clinical significance of these three miRNAs in breast cancer biology through bioinformatics analyses of their target genes and networks.
To further explore the differential expression of these three miRNAs in breast cancer, we assessed their expression levels in breast cancer tissues of different histological grades in comparison to benign breast fibroadenomas. A significant upregulation of the three miRNAs in grades II and III breast cancer was detected when compared with benign breast tumors ( Figure 1B). It is noteworthy that the levels of miR-141-3p, miR-181b1-5p, and miR-23b-3p expression were significantly upregulated (p < 0.0001) with 19.3, 10.2, and 6.17 folds, respectively, in grade III breast cancer compared with benign breast tissues.
Next, we studied the expression levels of these three miRNAs in different molecular subtypes of breast cancer compared to benign breast fibroadenomas. We defined the molecular subtypes using ER, PR, HER2, and Ki-67 as immunohistochemical surrogate markers as described by the St. Gallen's Consensus [17]. Our results demonstrate significant upregulation of miR-141-3p and miR-181b1-5p in all molecular subtypes of breast cancer when compared with benign breast tumors ( Figure 1C). For miR-23b-3p, although significant differences were detected between triple-negative breast cancer (TNBC), luminal A, and luminal B when compared with benign breast tumors, we did not find a significant difference with the HER2-positive molecular subtype. Furthermore, obvious high fold changes of miR-141-3p, miR-181b1-5p, and miR-23b-3p expression (39.1, 20.6, and 10.6 folds, respectively) were detected in TNBC in comparison with benign breast tissues.

Upregulation of miR-141-3p and miR-181b1-5p Expression Levels Is Strongly Associated with Highly Aggressive Breast Carcinomas
The high expression patterns of miR-141-3p, miR-181b1-5p, and miR-23b-3p in grade III and TNBC prompted us to pursue assessment of the three miRNAs in these two groups in comparison to other subtypes of breast cancer. The expression of miR-141-3p, miR-181b1-5p, and miR-23b-3p was then examined in breast cancer of various histological grades (grade II versus III). Significant high expression (p < 0.05) of miR-141-3p and miR-181b1-5p was associated with high-grade breast tumors (grade III) when compared with grade II breast cancer ( Figure 1B). Of note, no significant difference was detected in the expression level of miR-23b-3p when comparing breast tumors from different histological grades.
We next assessed the expression of the three miRNAs in TNBC tissues in comparison to other breast cancer molecular subtypes. Significant high expression (p < 0.05) of miR-141-3p and miR-181b1-5p in TNBC was detected compared with HER2-positive, luminal A, and luminal B subtypes ( Figure 1C). However, no significant difference could be identified in the expression level of miR-23b-3p when comparing TNBC with other molecular subtypes of breast cancer. These results suggested the prognostic ability of these miRNAs in discriminating grade III from grade II and TNBC from other molecular subtypes of breast cancer.

Prediction of miR-141-3p, miR-181b1-5p, and miR-23b-3p Target Genes in Breast Cancer
To predict the putative target genes of the three miRNAs, we used the miRDB online microRNAs target predication database [18]. All genes with a target score of ≥90 have been considered as differentially expressed genes (DEGs) [19]. Our analysis showed that 827 DEGs are target genes for miR-141-3p, miR-181b1-5p, and miR-23b-3p. In order to pinpoint genes associated with breast cancer from the identified DEGs, we employed the Pathway Studio ® web [20], an exhaustive resource of biological data to help better exploration of molecular interactions between molecules, cellular processes, and various diseases. Our analysis revealed that 310 target genes for these three miRNAs are linked to breast cancer (Supplementary Table S1).

Gene Ontology (GO) Functional Analysis and Pathway Enrichment Analysis for DEGs
To further understand the function and the biological mechanisms of the identified DEGs for miR-141-3p, miR-181b1-5p, and miR-23b-3p that are associated with breast cancer, we performed GO functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the freely accessible Database for Annotation, Visualization, and Integrated Discovery (DAVID) website [21,22]. The GO analysis classified the target genes of the studied miRNAs into three categories: biological process (BP), cellular component (CC), and molecular function (MF) [23]. Annotation results show that the identified DEGs are primarily enriched in the regulation of the phosphate metabolic process, regulation of phosphorylation, and regulation of cell motion. These biological processes are crucial for regulation of cancer cell metabolism, shape, and cell survival under stress. The five most significant (p < 0.05) enrichment terms of each GO category and the KEGG pathway enrichment analysis are listed in Table 4. Our findings from the KEGG pathway enrichment analyses reveal that miR-141-3p, miR-181b1-5p, and miR-23b-3p are mainly enriched in pathways in cancer, small-cell lung cancer, the ErbB-epidermal growth factor receptor (EGFR) signaling pathway, chronic myeloid leukemia, and mitogen-activated protein kinase (MAPK) signaling pathways. The complete GO analysis and KEGG pathway enrichment analysis for the identified DEGs are provided in Supplementary Table S2. These findings demonstrate that the DEGs, regulated by these three miRNAs, are involved in multiple biological processes and signaling pathways that may lead to invasiveness and metastasis of breast cancer cells.

Protein-Protein Interaction (PPI) Network Analysis, Identification of the Hub Genes, and Module Analysis
The identified 310 DEGs were then uploaded to the STRING online database to construct the PPI network [24]. Only PPI pairs with an association score of >0.4 were included in the created PPI network. The output data of the constructed PPI network (as .tsv) were afterwards imported to Cytoscape for further network analysis and identification of hub genes [25]. Our analysis using the Cytoscape plugin Molecular Complex Detection (MCODE) revealed 11 modules (clusters). Using the connectivity degree of ≥5 [26], we identified 23 hub genes in the top two significant modules: 10 in the first module and 13 in the second module (Supplementary Table S3).
To explore the PPI network of the identified hub genes, all genes in the two candidate modules were uploaded separately to the STRING database. Hub genes of module 1 form a single cluster of highly interconnected PPIs (Figure 3A), while genes of module 2 form two subclusters which are interconnected ( Figure 3B). Each network node represents proteins produced by a single, protein-coding gene locus, while edges represent specific and meaningful protein-protein associations that contribute to a shared function. Furthermore, GO annotation and pathway enrichment analysis show the BP, CC, MF, and all the KEGG pathways in which the hub genes of each module were associated (Supplementary Table S3). The three most significant enrichment terms from each category of GO annotation and KEGG pathway enrichment analysis are shown in Figure 3C,D. This analysis revealed that the hub genes of module 1 are significantly involved in the regulation of protein ubiquitination such as ubiquitin-protein transferase activity, ubiquitin-protein ligase activity, and ubiquitin conjugating enzyme activity. Additionally, five genes of this module are highly enriched (FDR = 3.70 × 10 −8 ) in the ubiquitin-mediated proteolysis pathway. On the other hand, hub genes in module 2 are highly involved in bioactive lipid receptor activity, lysophosphatidic acid receptor activity, and G-protein-coupled receptor binding. This module was also found to be significantly enriched with genes of the neuroactive ligand-receptor interaction pathway, pathways in cancer, and Ras-related protein 1 (Rap1) signaling pathways.   We next examined the prognostic value and clinical significance of miR-141, miR-181b, and miR-23b expression in breast cancer using the Kaplan-Meier (KM) plotter miRNA breast cancer online database [27]. Among datasets included in the KM plotter, we employed the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort dataset, which includes data from 1262 patients, to assess the association between the three studied miRNAs and OS in breast cancer patients [28]. In agreement with our results, for all 1262 patients, the OS is significantly lower in patients with high expression of both miR-141 (HR = 1.43, 95% CI = 1.17-1.74, p = 0.00037) and miR-181b (HR = 1.47, 95% CI = 1.19-1.82, p = 0.00029) compared with patients with low expressions ( Figure 4A,B). No statistically significant difference was observed in the OS time between breast cancer patients with low and high expression of miR-23b (HR = 1.14, 95% CI = 0.94-1.39, p = 0.18). Hence, the KM plotter results emphasized the prognostic significance of miR-141 and miR-181b in breast cancer.

Discussion
Recent biomarker research advances have demonstrated the great potential for developing miRNAs as novel biomarkers and therapeutic targets in breast cancer. In the current study, we demonstrated differential expression of miR-141-3p, miR-181b1-5p, and miR-23b-3p in breast cancer tissues compared with benign breast tumors using qRT-PCR. Interestingly, high fold changes of miR-141-3p and miR-181b1-5p were detected in tumors with known aggressive behavior-grade III and triple-negative breast tumors-when compared with benign control. Furthermore, our results demonstrated the potential diagnostic and prognostic utility of these three miRNAs individually or combined in breast cancer. Lastly, we explored the biological functions and clinical significance of these three miRNAs in breast cancer using bioinformatics analyses.
Together with other previously reported data, our findings suggest the oncogenic role of the three miRNAs, as their high expression is detected in breast cancer tissues in comparison to benign breast fibroadenomas. Our results are consistent with both in vitro and in vivo studies that

Discussion
Recent biomarker research advances have demonstrated the great potential for developing miRNAs as novel biomarkers and therapeutic targets in breast cancer. In the current study, we demonstrated differential expression of miR-141-3p, miR-181b1-5p, and miR-23b-3p in breast cancer tissues compared with benign breast tumors using qRT-PCR. Interestingly, high fold changes of miR-141-3p and miR-181b1-5p were detected in tumors with known aggressive behavior-grade III and triple-negative breast tumors-when compared with benign control. Furthermore, our results demonstrated the potential diagnostic and prognostic utility of these three miRNAs individually or combined in breast cancer. Lastly, we explored the biological functions and clinical significance of these three miRNAs in breast cancer using bioinformatics analyses.
Together with other previously reported data, our findings suggest the oncogenic role of the three miRNAs, as their high expression is detected in breast cancer tissues in comparison to benign breast fibroadenomas. Our results are consistent with both in vitro and in vivo studies that demonstrated the upregulation of miR-141, miR-181b1, and miR-23b expression in highly metastatic breast cancer cell lines [29][30][31][32]. Upregulation of miR-141 has been found to promote breast tumorigenesis through regulation of the phosphatidylinositol-4,5-bisphosphate 3-kinase/Protein kinase B (PI3K/AKT) signaling pathway, as well as the secretion of certain cytokines and growth factors [15]. Furthermore, it has been demonstrated that miR-181b overexpression has dual effects in inducing breast tumorigenesis and aggressiveness behavior, as it acts by (a) suppressing the expression of the proapoptotic Bim signal, causing dysregulation of the cell cycle, overgrowth, and the promotion of tumorigenesis [33,34], and (b) activating the most common signaling pathways of breast tumorigenesis, including IL6/STAT3 [35], transforming growth factor-β (TGF-β) [36], HIF-1 [34], and WNT/β-catenin [37]. Additionally, Jin et al. found that miR-23b is regulated by Her2/neu receptors and it can target tumor suppressor genes, leading to breast cancer initiation and progression [32].
Many previous studies are at odds with our results, as they demonstrated a tumor-suppressive role of miR-141, miR-181b, and miR-23 in different types of cancer such as lung, colon, prostate, and breast cancers [38][39][40][41]. For instance, other miRNA expression profiling studies reported downregulation of miR-141 and miR-23b in breast cancer tissues compared with matched surrounding tissues [39,42,43]. It was demonstrated by Wu et al. that miR-141 inhibits the growth and motility of hepatocellular carcinoma cells by targeting the Zinc finger E-box Binding homeobox 2 (ZEB2) gene [44]. Another study by Chen et al. reported the tumor suppressor role of miR-141 in renal cell carcinoma by inhibiting cell proliferation, migration, as well as invasion through targeting erythropoietin-producing hepatocellular A2 (EphA2) via modulating the EphA2/p-FAK/p-AKT/MMPs signaling cascade [45]. Similarly, miR-181b exerts its tumor suppressor effect in many cancers by promoting apoptosis and inhibiting cancer cell proliferation, migration, and invasion. These tumor suppressor effects of miR-181b are exerted by suppressing high-mobility group box-1 (HMGB1) and downregulation of NOVA alternative splicing regulator 1 [46,47]. Likewise, miR-23b was found to act as a tumor suppressor in bladder carcinoma by inhibiting tumor cell proliferation, colony formation, and migration and causing cell cycle arrest [48]. This discrepancy in miRNA expression in different studies can be explained by using different types of tissues and cell lines, different sample size, short-and long-term miRNA measurement together with employing different platforms of miRNA profiling techniques. Future studies with preclinical models will help solve this presumed contradiction.
To obtain further insights into the biological roles of miR-141-3p, miR-181b1-5p, and miR-23b-3p in breast cancer, bioinformatics and functional analyses were conducted using multiple bioinformatics platforms. Our KEGG enrichment analysis of the identified DEGs revealed that miR-141-3p, miR-181b1-5p, and miR-23b-3p regulated genes are involved in metabolic, biological, and cellular pathways in breast cancer. The most enriched pathways of these DEGs are cancer-related signaling pathways such as ErbB-epidermal growth factor receptor (EGFR) and mitogen-activated protein kinase (MAPK) signaling pathways. This is consistent with previous studies which demonstrated that miR-141 is an important modulator of the EGFR signaling pathway as it targets EGFR, preventing its translation in breast cancer cells [49,50]. Furthermore, Chiyomaru et al. reported that the miR-23b/27b cluster is involved in enhanced breast cancer cell proliferation, migration, and invasion through regulation of EGFR and c-Met signaling pathways [51]. In agreement with our KEGG pathway analysis, it was reported that dysregulation of MAPK leads to the loss of ER expression and plays a role in the development of the ER-negative breast cancer phenotype [52].
After a series of bioinformatic analyses, a total of 23 genes with a connectivity degree of ≥5 were considered as candidate hub genes that are regulated by miR-141-3p, miR-181b1-5p, and miR-23b-3p. The KEGG enrichment analysis of the identified 23 hub genes detected their involvement in critical cellular processes. The hub genes in module 1 are significantly involved in the ubiquitin-mediated proteolysis pathway, which has been previously proposed as a therapeutic target in different types of cancer, including breast cancer [53,54]. Interestingly, distinct breast cancer studies demonstrated that a number of fundamental proteins that act as oncogenes or tumor suppressor genes are significant participants in the ubiquitin-proteasome pathway [53,55]. Among the hub genes of the first module is Casitas B-lineage lymphoma-B (CBLB), a multifunctional adaptor protein and E3 ubiquitin ligase, which plays a role in regulating cell signaling. According to a previous in vitro study by Kang et al., the CBLB protein family promotes breast cancer development via inhibiting the tumor suppressor activity of TGF-β [56]. The hub genes in module 2 were significantly associated with neuroactive ligand-receptor interaction, pathways in cancer, and Ras-proximate-1 or Ras-related protein 1 (Rap1) signaling pathways. The neuroactive ligand-receptor interaction pathway involves many G-protein-coupled receptors and is one of the most common dysregulated pathways of malignancies, including breast cancer [57]. Rap1 is a small G protein that regulates diverse processes such as cell adhesion, cell-cell junction formation, and cell polarity, which are critical for cellular migration, invasion, and metastasis [58,59]. Moreover, the involvement of upregulated Rap1 in breast tumor formation and progression to malignancy has been reported in many studies [60,61]. Among the hub genes of module 2 is cannabinoid receptor-1 (CNR1), a component of neuroactive ligand-receptor interaction and Rap1 signaling pathways. Previous studies with human breast cancer cell lines demonstrated a positive correlation between CNR1 expression and the invasiveness of breast cancer [62]. The elevated expressions of Lysophosphatidic Acid Receptor 1 and 3 (LPAR1 and LPAR3), two other hub genes in module 2, are associated with poorly differentiated and advanced stages of breast cancer [63,64]. Taken together, these findings about the KEGG pathways and hub genes targeted by the three miRNAs provide further insight into the underlying biological processes and molecular mechanisms involved in initiation and progression of breast cancer. However, further analysis and validation of the molecular mechanisms are needed to confirm the clinical relevance of our findings.
Finally, the KM plotter was employed to evaluate the potential prognostic value of our studied miRNAs among 1262 patients with breast cancer. Our results confirmed that patients with elevated levels of miR-141 and miR-181b had shorter OS compared with patients with a lower expression of these miRNAs. These data strongly support the prognostic value of miR-141 and miR-181b in breast cancer. In alignment to our survival bioinformatics analysis, previous studies reported that high expression of miR-141 and miR-181 was markedly associated with shorter survival time of patients with breast cancer [65,66]. The worse OS associated with miR-141 and miR-181 can be explained by their suggested oncogenic role. Conversely, our results are completely at odds with Ping et al., who demonstrated that high expression of miR-141 was associated with better OS in TNBC patients [67]. One intrinsic limitation of our study is the small sample size, which did not allow us to study the associations between the three studied miRNAs and the clinical characteristics of the patients. In future research, it will be noteworthy to follow up on our patients' survival, relapse, and metastasis to independently confirm the findings of our bioinformatics analysis. Much work remains to be done in order to further explore and better understand the role of these three miRNAs in human breast cancers and their clinical significance.

Tissue Samples and Ethical Approval
A total of seventy breast cancer tissue specimens and thirty mammary gland fibroadenomas were obtained from female patients undergoing a surgical procedure at Surgical Department of Kasr Alainy Teaching Hospital, Faculty of Medicine, Cairo University between March 2018 and September 2019. Written consent was obtained from all patients prior to tissue samples collection. Both malignant and benign tissues were sent for histopathological processing, diagnosis, and grading in the Pathology Department, Faculty of Medicine, Cairo University. Formalin-fixed paraffin-embedded (FFPE) blocks were prepared and stored in an environment with controlled humidity and temperature. Histological grades of breast cancer samples were evaluated using the modified Bloom and Richardson criteria [68]. ER, PR, HER2, and Ki-67 statuses were assessed using immunohistochemical staining and scored according to the American Society of Clinical Oncology/College of American Pathologists (ASCO/CAP) guidelines [69]. Molecular subtype classification was carried out as previously described [70,71]. Breast cancer patients with the following characteristics were excluded: (i) received preoperative neoadjuvant chemotherapy, immunotherapy, or radiation therapy; (ii) diagnosed with inflammatory breast cancer; (iii) had a prior history of cancer. All procedures were performed in accordance with the ethical guidelines of the Helsinki Declaration after approval from the research ethics committee at Cairo University, Cairo, Egypt (IRB# BC 2136).

RNA Isolation and Quantitative RT-PCR
For RNA extraction, 3-5 sections of 10 µm thick slices were prepared from each FFPE breast tissue sample. Total RNA was extracted from the FFPE breast tumor samples, after removing paraffin with xylene using a RNeasy Fibrous Tissue kit (Qiagen, cat. No. 74704, Valencia, CA, USA) as described in the manufacturer's protocol. The purity and concentration of RNA were evaluated using a NanoDrop ND-100 Spectrophotometer (Thermo Scientific, Waltham, MA, USA). Then, we converted total RNA to complementary DNA (cDNA) using the miScript II RT kit (Qiagen, cat. No. 218161, Valencia, CA, USA) as described in the manufacturer's protocol. qRT-PCR was performed using the SYBR Green PCR Super Mix (Qiagen, cat. No. 218300, Valencia, CA, USA) and the specific microRNA Locked Nucleic Acid PCR primers (Qiagen, cat. No. 218300, Valencia, CA, USA). Specific miScript forward primers were as follows: mir-141-3p: 5 -UAACACUGUCUGGUAAAGA UGG-3 , mir-181b1-5p: 5 -AACAUUCAUUGCUGU CGGUGGG U-3 , miR-23b-3p: 5 -AUCACAUUGCCAGG GAUUACC-3 ; the reverse primers were universal. We used the BIO-RAD CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) to perform qRT-PCR. The PCR parameters were as follows: 95 • C for 15 min, followed by 40 cycles of 94 • C for 15 s, 55 • C for 30 s, and 70 • C for 30 s. Each sample was examined in triplicate and the expression of miRNAs was defined based on the threshold cycle (Ct); relative expression levels were normalized using the ∆∆C t equation with respect to RNU6B internal control [72].

Target Genes Prediction in Breast Cancer
To predict target genes of miR-141-3p, miR-181b1-5p, and miR-23b-3p, we used the miRDB (http://mirdb.org/) online microRNA target prediction database on 7 August 2020 [18]. We selected the genes with a target score of ≥90 for further filtration and network analysis [19]. In order to identify DEGs involved in breast cancer, the gene list with a target score of ≥90 generated from miRDB was imported to the Pathway Studio Web (Elsevier, Netherlands) (http://www.pathwaystudio.com) [20]. This software is an exhaustive resource of biological data extracted from published scientific research to help better exploration of molecular interactions between molecules, cellular processes, and various diseases. In our study, pathway studio analyses were conducted to find molecular interactions between the DEGs detected from miRDB and breast cancer.

GO Functional Analysis and Pathway Enrichment Analysis of Breast Cancer DEGs
To understand the potential biological functions and molecular mechanisms of the detected DEGs, we performed a GO functional analysis and a KEGG pathway enrichment analysis using the freely accessible DAVID database (http://david.ncifcrf.gov/) (version 6.8) [21,22]. GO functional analysis included assessing the BP, CC, and MF that might be affected by miR-141-3p, miR-181b-5p, and miR-23b-3p [23]. In addition, enrichment analysis was used to identify KEGG pathways associated with identified DEGs in breast cancer. Statistically significant difference was considered for p < 0.05.

PPI Network Construction, Identification of the Hub Genes, and Module Analysis
We submitted the retrieved DEGs that are targets of miR-141-3p, miR-181b1-5p, and miR-23b-3p to the STRING database (version 11.0) (https://string-db.org/) on 9 August 2020 to construct the PPI network [24]. STRING is a database that allows retrieval of known and predicted protein-protein interactions from a variety of sources such as genomic context predictions, high-throughput lab experiments, co-expression, automated text mining, and previous knowledge in databases [24]. The aim of PPI network is to summarize the predicted association for a group of proteins, where the nodes represent the protein and the edges represent the predicted functional associations. We used a medium confidence score of >0.4 as a threshold to construct the PPI network, where only interactions above this score were included in the constructed network. The output data for the PPI interaction network were next exported as tab separated values (.tsv format) and then imported to Cytoscape (version 3.8.0), an open source software platform, for further network analysis [25]. In Cytoscape, we applied the clustering algorithm MCODE to explore the gene-gene connectivity through focusing in highly interconnected nodes that often possess specific biological functions known as clusters or modules. The criteria were set as follows: MCODE score cutoff = 0.2, k-core = 2, max. depth from seed = 100, and degree cutoff = 2. We considered genes with a connectivity degree (MCODE score) of ≥5 as hub genes [26]. Subsequently, GO functional analysis and KEGG pathway enrichment analysis of the identified hub genes were carried out using the STRING database; p < 0.05 was considered to be significantly different.

OS Analysis of miR-141, miR-181b, and miR-23b in Breast Cancer Patients
We further used the KM plotter [27], an integrated online bioinformatics tool, to assess the prognostic significance of miR-141, miR-181b, and miR-23b in patients with breast cancer. In the current analysis, we used the METABRIC cohort online dataset via the KM plotter miRNA breast cancer online database, which is characterized as follows: it includes data from 1262 patients, long follow-up (median: 94 months), average characteristics (78% ER positive, 12% HER2-positive), and available treatment data [28]. To analyze the prognostic value of the three miRNAs, the patient samples were split into two groups according to the auto select best cutoff feature of each proposed miRNA, where the median survival was reported in months. The two patient cohorts were compared by a Kaplan-Meier survival plot, and the hazard ratio with 95% confidence intervals and log rank p-value were calculated. A p-value of <0.05 was considered statistically significant. The entire analysis of the KM plotter was performed using R 2.15.0 statistical software.

Statistical Analysis
At least three independent measurements of miRNAs were conducted. Our data analysis was performed using R version 3.6.1. All tests were done after the normality was tested with the Shapiro-Wilk normality test. Levels of miRNA expression in histological grades and molecular subtypes of breast cancer were described as median fold change in comparison with benign cases by the Kruskal-Wallis rank sum test and the Dunn (1 964) Kruskal-Wallis multiple comparison test. Levels of miR-141-3p, miR-181b1-5p, and miR-23b-3p expression were measured using qRT-PCR in the breast tissue specimens and normalized to RNU6B mRNA as the reference gene. The fold change was calculated as the ratio of miRNA expression to the internal control. Of note, log2 values were used to visualize the gene expression data for easier data interpretation. The correlations between the studied miRNA expressions were evaluated using Spearman's Rho correlation coefficient. To assess the diagnostic and prognostic accuracy of miR-141-3p, miR-181b1-5p, and miR-23b-3p, ROC curve analysis was conducted and the AUC was calculated [73]. p < 0.05 was considered to indicate a statistically significant difference.

Conclusions
To sum up, our results clearly emphasize the diagnostic and prognostic roles of miR-141-3p, miR-181b1-5p, and miR-23b-3p individually or combined in breast cancer. Upregulation of the three miRNA expressions was detected in breast cancer tissues compared with benign breast fibroadenomas. High expression levels of miR-141-3p and miR-181b1-5p are strongly associated with highly aggressive breast carcinomas. Furthermore, KM plotter online database analysis confirmed the association of elevated levels of miR-141 and miR-181b with shorter survival time of breast cancer patients. By using different tools of bioinformatics and target prediction analyses, we identified many potential miRNA regulatory pathways, which may improve our understanding of the molecular and biological mechanisms of breast cancer. Further work is needed using a larger cohort of patients in order to validate the potential clinical uses of these three miRNAs as biomarkers in breast cancer.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.