Exploring microRNA target genes and identifying hub genes in bladder cancer based on bioinformatic analysis

Background Bladder cancer (BC) is the second most frequent malignancy of the urinary system. The aim of this study was to identify key microRNAs (miRNAs) and hub genes associated with BC as well as analyse their targeted relationships. Methods According to the microRNA dataset GSE112264 and gene microarray dataset GSE52519, differentially expressed microRNAs (DEMs) and differentially expressed genes (DEGs) were obtained using the R limma software package. The FunRich software database was used to predict the miRNA-targeted genes. The overlapping common genes (OCGs) between miRNA-targeted genes and DEGs were screened to construct the PPI network. Then, gene ontology (GO) analysis was performed through the “cluster Profiler” and “org.Hs.eg.db” R packages. The differential expression analysis and hierarchical clustering of these hub genes were analysed through the GEPIA and UCSC Cancer Genomics Browser databases, respectively. KEGG pathway enrichment analyses of hub genes were performed through gene set enrichment analysis (GSEA). Results A total of 12 DEMs and 10 hub genes were identified. Differential expression analysis of the hub genes using the GEPIA database was consistent with the results for the UCSC Cancer Genomics Browser database. The results indicated that these hub genes were oncogenes, but VCL, TPM2, and TPM1 were tumour suppressor genes. The GSEA also showed that hub genes were most enriched in those pathways that were closely associated with tumour proliferation and apoptosis. Conclusions In this study, we built a miRNA-mRNA regulatory targeted network, which explores an understanding of the pathogenesis of cancer development and provides key evidence for novel targeted treatments for BC.


Background
Bladder cancer (BC) is a common urogenital malignancy and is associated with a high recurrence rate and high mortality [1]. Despite great progress in the treatment of BC, its recurrence rate is still high, and its prognosis is poor. Drug resistance and tumour metastasis often contribute to the poor prognosis of bladder cancer. Numerous studies have shown that miRNAs play an important role in these processes. It has been reported that miR-22-3p could enhance chemoresistance by targeting NET1. miRNA-124-3p suppressed cell migration and invasion by targeting ITGA3 signalling in bladder cancer [2], and the miR-626/EYA4 axis was associated with cancer proliferation and metastasis [3]. Accumulating evidence suggests that various miRNAs contribute to cancer initiation, development, and survival prognosis Open Access *Correspondence: yuhongyuan2019@163.com 1 Department of Urology, Taizhou Hospital of Zhejiang Province, Wenzhou Medical University, Taizhou 317000, Zhejiang, People's Republic of China Full list of author information is available at the end of the article and have potential applications in cancer therapy in numerous malignancies, including BC [4][5][6].
MicroRNAs (miRNAs) are a series of small endogenous single-stranded noncoding RNAs that are usually approximately 21-25 nucleotides in length. They regulate gene expression post-transcriptionally by binding to the 3'-untranslated region of target mRNAs [7]. Many studies have shown that miRNAs regulate tumour proliferation and apoptosis mainly by regulating mRNA expression [8][9][10]. miR-133 was reported to play an essential role in muscle development by regulating many genes and has been identified as an important factor in cancer development [11]. miRNA-217 could accelerate the proliferation and migration of BC by inhibiting KMT2D [12]. Thus, there is a regulatory network between miRNAs and mRNAs that can regulate the occurrence, differentiation and development of tumours.
In this study, we screened differentially expressed microRNAs (DEMs) and differentially expressed genes (DEGs) between cancer tissues and normal tissues via bioinformatic analyses. The overlapping common genes (OCGs) between DEGs and DEM-targeted genes were extracted, classified, and extensively analysed. We constructed a miRNA-targeted gene network and investigated a protein-protein interaction (PPI) network along with its hub genes and significant module. Finally, we determined which hub genes and miRNAs might promote the progression of BC, which might help researchers examine molecular mechanisms involved in disease diagnosis, pathogenesis, and prognosis, thus providing information on precise gene therapy for BC research.

Data acquisition and processing
We downloaded the gene expression profiles GSE112264 (miRNA) and GSE52519 (mRNA) from the National Center for Biotechnology Information (NCBI) GEO database (https:// www. ncbi. nlm. nih. gov/ geo). A flowchart was created to illustrate our study programme (Fig. 1). The dataset GSE112264 series contained 50 BC samples and 41 normal controls. The GSE52519 series included 9 cancer tissues and 3 normal controls.

Identification of DEMs and DEGs
The DEMs and DEGs were obtained from the GEO database using the R language 3.6.1 version limma package. |log fold change (FC)|> 1 and p-value < 0.05 were considered statistically significant for the exploration of DEGs and DEMs.

PPI network construction and analysis
We predicted the miRNA-targeted genes from the Fun-Rich software database [13,14] and identified OCGs between miRNA-targeted genes and DEGs. The PPI network was built according to these OCGs using the STRING database (http:// string-db. org). Finally, we used Cytoscape 3.7.2 [15] and its plug-in cytoHubba to visualize the PPI network.

GO enrichment analyses of the OCGs
We predicted the potential biological functions of those network genes by Gene Ontology (GO) analysis through the "cluster Profiler" and "org.Hs.eg.db" R packages. A p-value < 0.05 was considered statistically significant.

Hub gene screening and analysis
We used Cytoscape 3.7.2 and its plug-in cytoHubba to find hub genes. Hierarchical clustering of hub genes was built using the UCSC Cancer Genomics Browser (http:// genome-cancer. ucsc. edu). The differential expression analysis of hub genes was performed via the GEPIA database (http:// gepia. cancer-pku. cn/). To further investigate the potential functions of hub genes, we performed gene set enrichment analysis (GSEA) and multiple GSEA.

Statistical analysis
All statistical analyses were performed using R version 3.6.1, and a p-value < 0.05 was accepted for statistical significance.

Identification of DEMs and DEGs
Heat map analyses were performed to identify DEGs and DEMs. There were 406 differentially expressed miR-NAs, including 127 upregulated and 279 downregulated miRNAs, from the microRNA dataset GSE112264 (Fig. 2). In addition, 370 upregulated genes and 253 downregulated genes were screened from the gene microarray dataset GSE52519.

Overlapping common gene analysis between DEM-targeted genes and DEGs
We screened 5218 DEM-targeted genes according to the FunRich database. After determining the overlap between DEM-targeted genes and DEGs through a Venn diagram, we found 166 OCGs (Fig. 3a) and used the STRING database to visualize the protein relationship. Then, we built a protein-protein interaction (PPI) network (Fig. 3b) and calculated the protein node degrees through cytoHubba (Fig. 3c).

GO enrichment analyses of OCGs
We predicted the potential biological functions of these genes through GO enrichment analysis. GO analysis results were listed in Table 1. The top ten terms of GO annotation in biological process (BP), cellular component (CC) and molecular function (MF) were displayed. The BP-enriched functions were mainly associated with muscle organ development and muscle cell proliferation ( Fig. 4a, b). The CC biological functions were mainly enriched in focal adhesion and cell substrate adherent junctions (Fig. 4c, d). The MF biological functions were mainly enriched in actin binding and actin filament binding (Fig. 4e, f ).  The PPI network of OCGs was constructed using Cytoscape. Upregulated genes are marked in yellow; downregulated genes are marked in green; c The node degrees of OCGs were calculated and visualized using the cytoHubba plug-in of Cytoscape. The redder the colour, the higher the node degree

Hub gene expression and functional analysis
Hierarchical clustering analysis indicated that the hub genes were significantly differentially expressed between the normal and tumour tissues (Fig. 6). GEPIA database analysis indicated that those hub genes were significantly upregulated in the tumour tissues; however, VCL, TPM2, and TPM1 were significantly downregulated (Fig. 7). This result was consistent with hierarchical clustering analysis. Therefore, the results indicated that these hub genes were oncogenes, but VCL, TPM2, and TPM1 were tumour suppressor genes. To further investigate the potential functions of the hub genes, we performed GSEA. The high enrichment plot of those hub genes was most enriched in the "cell cycle" and "nucleotide excision repair" pathways. VEGFA was most enriched in bladder cancer, and VCL was highly related to vascular smooth muscle contraction (Fig. 8).
Multiple GSEA also showed that these hub genes with lower expression levels were most enriched in the "vascular smooth muscle contraction", "snare interactions in vesicular transport" and "abc transporters" pathways. The hub genes with higher expression levels were most enriched in the "cell cycle", "nucleotide excision repair" and "cell apoptosis" pathways ( Fig. 9). These pathways were closely associated with tumour proliferation and apoptosis.

Discussion
miRNAs regulate the progression of tumours by regulating their target genes, and some miRNAs have been identified as being involved in several types of cancer [16,17]. Therefore, it is of great significance to study relationships between miRNAs and targeted genes in BC. In our study, we screened the DEGs and DEMs through GSE52519 and GSE112264 expression profiles, respectively. There were 127 upregulated and 279 downregulated miRNAs that were differentially expressed in BC patients. We found that hsa-miR-103a-3p, hsa-miR-107, hsa-miR-130a-3p, hsa-miR-133b, hsa-miR-142-3p, hsa-miR-193b-3p, hsa-miR-22-3p, hsa-miR-29b-3p, hsa-miR-367-3p,  hsa-miR-4295, hsa-miR-4500 and hsa-miR-503-5p were potential core miRNAs that play an important role in the development and prognosis of BC. hsa-miR-29b-3p was the most significantly expressed in the tumour tissues and regulated two hub genes, VEGFA and TPM1. hsa-miR-29b-3p plays an important role in many pathophysiological processes, such as cell proliferation, apoptosis and metastasis [18][19][20], and these processes have been reported in a series of tumour types [21][22][23]. hsa-miR-29b-3p could control tumour proliferation, invasion and angiogenesis by regulating VEGFA [24]. Chou et al. [25] reported that the hsa-miR-29b-VEGFA axis was closely related to metastasis and deterioration in breast cancer. Szczyrba et al. [23] showed that induced VEGF protein expression regulated by hsa-miR-29b was involved   in cancer development. It was shown that hsa-miR-103a-3p played an important role in bladder functions and was significantly correlated with patient survival [26,27]. Moreover, some studies showed that miR-103a-3p and hsa-miR-107 were significantly upregulated in tumour tissues of BC [28], and they were positively associated with tumour stages [29].
hsa-miR-133b was also associated with the prognosis of patients [30], and hsa-miR-142-3p was upregulated in bladder urothelial carcinoma patients [31]. miR-193b-3p could downregulate the cell cycle to inhibit cell proliferation [32]. Promoting miR-4295 transcription could reduce p63α translation and enhance urothelial transformation for tumorigenesis [33]. As a result, these miRNAs have been reported to be closely correlated with the occurrence and development of BC. We selected 10 hub genes according to node degrees, such as VEGFA, VCL, TPM2, TPM1, RRM2, H2AFX, CHEK1, CEP55, CDCA8 and AURKA. In the present study, VEGFA, RRM2, H2AFX, CHEK1, CEP55, CDCA8 and AURKA were significantly upregulated; however, VCL, TPM2 and TPM1 were significantly downregulated in cancer samples of BC. Furthermore, GEPIA database analysis showed that the result was consistent with the hierarchical clustering analysis and indicated that VCL, TPM2 and TPM1 were tumour suppressor genes and the others were oncogenic genes. Our study showed that VEGFA was the core of these hub genes, which further validated its important role in tumour deterioration and development. VEGFA is a member of the VEGF growth factor family, which can promote the proliferation and migration of vascular endothelial cells. VEGFA is significantly upregulated in many cancers and is essential for tumour proliferation, invasion and metastasis [34,35]. VCL is a cytoskeletal protein associated with cell-cell and cell-matrix junctions, where it is thought to function as one of several interacting proteins involved in anchoring F-actin to the membrane. Recent research has shown that morphological and mechanical stability are directly related to the actin filament organization used by cancer cells to adapt to altered laminin-rich microenvironments [36]. VCL might affect cancer development through the focal adhesion pathway [37]. TPM2 encodes beta-tropomyosin, a member of the actin filament binding protein family that is poorly expressed in high-grade, relapsed, and metastatic prostate tumours and may be a potential prognostic biomarker in prostate cancer [38]. TPM1 is also a member of the tropomyosin family of highly conserved, widely distributed actin-binding proteins involved in the contractile system of striated and smooth muscles and the cytoskeleton of non-muscle cells. The upregulation of TPM1 inhibited cell proliferation, induced apoptosis, and inhibited cancer growth in BC cells [39]. Moreover, the migratory and invasion ability of cancer cells was enhanced by the destruction of stress fibres and by TPM1-mediated relevant adhesive structures [40,41]. Additionally, Thorsen et al. [42] also revealed that the mRNA expression level and protein expression level of TPM1 were downregulated in BC. These results supported our idea of TPM1 as a tumour suppressor. To further investigate the potential functions of these hub genes, we performed GSEA and multiple GSEA. The results showed that these hub genes were most enriched in those pathways that were closely associated with tumour proliferation and apoptosis.

Conclusions
In conclusion, we constructed a miRNA-mRNA regulatory network and showed the connections between hub genes and miRNAs in the regulatory network. This study indicated the potential mechanisms of the miRNA-mRNA regulatory network in the physiological and pathological processes of BC, and it might provide insight into therapy and improve prognosis in the future. Since the number of gene chips included in this study is not large, there may be some shortcomings and contradictions, such as data bias and differences in detection methods. Thus, we are planning to further perform a series of experiments to verify the results of this study and explore the molecular mechanism of the regulatory network.