miR-130b Acts as an Oncogene and Prognostic Biomarker of Breast Cancer: A Bioinformatic Analysis

Background: Breast cancer (BC) is the most prevalent cancer among females globally. microRNAs (miRNAs) could regulate the expression levels of cancer-related genes through binding with target mRNAs. In various cancers, the abnormal expression of miR-130b has been detected. We aims to investigate the molecular mechanism and biological function of miR130b in breast cancer. Methods: We obtained two microRNA expression proles from the Gene Expression Omnibus (GEO) database, including GSE45666 and GSE26659. We identied differentially expressed miRNAs (DE-miRNAs) between BC tissue and normal breast tissue based on the GEO2R web tool. DE-miRNAs were ltered by signicant prognostic value resulting from Kaplan–Meier plotter. We used the JASPAR database to explore upstream regulators of miR-130b. The potential molecular mechanisms of miR-130b correlation genes were revealed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in WebGestalt. Protein–protein interaction (PPI) network of miR130b target genes was constructed by STRING. Cytoscape software was used to visualize the PPI network and hub genes. Results: miR-130b was highly expressed in breast cancer tissues, which positively correlates with poor prognostic. JASPAR revealed THAP11 might be the upstream regulator of miR-130b. In addition, GO, and KEGG pathway revealed that miR-130b positively regulated PFKP, STAT1, SRC, and NOTCH2, participating in the Thyroid hormone signaling pathway. The PPI network further identied that AR, KIT, and ESR1 as hub genes in BC development. Conclusion: miR-130b, which is regulated by THAP11, acts as an oncogene and prognostic biomarker in BC by mediating the Thyroid hormone signaling pathway and potential target genes. miR-130b might be a novel therapeutic target for BC treatment.


Introduction
A recent study showed that 276480 American women were diagnosed with breast cancer (BC) in 2020, ranked the most prevalent cancer in female cancers, and caused the death of 42170 women ). The incidence and death rate of BC increase exponentially with aging ). Despite recent progress in prevention and treatment, still many people died from BC each year (Kan et al. 2018). Therefore, identifying speci c molecules that can serve as tumor biomarkers for detection and molecular pathways to provide potential therapeutic targets is critical for BC management. 2019) suggested that TH regulates BC's metastasis by controlling the signaling pathway, which induced integrin αvβ3 to FAK/paxillin/cortactin/N-WASP/Arp2/3 that converge in cell motility. As stated above, it has been proven that the Thyroid hormone signaling pathway is associated with BC cell. However, the relationship between miRNAs and the TH signaling pathway, which promotes BC's evolution and expansion, is still unclear.
This study demonstrated that miR-130b was upregulated in BC tissues by analyzing two microRNA expression pro les in Gene Expression Omnibus (GEO). To get a deep insight into underlying mechanisms contributing to BC, we evaluate the relationship between miR-130b expression and prognostic value. Furthermore, we explored the correlation genes. THAP11, as a transcription factor, was identi ed to be an upstream regulator of miR-130b. The results of GO and KEGG pathway analysis demonstrated that miR-130b might regulate the TH signaling pathway in BC. In addition, a Proteinprotein interaction (PPI) network was constructed to investigate hub genes. These results provide a novel therapeutic target for BC and ameliorating the prognosis.

Survival analysis
The Kaplan-Meier plot (http://kmplot.com/analysis/index.php?p=background) was used as a web tool, containing survival data of 21 cancers based on GEO, EGA, and TCGA databases. That investigate and verify the prognostic value of 54k genes (mRNA, miRNA, protein). miRNA BC database was applied to evaluate the prognostic values of DE-miRNAs. Numerous quantile expressions of DE-miRNAs were used to divide the BC patients into two groups. DE-miRNAs with P-value<0.05 were considered to be statistically signi cant.

Prediction of upstream regulators
We investigated the relationship between miR-130b and genes obtained from the TCGA-BRCA cohort by using the LinkedOmics database (http://www.linkedomics.org/admin.php). This database involves multiomics data within or across 32 types of cancer, including breast cancer. Results were identi ed by Pearson's correlation analysis. JASPAR (http://jaspar.genereg.net) is an online database that stores transcription factor (TF) binding pro les, was used to predict the TFs of CDCA8.
Extracted over-expressed genes of breast cancer from the TCGA database GEPIA (http://gepia.cancer-pku.cn) is a web server that supplies multiple capacities, including tumor, normal differential expressed investigation, and interactive gene analysis. GEPIA was used to compare the median expression of miR-130b in BC and normal samples and to analyze the correlation between CDCA8 and TFs. Furthermore, overexpressed genes in BC were explored by GEPIA.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis
The correlated genes with miR-130b from miRWalk http://miRwalk.umm.uni-heidelberg.de), which were input into an open-access web tool, named WebGestalt (http://www.webgestalt.org), for GO and KEGG pathway enrichment analysis. The functional enrichment network graph was drawn. GO analysis contains three terms, including biological process (BP), cellular component (CC), and molecular function (MF). We performed a KEGG pathway analysis to identify possible pathways that involved miR-130b correlated genes,. False discovery rate (FDR) ≤0. 05 were regarded as statistically signi cant.

PPI network construction and hub genes identi cation
To further explore the downstream targets of miR-130b, a PPI network was generated via the Search Tool for the Retrieval of Interacting Genes database (STRING) (https://string-db.org/). Genes in the network were reduced by searching overlapping genes between potential target genes of miR-130b and overexpressed genes in BC. Cytoscape software (www.cytoscape.org/) was selected to construct and visualize the PPI network. We used one plugin in cytoscape (CytoHubba) to compute each overlapping genes, hub genes in our study's hub genes were the top three genes.

Statistical analysis
Prism software (GraphPad, CA, USA) was used to perform statistical analysis. t-test was used to analyze the differences between and among groups. * P < 0.05; * * P < 0.01 and * * * P < 0.001 were considered statistically different

DE-miRNAs Investigation
Two microRNA expression pro les (GSE45666 and GSE26659) were selected to investigate DE-miRNAs between BC tissue and normal breast tissue. As shown in Table 1, the GSE45666 dataset included 101 BC samples and 15 normal breast samples, GSE26659 dataset was composed of 77 BC samples and 17 normal breast samples. After using the GEO2R program to analyze microRNA expression, we obtained 129 and 257 DE-miRNAs from GSE45666 and GSE26659 datasets ( Fig. 1). 10 common upregulated DE-miRNAs are shown in the Venn diagram (Fig. 2).

The prognostic value of miR-130b in BC
The prognostic value of ten common DE-miRNAs in BC was obtained from Kaplan-Meier plotter ( Fig.   3a-J). After evaluating the overall survival of the DE-miRNAs in BC, two out of ten microRNAs showed statistical signi cance, including miR-130b (p=0.012) and miR-21 (p= 0.004). Previously numerous studies are reported on the relationship of miR-21 and BC, however, these studies did not focus on the molecular mechanisms of miR-130b in BC, so we selected miR-130b as our focal point. The results demonstrated that the median survival time of highly-expressed miR-130b patients (124.53 months) was dramatically shorter than low-expressed miR-130b patients (215.20 months). miR-130b was highly expressed in BC tissue than normal breast tissue (Fig. 3k, l). miR-130b is closely related to CDCA8 Next, we explored the correlation genes with miR-130b (Fig. 4a). Top 50 positively and negatively correlated genes are visualized in Fig. 4b and Fig.4c. As shown in Fig. 4b and Fig. 5a, CDCA8 performed the highest pearson-correlation with miR-130b (Pearson-correlation = 0.6502, P < 0.01). CDCA8 expressed signi cantly higher in BC tissue compared to normal breast tissue based on GEPIA (Fig. 5b), which is in accordance with miR-130b.
THAP11 plays a key role in the upstream regulation of miR-130b The -0516 ~ -2615 bp region was de ned as a potential promoter region of the CDCA8 gene. By analyzing JASPAR data, THAP11 had the highest match score with CDCA8 among all predicted transcription factors (Fig. 6a, b). To further understand the correlation between THAP11 and CDCA8, we analyzed the data on GEPIA. The results demonstrated that CDCA8 had a positive correlation with THAP11 (P < 0.01) (Fig. 6c). Collectively, these data indicated that THAP11, as a TF, is a potential regulator that binds to and initiates the transcription process of CDCA8.
Functional and pathway analysis of miR-130b correlated genes Based on the GO term analysis, we used Go Slim to provide a general overview (Fig. 7). Gene Set Enrichment Analysis (GSEA) was conducted to further elaborate the GO terms, showing that genes correlated with miR-130b were predominantly enriched in serotonin receptor signaling pathway  Table 2). As shown in Fig. 8a, the KEGG pathway enrichment analysis results demonstrated that the correlated genes were signi cantly involved in several pathways, including the Thyroid hormone signaling pathway, ribosome, human immunode ciency virus 1 infection. TH signaling pathway is shown in Fig. 8b. Red plots indicated the miR-130b correlated genes. The signi cant KEGG pathway enrichment plots are shown in Fig. 9.
AR, KIT and ESR1 were hub genes in the PPI network We selected 323 common genes as the potential target genes of miR-130b (Fig. 11a). Based on the GEPIA database 3554 predicted over-expressed genes in BC were obtained, among which 54 genes were validated in the overlapping regions (Fig. 11b). Protein-protein interaction network is shown in Fig.11c. We identi ed hub genes with a connectivity degree in the PPI network ( Table 3). The top three genes, AR (degree=9), KIT (degree=7) and ESR1 (degree=7) were de ned as hub genes (Fig. 11d), followed by ZEB2 (degree=6), HSPA8 (degree=5), and MET (degree=4). These hub genes were all upregulated in BC.

Discussion
Dysregulation of miRNAs plays key roles in the initiation, differentiation and migration of several types of cancer, including BC (Bertoli et al. 2015). Our study aims to improve BC prognosis by validating DE-miRNAs in BC tissue compared with normal tissue and elucidate the mechanisms regulating particular target genes and pathways. In the current study, miR-130b was identi ed to be signi cantly overexpressed in human BC tissue than adjacent normal tissue. After verifying the expression, we investigated the overall survival of miR-130b in BC by using the Kaplan-Meier Plot, which indicated that miR-130b was positively correlated with a worse prognosis. These results elucidated that miR-130b might function as an oncogene when overexpressed in BC. Ana The LinkedOmics database results indicated a strong positive correlation between CDCA8 and miR-130b. We found THAP11 as a transcription factor regulating CDCA8. Accordingly, we hypothesized that THAP11 played key roles in the upstream regulatory mechanism of miR-130b, and regulated the expression level of miR-130b. THAP11, also known as thanatos associated protein 11, belongs to a unique family of TFs, which recognized particular DNA sequences through atypical zinc nger motif and regulated cell cycle and cell growth (Cukier et al. 2016 Next, we attempted to explore the downstream regulatory mechanism of miR-130b in breast cancer via predicting and analyzing genes correlated with miR-130b. A series of functional analysis was used to identify the biological functions and related signaling pathways of miR-130b in BC. Based on the data of GSEA, miR-130b was positively upregulated genes in the TH signaling pathway, such as PFKP, STAT1, SRC, and NOTCH2. Our results suggested that miR-130b was over-expressed in BC, which upregulated oncogenes, including PFKP, STAT1, SRC and NOTCH2, and played an essential role in BC progression. In our study, PFKP was analyzed as the most positively upregulated gene target in the TH signaling pathway with miR-130b (Pearson-correlation = 0.3212, p < 0.01) in BC. PFKP, also known as phosphofructokinase platelet, is the platelet isoform, which encodes a phosphofructokinase A family protein. PFKP might be a signi cant mediator for cell metabolism, participating in cancer metastasis and initiation (Lang et al. 2019). Recently, PFKP is recognized for high epidemic property in several kinds of aggressive tumors such as glioblastoma and BC (Lee et al. 2017;Moon et al. 2011). In addition, PFKP is closely related to cell survival of muscle-invasive bladder cancer (Sun et al. 2016). Wang J et al. (Wang J 2016) reported that high mRNA expression levels of PFKP were correlated with amassment of apoptotic proteins and cell proliferation in renal carcinoma, indicating the relationship between PFKP and the growth of renal carcinoma. Moreover, lots of evidence suggested that PFKP act as an oncogene in lung cancer and hepatocellular carcinoma (Park et al. 2013;Wang et al. 2015). Our study disclosed a novel mechanistic interaction between miR-130b and PFKP, supposed that the BC progression may be reduced by down-regulating miR-130b, which reduces the expression level of PFKP and constrains the cell initiation and metastasis.
Apart from PFKP, the other three oncogenes associated with miR-130b, including STAT1, SRC and NOTCH2 were worthy of conducting further investigations. Previous studies reported that most of them play a signi cant role in cell proliferation and differentiation by up-regulating in various human tumors and are closely associated with their overall survival. The expression of STAT1 is related to cell growth, regulation, immune evasion and metastasis (Ryan et  To get a deeper insight into the downstream target genes of miR-130b, we analyzed the databases from DIANA-mT, miRDB, miRWalk, Targetscan and GEPIA. 54 common genes were selected, and the PPI network suggested that AR, KIT and ESR1 were hub genes, indicating the pivotal molecular function of miR-130b. The androgen receptor (AR) played a key role in mediating androgens biological effect, that act as a prognostic indicator in tumor growth, metastasis, relapse, hormonal therapy and chemotherapy (Feng et al. 2017). Interestingly, many studies have suggested that the expression level of AR is related to higher prognostic in BC (Qu et  In summary, we revealed that the expression level of miR-130b was upregulated in BC tissue compared to normal tissue, which was signi cantly correlated with a poor prognosis of BC, signify that miR-130b functions as a signi cant oncogene in BC. Moreover, our study demonstrated that the expression level of miR-130b was regulated by THAP11 and further facilitated BC's development by regulating crucial downstream genes of the TH signaling pathway, including PFKP, STAT1, SRC, and NOTCH2. On the other hand, we predicted the potential downstream target hub genes of miR-130b, AR, KIT, and ESR1. This study for the rst time, reported the novel molecular mechanisms between miR-130b and BC to the best of our knowledge. Therefore, we identi ed that miR-130b might act as a prognostic biomarker in BC. However, further in-depth research on molecular mechanisms of miR-130b will be helpful to provide a novel therapeutic approach for the treatment of BC.

Conclusion
In conclusion, our study elucidated that miR-130b might be regulated by THAP11 and act as a potential oncogene by mediating the TH signaling pathway and the downstream target hub genes. However, further experiments are still needed to validate the molecular mechanisms of miR-130b in BC.  Figure 1 Differentially expressed miRNAs in GSE45666, and GSE26659 a-b Volcano plots show DE-miRNAs between breast cancer tissue and normal breast tissue in two microarray gene pro ling datasets. Based on |log2FC (fold change)| ≥1.5 and adjusted P-value < 0.05, the red dots imply upregulated miRNAs, and blue dots imply downregulated miRNAs, gray dots means normally expressed miRNAs.

Figure 6
Transcription factor of CDCA8 a, Prediction TF of CDCA8 from JASPAR. b, The relative score and binding site of THAP11. c, Pearson-correlation between THAP11 and CDCA8 Figure 7 Gene Ontology annotations of genes correlation with miR-130b in BC a, biological process categories. b, cellular component categories c, molecular functions categories Figure 8 Kyoto Encyclopedia of Genes and Genomes pathway a, the genes correlated with miR-130b were signi cantly involved in various pathways. b, Thyroid hormone signaling pathway diagram