Role of downregulated miR-133a-3p expression in bladder cancer: a bioinformatics study.

It has been discovered that miR-133a-3p acts as a tumor suppressor in bladder cancer (BC). Nevertheless, the function of miR-133a-3p in BC remains unclarified. Thus, we carried out this study to validate the expression of miR-133a-3p in BC and provide insights into the molecular mechanism underlying it. To assess the expression of miR-133a-3p in BC, we searched eligible studies from literature and Gene expression Omnibus (GEO) to perform a meta-analysis. We also plotted the summary receiver operating characteristic (SROC) curve to evaluate the diagnostic ability of miR-133a-3p in BC. Additionally, the potential target genes of miR-133a-3p were acquired from 14 online software programs and GEO database. Protein-protein interaction (PPI) network was created to identify the hub genes. Then, Gene Ontology (GO) functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were carried out to investigate the regulatory network of the target genes. From the meta-analysis, miR-133a-3p was remarkably downregulated in BC tissues compared with that in non-cancer tissues (standard mean difference =-3.84, 95% confidence interval =-6.99-0.29). Moreover, results from SROC suggested that miR-133a-3p exhibited the ability to diagnose BC (area under curve =0.8418). As for the bioinformatics study, 488 genes were chosen as the potential targets of miR-133a-3p in BC, among which 10 genes were defined as hub genes (all degrees >5). Further GO and KEGG pathway analysis indicated that the target genes of miR-133a-3p aggregated in specific biological process and pathways. In conclusion, miR-133a-3p possessed great diagnostic potential with its downregulation in BC, and miR-133a-3p might serve as a novel biomarker for BC.

Introduction

Bladder cancer (BC) is the most common malignant tumor in the urinary tract neoplasms with a high prevalence in the world, posing a threat to patients.1–10 According to the statistical data from Siegel et al, there were estimated 79,030 new cases of BC and 16,870 BC deaths in 2017.11 BC can be divided into two types: non-muscle-invasive BC and muscle-invasive BC. Treatment for non-muscle-invasive tumor consisted of transurethral resection of bladder tumor (TURBT) and post-TURBT adjuvant therapy based on risk level, while radical cystectomy and cisplatin-based chemotherapy constitute major treatments for muscle-invasive BC.12–14 However, muscle-invasive BC had high recurrence and distant metastasis in patients undergoing the aforementioned treatment.15 Therefore, there is a urgent need to seek a novel biomarker for BC.

MicroRNA (miRNA) is an endogenous non-coding small single-stranded RNA in eukaryotes, about 19–24 nucleotides in length. This small, non-coding molecule inhibited the expression of downstream mRNAs either by degrading mRNAs or suppressing the translation process by binding to the 3′untranslated site of target mRNAs.16–21 Previous researches revealed that miRNAs are involved in fundamental biological process such as proliferation, differentiation, apoptosis, metabolism and immunity, playing roles of oncogenes or tumor suppressor genes in a wide type of human cancers.22–24

MiR-133a-3p, a member of the miRNA family, plays a pivotal role in the oncogenesis of various cancers, such as ovarian cancer, colorectal cancer, breast cancer and prostate cancer.25–28 Although several studies reported that miR-133a-3p was downregulated in BC and acted as a tumor suppressor by interacting with related molecules including IGF1R, MMP-9 and EGFR,29–33 the exact mechanism of miR-133a-3p in the pathogenesis of BC was far from elucidated. Therefore, in the present study, we assessed the expression and diagnostic significance of miR-133a-3p in BC by conducting meta-analysis of literature and GSE datasets. Furthermore, bioinformatics study was performed to probe into the underlying molecular mechanism of miR-133a-3p in BC.

Materials and methods

The extraction of literature studies

To validate the expression of miR-133a-3p in BC, we conducted a systematic search in the following database: PubMed, Chinese VIP, CNKI, WanFang database, SinoMed, Embase, Web of science, Science Direct and Wiley Online Library. The search was conducted with the following keywords: (miR-133a OR miRNA-133a OR microRNA-133a OR miR-133a OR miRNA-133a OR microRNA-133a OR “miR-133a” OR “miRNA-133a” OR “microRNA-133a”) AND (malignan* OR cancer OR tumor OR tumour OR neoplas* OR carcinoma) AND (bladder OR urinary OR BC). Studies that met the following inclusion criteria were included in our studies: 1) studies provided the expression data of miR-133a-3p between BC and non-cancer tissues, 2) the subjects of the experiment were human beings and 3) the studies should be published in English or Chinese. In addition, studies were excluded according to the following exclusion criteria: 1) the expression of miR-133a-3p was detected in cell lines, 2) there were no non-cancer tissues as the control for the BC tissues and 3) the articles belonged to the types of review, meta-analysis, letters, commentaries and conference abstracts. After the selection of qualified literature studies, the basic information and data were extracted from each of the included studies: first author, ID, year of publication, country, ethnicity of the samples, sample sizes of the experimental and control group, male and female population in experimental and control group, types of specimen, types of experiment, platform and the mean (M) and standard deviation (SD) of miR-133a-3p expression in experimental and control group. The experimental design of the paper is presented in Figure 1.

Figure 1

The experimental design of the paper.

Note: The flowchart illustrates that we validated miR-133a-3p expression in BC and the underlying molecular basis via systematic review and bioinformatics study.

Abbreviations: BC, bladder cancer; GO, gene ontology; GEO, Gene expression Omnibus; PPI, protein-protein interaction; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Meta-analysis

Chip selection and data extraction

Apart from literature database, Gene expression Omnibus (GEO) database proved to be a powerful tool to mine the expression of miRNAs in human cancers. Thus, we also searched GEO database to investigate miR-133a-3p expression in BC. Eligible chip data for the meta-analysis were searched with the following strategies: (miR-133a OR miRNA-133a OR microRNA-133a OR miR-133a OR miRNA-133a OR microRNA-133a OR “miR-133a” OR “miRNA-133a” OR “microRNA-133a”) AND (malignan* OR cancer OR tumor OR tumour OR neoplas* OR carcinoma). The inclusion and exclusion criteria of chip data were the same as those of literature studies. The selected chip data were also extracted for the following basic information and data: first author, ID, year of publication, country, ethnicity of the samples, sample sizes of the experimental and control group, male and female population in experimental and control group, types of specimen, types of experiment, platform and the mean (M) and standard deviation (SD) of miR-133a-3p expression in experimental and control groups. The expression patterns of miR-133a-3p in BC and non-cancer tissues from literature studies and GEO microarray chips were visualized by GraphPad Prism 5.

Statistical analysis

Pooled standard mean difference (SMD) with 95% confidential interval (CI) of miR-133a-3p expression in BC was summarized for literature studies and GEO microarray data. Cochran’s Q test and Higgins I2 statistic were employed to check the heterogeneity of chip data. The pooled effect was calculated by random-effects model when great heterogeneity existed between chip data (P>0.05 or I2<50%); otherwise, fixed-effect model was applied (P<0.05 or I2>50%).34 To detect the source of heterogeneity, we performed subgroup analysis and sensitivity analysis to assess the influence of a single chip data on the overall SMD of the whole study cohorts. Additionally, publication bias was identified by the funnel plot from Begg’s and Egger’s test. STATA v.12.0 was used for all the statistical analysis for meta-analysis.

The diagnostic significance of miR-133a-3p

We employed MedCalc to extract the sensitivity and positivity of studies included in our meta-analysis by creating receiver operating characteristic curve for studies that provided expression value of miR-133a-3p in each specimen of BC and normal tissues. Then, we calculated true positivity (TP), false positivity (FP), false negativity (FN) and true negativity (TN) based on sensitivity and specificity from each study. Summary receiver operating characteristic (SROC) curve was generated to evaluate the diagnostic capacity of miR-133a-3p in BC by putting TP, FP, FN and TN of each study into MetaDiSc v.1.4.

Potential target genes of miR-133a-3p

To obtain the potential targets of miR-133a-3p, we used online prediction software and GEO dataset in combination. A total of 14 online prediction software, miRWalk, Microt4, miRanda, mirbridge, miRDB, miRMap, miRNAMap, Pictar2, PITA, RNA22, RNAhybrid, Targetscan, mirTarbase and PolymiRTS, were employed to collect the target genes of miR-133a-3p. Predicted genes that appear in more than five of all the online software were screened out. The differentially expressed genes (DEGs) from GSE datasets with the transfection of miR-133a-3p were another source of the targets of miR-133a-3p. Interpretation, normalization and log2 transformation were applied to process the primitive GSE data through GCBI (https://www.gcbi.com.cn). DEGs between experiment and control samples of the selected datasets were also identified via GCBI according to the threshold of fold-change >1.5 and P-value <0.05. Finally, the intersection of predicted target genes and DEGs were considered as the potential target genes of miR-133a-3p.

The construction of PPI network

The potential target genes of miR-133a-3p were put into the Search Tool for the Retrieval of Interacting Genes (STRING) to analyze the interactions between target genes. The interaction pairs with a combined score of more than 4 were illustrated by the protein-protein interaction (PPI) network composed of nodes and edges. Nodes and edges in the network represented the target genes and the interactions between them, respectively. Furthermore, hub genes were selected by calculating the degrees of the nodes.

Expression of hub genes from Gene Expression Profiling Interactive Analysis

Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) is a web-based tool that provides comprehensive expression analyses of genes on the basis of TCGA and GTEx data, which facilitated the data mining from cancer genomics.35 In this study, box plots of hub genes from PPI network in BC and normal tissues were downloaded from GEPIA for the purpose of indirectly verifying whether these genes were targeted by miR-133a-3p.

Functional and pathway enrichment analysis for target genes

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were conducted to explore the functions and metabolic pathways that the potential target genes of miR-133a-3p participated in. The functions of the target genes were aggregated according to GO terms of biological process (BP), cellular component (CC) and molecular functions (MF). Both GO and KEGG pathway analyses were performed in Database for Annotation, Visualization and Integrated Discovery. GO terms and pathways with P<0.05 were of significance. Furthermore, interactions between the function of target genes were visualized by three GO maps in Cytoscape v.5.3.0.

Consequence

The selection process for eligible studies from literature was illustrated by the flowchart in Figure 2. Finally, six studies consisting of 322 BC tissues and 95 non-cancer tissues were included to analyze the expression of miR-133a-3p in BC and non-cancer tissues.32,33,36–39 Due to the lack of sufficient data for pooling estimates, meta-analysis was not conducted for literature studies. Basic information of the selected literature studies is listed in Table 1. All the samples of the included studies were tissues from BC or non-cancer bladder. In the study by Pignot et al, samples with small and large sizes were included to evaluate the expression of miR-133a-3p through qRT-PCR.36 Yoshino et al also detected miR-133a-3p expression in samples with small and large sizes by miRNA expression signature and qRT-PCR, respectively.37 Expression data of miR-133a-3p were only available from the experiment of miRNA expression signature in Yoshino et al’s study. Similarly, Song et al assessed miR-133a-3p expression via microarray and qRT-PCR in the same set of sample.32 We only recorded expression data of miR-133a-3p from the microarray experiment considering the accessibility of the data. MiR sequencing and digital gene expression sequencing were applied to examine miR-133a-3p expression in Han et al’s study and Zhu et al’s study, respectively.38,39 As for Ichimi et al’s study, miRNA expression signatures and qRT-PCR were performed to examine miR-133a-3p expression in the same samples.33 MiR-133a-3p expression data from the group of miRNA expression signature was excluded because cell lines were included into the samples to detect miR-133a-3p expression. Data extracted from the four studies indicated that miR-133a-3p expression was obviously lower in BC tissues than in non-cancer tissues.

Figure 2

The flow chart.

Notes: The flow chart describes the process of selection of literature studies in detail. A total of six studies were included.

Table 1

The basic information and data of the included studies

First author	ID	Year	Country	Ethnicity	Sample size		Male		Specimen	Experiment type	Platform	Control M	Control SD	Exp M	Exp SD
					Control	Case	Control	Case
Pignot et al36	PMID: 23169479	2013	France	European	7/11	14/152	NA	0.83	Tissue	PCR		158.6/1.00*		0.86/0.01*
Yoshino et al37	PMID: 21304530	2011	Japan	Asian	5/10	11/23	NA	0.81	Tissue	miRNA expression signatures and PCR	TaqMan LDA Human microRNA Panel v2.0	1.17E–01		2.48E–03
Song et al32	PMID: 21133599	2010	China	Asian	25	25	0.6	0.6	Tissue	Microarray and PCR		479.04		15.42
Han et al38	PMID: 21464941	2011	China	Asian	9	9	9	9	Tissues	MiR sequencing	Illumina Cluster Station and Genome Analyze (Illumina)	1,067.726		14.5754
Zhu et al39	PMID: 22140553	2011	China	Asian	9	9	9	9	Tissues	Digital gene expression sequencing	Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA)	18.30194	10.2901	1.677944	2.326551
Ichimi et al33	PMID: 19378336	2009	Japan	Asian	31	104	0.87	0.71	Tissue	miRNA expression signatures and PCR	GeneSpring GX V7.3.1	0.96	0.196	0.048	0.016
Ratert et al70	GSE36121	2014	Germany	European	8	16	8	14	Tissue	Non-coding RNA profiling by array	GPL8227	−0.69125	0.896317	−10.2781	1.643765
Xu et al71	GSE39093	2015	China	Asian	10	10	10	10	Tissue	Non-coding RNA profiling by array	GPL8786	20.38886	4.839885	17.16362	5.542755
Meiri et al72	GSE20414	2012	Israel	Asian	2	2	NA	NA	Tissue	Non-coding RNA profiling by array	GPL10067	12.73932	0.324303	9.07	2.008183
Lu et al24	GSE2564	2012	USA	American	2	7	NA	NA	Tissue	Non-coding RNA profiling by array	GPL1986	5.865	1.223295	5.76	0.861762

Note:

Samples with larger size.

Abbreviations: Exp, experiment; miRNA, microRNA; NA, not available.

Having searched eligible chip data, a total of 4 chips, GSE36121, GSE39093, GSE20414 and GSE2564, that contained miR-133a-3p expression data between BC and non-cancer tissues were included into the meta-analysis (Table 1). As shown in Figure 3, miR-133a-3p expression was significantly lower in BC tissues than in normal tissues in the study by Zhu et al (Figure 3A) and GSE36121 (Figure 3B).39 All the expression data of miR-133a-3p from two literature studies and all the GEO microarray data yielded an overall SMD of −3.64 (−6.99–0.29) (Figure 4),38,39 which indicated that miR-133a-3p expression was downregulated in BC tissues compared with that in non-cancer tissues. Due to the significant heterogeneity between chips, random-effects model was used to pool the estimates. According to the characteristics of the chip data, subgroup analysis and sensitivity analysis were chosen for tracking the source of heterogeneity. Since only six studies were included to calculate the overall SMD, meta-regression analysis was inappropriate for identifying the source of heterogeneity.40,41 By excluding one study each time to compare pooling effects of the remaining studies with the initial one, we found that no study exerted significant influence to the whole study cohorts (Figure 5). Then, we expected to find out the source of heterogeneity by conducting subgroup analysis. We considered ethnicity and experiment type as the possible sources of heterogeneity. Unfortunately, both the subgroups of ethnicity and experiment type failed to explain heterogeneity between studies. As shown in Table 2, 95% CI of SMD for subgroups of Asian, non-Asian, other experiment and non-coding RNA profiling by array was 0, which indicated no statistical significance. With regard to the publication bias, the symmetrical funnel plot created by Begg’s test revealed that no publication bias was detected (P>0.05) (Figure 6). Now that, only one literature study and all the four GEO microarray data provided detailed data of miR-133a-3p expression in each specimen, SROC was plotted based on the five qualified studies.36 The area under curve (AUC) value (AUC =0.8418) (Figure 7) indicated that miR-133a-3p showed significant ability to differentiate BC tissues from non-cancer tissues.

Figure 3

The pattern of miR-133a-3p expression in BC tissues and normal tissues.

Note: MiR-133a-3p expression was remarkably lower in BC tissues than in normal tissues from Zhu et al’s study39 (A) and GSE36121 (B).

Abbreviation: BC, bladder cancer.

Figure 4

The forest plot of meta-analysis.

Notes: The pooled SMD of −3.64 (−6.99 to −0.29) with great heterogeneity (I2=97.0%, P<0.001) suggested that miR-133a-3p expression was significantly lower in BC tissues than in normal tissues. Weights are from random effect analysis.

Abbreviations: SMD, standard mean difference; BC, bladder cancer; 95% CI, 95% confidence interval.

Figure 5

The sensitivity analysis.

Note: The sensitivity analysis revealed that no study exerted significant influence on the overall pooling result.

Abbreviation: CI, confidence interval.

Table 2

Subgroup analysis

Subgroup	I2	SMD	95% CI
Ethnicity
Asian	97.90%	−3.79	(−8.29, 0.71)
Non-Asian	95.70%	−3.32	(−9.70, 3.05)
Experiment type
Non-coding RNA profiling by array	89.60%	−2.36	(−5.02, 0.30)
Other experiment types	98.60%	−5.96	(−13.27, 1.35)

Notes: Non-Asian subgroup contained European and American; other experiment types contained digital gene expression sequencing and miRNA expression signatures and PCR.

Abbreviations: SMD, standard mean difference; 95% CI, 95% confidential interval; miRNA, microRNA; PCR, polymerase chain reaction.

Figure 6

The publication bias.

Note: The symmetrical funnel plot indicates that no publication bias was detected.

Abbreviation: SMD, standard mean difference.

Figure 7

The SROC curve.

Note: An AUC value of 0.8418 from SROC curve revealed that miR-133a-3p was of significant diagnostic value for BC.

Abbreviations: SROC, summary receiver operating characteristic; AUC, area under curve; BC, bladder cancer.

Potential target genes of miR-133a-3p and PPI network

From the result of 14 online software, a total of 18,103 genes were predicted to be the target genes of miR-133a-3p and 2,722 genes that appeared in more than five of the online software were screened out. Two GSM microarray chips, GSM492572 and GSM492574, were included to obtain DEGs between control cell lines and cell lines transfected with miR-133a-3p. After removing the duplicate, a total of 1976 DEGs were defined as the target genes of miR-133a-3p. Eventually, a total of 488 genes overlapping in the online software and GSM microarray chips were selected as the potential targets of miR-133a-3p. Then, we generated a PPI network comprising 110 edges and 75 nodes by putting the 488 potential target genes into STRING (Figure 8). Furthermore, we identified 10 hub genes (all degrees >5), PHLPP1, ESR1, ABL1, RPS6KB1, PAK2, STK38, IGF1R, UBE2I, CDK16 and CSK, from the PPI network by calculating the degree of each node.

Figure 8

The PPI network of the target genes of miR-133a-3p.

Notes: A total of 75 nodes and 110 edges constituted the network, from which 10 hub genes were identified according to the value of degree. Nodes with different colors represent query proteins and first shell of interactors.

Abbreviation: PPI, protein-protein network.

Expression of hub genes from GEPIA

Among all the 10 hub genes from PPI network, four genes, PAK2, UBE2I, CDK16 and CSK, were observed to present a higher expression in 404 BC tissues than in 28 normal tissues (Figure 9).

Figure 9

Expression of hub genes in BC and normal tissues from GEPIA.

Notes: Expression of the hub genes were detected in 404 BC tissues (T) and 28 normal tissues (N) from GEPIA. Four of the genes, PAK2 (A), UBE2I (B), CDK16 (C) and CSK (D), were overexpressed in BC tissues than in normal tissues. Axis units are Log2(TPM+1).

Abbreviations: BC, bladder cancer; GEPIA, Gene Expression Profiling Interactive Analysis.

According to the results from GO analysis, there were 20, 16 and 5 GO terms significantly enriched by target genes in BP, CC and MF (Table 3). The 3 most significant GO terms of BP assembled by miR-133a-3p were membrane protein ectodomain proteolysis, collagen biosynthetic process and response to oxidative stress. In CC, target genes most largely gathered in the following GO terms: extracellular exosome, focal adhesion and cytosol. As for MF, ATP binding, RNA polymerase II transcription factor activity, sequence-specific DNA binding and protein serine/threonine kinase activity were the three most significant GO terms. Moreover, three GO maps in Figures S1–S3 illustrated the interactions between the functions of target genes. Each node and link in the map represented specific GO term and interactions between them. A darker color of the node reflected a greater significance of the GO term, and sizes of the nodes indicated the number of genes assembled in each GO term.

Table 3

GO enrichment analysis of the target genes of miR-133a-3p

Category	Term		Count	%	P-value
GOTERM_BP_DIRECT	GO:0006509	Membrane protein ectodomain proteolysis	6	0.007358532	9.84E–05
GOTERM_BP_DIRECT	GO:0032964	Collagen biosynthetic process	3	0.003679266	0.004966407
GOTERM_BP_DIRECT	GO:0006979	Response to oxidative stress	7	0.008584954	0.007494953
GOTERM_BP_DIRECT	GO:0031623	Receptor internalization	5	0.00613211	0.009738182
GOTERM_BP_DIRECT	GO:0007264	Small GTPase mediated signal transduction	13	0.015943486	0.010401849
GOTERM_BP_DIRECT	GO:0046777	Protein autophosphorylation	9	0.011037798	0.015688851
GOTERM_BP_DIRECT	GO:0048280	Vesicle fusion with Golgi apparatus	3	0.003679266	0.016393286
GOTERM_BP_DIRECT	GO:0045931	Positive regulation of mitotic cell cycle	4	0.004905688	0.01978581
GOTERM_BP_DIRECT	GO:0032456	Endocytic recycling	4	0.004905688	0.022591642
GOTERM_BP_DIRECT	GO:0030217	T cell differentiation	4	0.004905688	0.02560231
GOTERM_BP_DIRECT	GO:2000036	Regulation of stem cell population maintenance	3	0.003679266	0.027032891
GOTERM_BP_DIRECT	GO:0051292	Nuclear pore complex assembly	3	0.003679266	0.027032891
GOTERM_BP_DIRECT	GO:0000910	Cytokinesis	4	0.004905688	0.032235851
GOTERM_BP_DIRECT	GO:0030048	Actin filament-based movement	3	0.003679266	0.0331441
GOTERM_BP_DIRECT	GO:0072659	Protein localization to plasma membrane	5	0.00613211	0.034128257
GOTERM_BP_DIRECT	GO:0015031	Protein transport	6	0.007358532	0.037343901
GOTERM_BP_DIRECT	GO:2001241	Positive regulation of extrinsic apoptotic signaling pathway in absence of ligand	3	0.003679266	0.039735772
GOTERM_BP_DIRECT	GO:0030330	DNA damage response, signal transduction by p53 class mediator	3	0.003679266	0.046774564
GOTERM_BP_DIRECT	GO:0071560	Cellular response to transforming growth factor beta stimulus	4	0.004905688	0.047906842
GOTERM_CC_DIRECT	GO:0070062	Extracellular exosome	96	0.117736515	1.58E–06
GOTERM_CC_DIRECT	GO:0005925	Focal adhesion	23	0.028207707	3.54E–05
GOTERM_CC_DIRECT	GO:0005829	Cytosol	42	0.051509726	0.002077235
GOTERM_CC_DIRECT	GO:0043209	Myelin sheath	10	0.01226422	0.006631547
GOTERM_CC_DIRECT	GO:0016020	Membrane	37	0.045377615	0.008847182
GOTERM_CC_DIRECT	GO:0001931	Uropod	3	0.003679266	0.015496611
GOTERM_CC_DIRECT	GO:0005654	Nucleoplasm	50	0.061321102	0.016133224
GOTERM_CC_DIRECT	GO:0001726	Ruffle	6	0.007358532	0.017212833
GOTERM_CC_DIRECT	GO:0030175	Filopodium	5	0.00613211	0.020498143
GOTERM_CC_DIRECT	GO:0019898	Extrinsic component of membrane	6	0.007358532	0.0245113
GOTERM_CC_DIRECT	GO:0005739	Mitochondrion	33	0.040471927	0.025783112
GOTERM_CC_DIRECT	GO:0005902	Microvillus	5	0.00613211	0.026546689
GOTERM_CC_DIRECT	GO:0005802	Trans-Golgi network	8	0.009811376	0.02787919
GOTERM_CC_DIRECT	GO:0016324	Apical plasma membrane	10	0.01226422	0.028054824
GOTERM_CC_DIRECT	GO:0005794	Golgi apparatus	22	0.026981285	0.029322169
GOTERM_CC_DIRECT	GO:0005884	Actin filament	5	0.00613211	0.038776107
GOTERM_MF_DIRECT	GO:0005524	ATP binding	52	0.063773946	0.001227934
GOTERM_MF_DIRECT	GO:0000981	RNA polymerase II transcription factor activity, sequence-specific DNA binding	12	0.014717064	0.004972384
GOTERM_MF_DIRECT	GO:0004674	Protein serine/threonine kinase activity	14	0.017169909	0.007071712
GOTERM_MF_DIRECT	GO:004482	2poly(A) RNA binding	34	0.041698349	0.021519286
GOTERM_MF_DIRECT	GO:000154	Beta-amyloid binding	4	0.004905688	0.024757557

Abbreviation: GO, gene ontology.

We also analyzed metabolic pathways of the target genes via KEGG pathway analysis. As shown in Table 4, 13 significant pathways were recorded, among which endocytosis, AMPK signaling pathway and proteoglycans in cancer were the three most significant ones.

Table 4

KEGG pathway analysis of the target genes of miR-133a-3p

Category	ID	Term	Count	%	P-value
KEGG_PATHWAY	hsa04144	Endocytosis	19	0.023302019	4.84E–04
KEGG_PATHWAY	hsa04152	AMPK signaling pathway	12	0.014717064	8.00E–04
KEGG_PATHWAY	hsa05205	Proteoglycans in cancer	14	0.017169909	0.005393183
KEGG_PATHWAY	hsa05100	Bacterial invasion of epithelial cells	8	0.009811376	0.007318134
KEGG_PATHWAY	hsa04071	Sphingolipid signaling pathway	10	0.01226422	0.008048532
KEGG_PATHWAY	hsa04810	Regulation of actin cytoskeleton	14	0.017169909	0.008355754
KEGG_PATHWAY	hsa00510	N-glycan biosynthesis	6	0.007358532	0.013382535
KEGG_PATHWAY	hsa00230	Purine metabolism	12	0.014717064	0.013384166
KEGG_PATHWAY	hsa04520	Adherens junction	7	0.008584954	0.016781642
KEGG_PATHWAY	hsa04530	Tight junction	10	0.01226422	0.018114206
KEGG_PATHWAY	hsa04261	Adrenergic signaling in cardiomyocytes	10	0.01226422	0.026173267
KEGG_PATHWAY	hsa05131	Shigellosis	6	0.007358532	0.03783327
KEGG_PATHWAY	hsa03015	mRNA surveillance pathway	7	0.008584954	0.048738091

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Discussion

Although miR-133a-3p has been reported to inhibit the deterioration of BC and exhibit downregulated expression in BC,37,42,43 the molecular mechanism of miR-133a-3p in BC was far from elucidated. Therefore, we validated the expression and diagnostic significance of miR-133a-3p in BC and further explored the underlying molecular mechanism through bioinformatics study.

Firstly, we searched literature studies and employed high throughput database, GEO, to mine the expression of miR-133a-3p in BC and non-cancer tissues, which is one of the highlights of our study. Meta-analysis of literature studies and GEO microarray data validated that miR-133a-3p presented significantly lower expression in BC tissues than in non-cancer tissues. Nevertheless, there was distinct heterogeneity between the included studies. Even sensitivity analysis and subgroup analysis failed to detect the source of heterogeneity. We assumed that the significant heterogeneity might originate from the different numbers of samples in BC and non-cancer bladder groups between studies or the different platforms used for the evaluation of miR-133a-3p expression. The number of samples in GEO database was obviously less than that from literature studies. Particularly, only two BC tissues and two non-cancer bladder tissues were recruited in GSE20414. Additionally, the platform for examining miR-133a-3p expression differed from one another.

Though the expression of miR-133a-3p in BC has been researched by several studies, no study up to date has evaluated the diagnostic significance of miR-133a-3p in BC. Therefore, we also explored the diagnostic ability of miR-133a-3p in BC by drawing SROC curve. The AUC value of SROC implied that miR-133a-3p might serve as a future diagnostic target for BC. Nevertheless, we collected limited number of samples; future studies with larger sample size are necessary to further assess the diagnostic ability of miR-133a-3p in BC.

Now that we have found out the diagnostic significance of miR-133a-3p expression in BC, we endeavor to clarify the underlying molecular mechanism of miR-133a-3p in BC by conducting bioinformatics study. Although previous studies have elaborated on the interaction between miR-133a-3p and several downstream targets in the pathogenesis of BC,37,44,45 the miR-133a-3p-centered regulation network in BC remained elusive on account that a single miRNA could target multiple downstream genes involved in diversified biological process and signaling pathways. Thus, we firstly identified potential targets of miR-133a-3p. Since predicted genes of miR-133a-3p from online software might express in wide types of human tissues, potential targets of miR-133a-3p in BC were obtained from the intersection of DEGs in GEO microarray chips transfected with miR-133a-3p and the selected predicted genes from online software, which enhanced the reliability of our result. As a result, a total of 488 genes were identified as the target genes of miR-133a-3p, and we focused on the hub genes to facilitate our understanding of the functions of target genes in regulating the development of BC. According to the results from GEPIA, four of the 10 hub genes, PAK2, UBE2I, CDK16 and CSK, had a higher expression in BC tissues than in normal tissues. Since we have confirmed downregulation of miR-133a-3p in BC tissues, overexpression of the four hub genes in BC tissues proved indirectly that these genes might be targeted by miR-133a-3p. Among all the 10 hub genes, some of them such as PHLPP1, IGF1R and CSK played essential roles in the regulation of various cellular activities such as cell proliferation, survival, migration and apoptosis.46–48 Furthermore, IGF1R has been validated to be the target of miR-133a-3p in various cancers such as osteosarcoma, hepatocellular carcinoma and ovarian cancer.25,29,49 In BC, IGF1R was reported to be targeted by miR-145,50 which strengthened the possibility of IGF1R being targeted by miR-133a-3p. We conceived that PHLPP1, IGF1R and CSK might influence the progression of BC in several malignant aspects. Additionally, the hub genes RPS6KB1, STK38, and CDK16 engaged in specific biological processes. RPS6KB1 is a serine/threonine kinase with important roles in integrating signals related to cell growth and proliferation as well as in mediating protein synthesis.51–53 STK38 is a member of the AGC kinase family and its phosphorylation is pivotal for the regulation of apoptosis and cell cycle-related processes.54,55 CDK16 belongs to the PCTAIRE family that plays an indispensable role in cancer cell proliferation and anti-apoptosis.56–59 The growth, proliferation and apoptosis of cancer cells in BC might be affected by these genes as potential targets of miR-133a-3p. Although the genes ESR1 and ABL1 functioned mainly in breast cancer and chronic myeloid leukemia, they were also listed as the hub genes.60,61 Moreover, ESR1 expression in NMIBC was significantly different between stage and grade.62 ABL1 was considered as a candidate tumor suppressor in the study by Amira et al with its significantly reduced expression in BC tissues compared with normal bladder tissues.73 These findings suggested that ESR1 and ABL1 might also participate in the development of BC through some unknown interactions with miR-133a-3p. Further studies are needed to validate the interactions between miR-133a-3p and these hub genes via experiment.

To achieve a comprehensive understanding of the functions of the potential target genes, we further carried out GO annotation analysis and KEGG pathway analysis. It can be deduced from GO analysis that the target genes of miR-133a-3p exerted an influence on the progression of BC significantly by regulating biological processes such as membrane protein ectodomain proteolysis, collagen biosynthetic process and response to oxidative stress as well as MF such as ATP binding, RNA polymerase II transcription factor activity, sequence-specific DNA binding and protein serine/threonine kinase activity. Results from KEGG pathway analysis also indicated how target genes of miR-133a-3p functioned in BC. Significant pathways such as endocytosis, AMPK signaling pathway, proteoglycans in cancer and bacterial invasion of epithelial cells enlightened us that the target genes of miR-133a-3p might participate in these pathways to impact the initiation and development of BC. Among the 14 significant signaling pathways, several of them such as proteoglycans in cancer, AMPK signaling pathway, adherens junction and tight junction closely correlated with human cancers. AMPK signaling pathway plays crucial roles in cellular metabolism, preserving cellular energy homeostasis and cell apoptosis.63,64 There is abundant evidence supporting that AMPK is critical for the survival and growth of tumor cells.65,66 In BC, AMPK signaling pathway serves as a therapeutic target for certain drugs such as metformin and gefitinib, which reveal the importance of AMPK signaling pathway in BC. Adherens junction and tight junction are two pathways related to epithelial-mesenchymal transition (EMT), a process that confers invasive phenotype to various cancers.67 EMT is also reported to accelerate the deterioration of BC.68,69 We hypothesized that the target genes of miR-133a-3p might modulate the process of EMT to influence the progression of BC.

Although this study investigated the clinical pathological significance of miR-133a-3p and explored the potential underlying molecular mechanism, it should be noted that we only hypothesized the MF of the target genes, and greater efforts are needed to verify the relationship between miR-133a-3p and the hypothesized hub genes and signaling pathways.

In summary, the present study confirmed that miR-133-a-3p acts as a tumor suppressor in BC with its downregulated expression in BC tissues. MiR-133a-3p might perform its inhibitory function in BC by regulating network of target genes through specific signaling pathways.

The GO maps for BP of target genes of miR-133a-3p.

Notes: There are 69 nodes and 101 arrows in the map. GO terms of BP are represented by nodes, and the relationships between the terms are manifested by arrows. The depth of color indicates the significance of the corresponding term.

Abbreviations: GO, Gene Ontology; BP, biological process.

The GO maps for CC of target genes of miR-133a-3p.

Notes: There are 75 nodes and 130 arrows in the map. GO terms of CC are represented by nodes, and the relationships between the terms are manifested by arrows. The depth of color indicates the significance of the corresponding term.

Abbreviations: GO, Gene Ontology; CC, cellular component.

The GO maps for MF of target genes of miR-133a-3p.

Notes: There are 36 nodes and 41 arrows in the map. GO terms of MF are represented by nodes, and the relationships between the terms are manifested by arrows. The depth of color indicates the significance of the corresponding term.

Abbreviations: GO, Gene Ontology; MF, molecular functions.