Functional Genomics Profiling of Bladder Urothelial Carcinoma MicroRNAome as a Potential Biomarker.

Though bladder urothelial carcinoma is the most common form of bladder cancer, advances in its diagnosis and treatment have been modest in the past few decades. To evaluate miRNAs as putative disease markers for bladder urothelial carcinoma, this study develops a process to identify dysregulated miRNAs in cancer patients and potentially stratify patients based on the association of their microRNAome phenotype to genomic alterations. Using RNA sequencing data for 409 patients from the Cancer Genome Atlas, we examined miRNA differential expression between cancer and normal tissues and associated differentially expressed miRNAs with patient survival and clinical variables. We then correlated miRNA expressions with genomic alterations using the Wilcoxon test and REVEALER. We found a panel of six miRNAs dysregulated in bladder cancer and exhibited correlations to patient survival. We also performed differential expression analysis and clinical variable correlations to identify miRNAs associated with tobacco smoking, the most important risk factor for bladder cancer. Two miRNAs, miR-323a and miR-431, were differentially expressed in smoking patients compared to nonsmoking patients and were associated with primary tumor size. Functional studies of these miRNAs and the genomic features we identified for potential stratification may reveal underlying mechanisms of bladder cancer carcinogenesis and further diagnosis and treatment methods for urothelial bladder carcinoma.

Introduction

Bladder cancer is the ninth most common cancer worldwide, accounting for 430,000 new cases and 165,000 deaths per year [1]. Over 90% of bladder cancers are urothelial carcinoma, a transitional cell cancer [2]. An estimated 25% of patients are diagnosed with muscle invasive carcinoma, which has a much higher mortality rate than nonmuscular invasive bladder cancer [2]. However, the efficacy of current diagnosis and treatment methods is under question. A 2012 survey revealed that out of 4790 diagnosed patients, only one received the full scope of recommended treatments, most likely because a comprehensive treatment is overly complex and uncomfortable for patients [3]. For patients who are treated for muscle invasive cancer, poor prognosis and recurrence of tumor are common despite radical surgery or radiation treatments [2]. The discomfort involved in screenings for urothelial carcinoma, including upper urinary tract imaging and cystoscopy, serves as a barrier for early diagnosis of bladder cancer [3]. An understanding of the molecular pathogenesis of bladder urothelial carcinoma can increase the standard of care for patients by providing biomarkers that improve diagnosis efficiency, reduce treatment duration, and increase the success rate of treatments.

MicroRNAs (miRNAs) are a class of noncoding RNAs 19 to 22 nucleotides long that regulate protein-coding target genes through RNA interference and gene silencing. Because of their widespread presence and detectability, miRNAs may be promising biomarkers that can address current needs in diagnosis and treatment. miRNAs bind to complementary sequences of 3′-untranslated regions of mRNAs to block translation or to mark the mRNAs for degradation [4]. The dysregulation of miRNA and other classes of noncoding RNAs, such as long noncoding RNAs, is implicated in the carcinogenesis of a variety of cancers, including bladder urothelial carcinoma, and can act as either tumor suppressors or oncogenes [5]. Currently, patients must undergo painful cystoscopy to test for the presence of bladder carcinoma and recurrence of cancer after treatment because of the low sensitivity of less invasive urinary cytology tests [6]. However, ribonucleotides are stable and detectable in body fluids such as urine, and most circulating RNAs are reported to be miRNAs or long noncoding RNAs, as opposed to mRNAs, making miRNA potential diagnostic markers.

A number of studies have previously identified specific miRNAs upregulated or downregulated in bladder cancer [5], [7]. However, these studies are often severely limited by the number of patients involved. Most examined dysregulation of individual miRNAs in isolation from other dysregulated RNAs, with an inconsistent method of analysis employed for the individual studies [5]. To perform a more comprehensive analysis of miRNA dysregulation, we utilized RNA-sequencing data of bladder urothelial carcinoma samples from 409 patients available from The Cancer Genome Atlas (TCGA) database. A panel of six dysregulated miRNAs was identified to be differentially expressed and statistically associated with patient survival rates. These genes were further tested for association with clinical variables, mutations, and copy number variations (CNVs) to assess their clinical significance. Finally, we utilized an information theory-based bioinformatics tool to perform more rigorous correlations of miRNA expressions with genomic alterations. Besides providing potential insights into the genomic elements responsible for miRNA dysregulation, these correlations can allow for the stratification of patients into groups with similar genomic profiles through miRNA expressions to aid with diagnosis and treatment. Because tobacco smoking is a major risk factor for bladder cancer, we also examined miRNA dysregulation in smokers compared to nonsmokers so that etiology may be factored into a comprehensive miRNA panel for characterizing bladder cancer.

Material and Methods

miRNA-Sequencing Datasets and Clinical Data

Level 3–normalized miRNA expression datasets and clinical data for 409 bladder urothelial carcinoma patients were obtained on 17 April 2017 from TCGA (https://tcga-data.nci.nih.gov/tcga). Datasets from 19 solid normal tissue samples of bladder urothelial patients were also obtained. All cases were muscle invasive urothelial carcinoma. Most samples (306 patients) were obtained using transurethral resection. Forty-six samples were obtained through endoscopic biopsy. Other samples were obtained through cystectomy or cystoprostatectomy. No patients received chemotherapy or radiation therapy prior to sample collection. After sample collection, 10 patients were treated with radiation therapy, and 114 patients received chemotherapy. One hundred and nine patients died at the time of data upload. Eleven patients had distant metastasis. Patients were followed until death or to a maximum of 4684 days in this TCGA study. We separated patients into two pairs of cohorts: cancer and normal cohorts as well as smoker and nonsmoker cohorts.

miRNA Differential Expression Analyses

MiRNA read counts were extracted from the TCGA Level 3 gene expression files. The read count tables were imported into edgeR v3.5 (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html), and lowly expressed miRNAs (counts per million <1 in more than the number of samples in the smaller cohort) were filtered from the analysis. Following TMM (trimmed mean of M-values) normalization, pairwise designs were applied to identify significantly differentially expressed miRNAs in 1) bladder urothelial carcinoma tumor tissues versus normal tissues and 2) tumor tissues from lifelong nonsmoking patients versus tumor tissues from patients with history of smoking. All miRNAs identified as differentially expressed in each edgeR comparison were retained as candidates.

Association of Candidate miRNAs with Patient Survival and Clinical Variables

Survival analyses were performed using the Kaplan-Meier model, with miRNA expression in bladder urothelial carcinoma tumors designated as a binary variable based on expression above or below the median. The time to censoring or event was designated as days to last follow-up or days to death, respectively. Univariate Cox regression analysis was used to identify candidates significantly associated with patient survival (P < .05). Survival-correlated miRNAs between cancer and normal cohorts were evaluated for clinical significance. Employing the Kruskal-Wallis test, we investigated miRNA association with neoplasm grade, pathologic stages, and nodal extracapsular spread using clinical data and miRNA expression values (counts per million) from bladder urothelial carcinoma patients. The UICC/AJCC TNM staging system was used for pathological stages. In pathological T stage analysis, patients with stages T1a and T1b were grouped into stage T1, and likewise for stages T2, T3, and T4. Patients with no available information for a given characteristic were filtered from analyses involving that variable. Multivariate Cox analysis was then performed to evaluate whether correlations were independent of clinical variables such as smoking history, gender, pathological stages, and presence of lymphovascular invasion.

Association of Candidate Survival-Related miRNAs with Common Genomic Alterations

Mutation and CNV data for the bladder urothelial carcinoma tumors were obtained from mutation and CNV annotation files generated by the Broad Institute GDAC Firehose on 16 May 2017. We focused our analysis on the 26 most frequently mutated genes in urothelial bladder carcinoma and 39 loci of frequent CNV, as determined by Weistein et al. [7]. Wilcoxon rank sum tests were employed to test for significant associations between miRNA expression level (counts per million) and mutational and CNV status.

Information-Coefficient–Based Correlation of miRNA Expression with Genomic Alterations

The annotation files for mutation and CNVs were compiled into a binary input file for the program REVEALER (repeated evaluation of variables conditional entropy and redundancy), designed to computational identify a set of specific CNVs and mutations most likely responsible for the change in activity of a target profile [8]. The target profile was defined in our study to be miRNA expression. In order to identify a set of most relevant genomic alterations, REVEALER runs multiple iterations of the correlation algorithm, with the genomic feature exhibiting the strongest correlation in each run serving as a seed for the successive run. We set the maximum number of iterations to three. A seed is a particular mutation or copy number gain or loss event that most likely accounts for the target activity. When given a seed, REVEALER will focus correlation on only patients with altered target activity not accounted for by the seed. Since we do not know which genomic alteration is responsible for the dysregulation of each miRNA, we set the seed for the first iteration to null. We set the threshold of genomic features to input to features present in less than 75% of all samples.

Results

Identifying Dysregulated miRNA in Bladder Urothelial Carcinoma

To determine the extent of miRNA dysregulation in bladder urothelial cancer patients, miRNA read-counts data of cancer genomes from 409 patients were downloaded from TCGA. These miRNA expression data were compared against the miRNA read counts of 19 normal tissue genomes, with each normal sample paired with a cancer sample from the above 409 patients. After performing differential expression analysis using the negative-binomial model, 236 miRNAs were identified to be differentially expressed in cancer tissue compared to normal tissue ( < .05; Figure 1). The P value was corrected for multiple hypothesis testing using the Bonferroni method. We focused our analysis on 218 miRNAs that had a fold change greater than 2 and FDR < 0.0005.

Figure 1

Summary of differential expression results. (A) Heatmap of significantly dysregulated miRNAs (FDR< 0.05) in normal tissues (n=19) versus in bladder urothelial carcinoma tissues (n=409). Low expression is represented by solid green; dark green to black represents medium expression; dark red to solid red represents higher expression. (B) Circular representation (Circos Plot) of miRNA differential expression between normal tissue cohort (inner ring) and cancer tissue cohort (outer ring). Each red bar represents the averaged expression within a cohort of a particular miRNA, arranged according to its position within the genome.

A large number of previously identified dysregulated miRNAs in bladder carcinoma were confirmed to be differentially expressed in our analysis, including miR-1 [9], miR-96 [10], mir-143 [11], miR-133a [9], mir-490 [12], mir-195 [13], and miR-99a [14], in order of decreasing fold change. All of the above miRNAs showed greater than four-fold difference in our analysis, and their dysregulation (upregulation or downregulation) corresponded with previously published results. However, most dysregulated miRNAs we identified have not been shown to be associated with bladder cancer in published literature.

Assessing the Correlation of Dysregulated miRNA with Patient Survival

The expression levels of the 218 candidate miRNAs were analyzed for correlation with patient survival using the Kaplan-Meier estimate model. The Kaplan-Meier method uses right-censored observations to account for patients that have yet to reach the outcome of interest, which is death in survival analysis, at the time of data collection. Expression levels of the candidate miRNAs were divided into high and low expressions at the median expression level, and a Kaplan-Meier regression curve was plotted for each level (Figure 2). The Cox proportional-hazards regression model was applied to determine whether the differential survival rate for the two curves was statistically significant. The expression levels of 44 miRNAs were found to have significant correlation with patient survival rates (P<.05). Six miRNAs exhibited correlation with patient survival consistent with the direction of dysregulation (i.e., high expression of upregulated miRNAs in bladder cancer samples corresponded with poorer survival prognosis, while low expression of downregulated miRNAs corresponded with poorer prognosis) (Table 1). The upregulation of four miRNAs, miR-483, miR-33b, miR-33a, and miR-3065, correlated with lower patient survival rates, suggesting the potential oncogenic roles of these miRNAs. The downregulation of two miRNAs, miR-378c and miR-30e, correlated with lower patient survival rates, suggesting tumor-suppressing roles. All six miRNAs have not been implicated in bladder urothelial carcinoma in previous studies. However, with the exception of miR-483, these miRNAs have well-documented correlations with the development of other cancers.

Figure 2

Kaplan-Meier survival plots depicting survival outcomes of bladder urothelial carcinoma patients associated with six differentially expressed miRNAs. The two curves in each plot denote high and low miRNA expression, respectively.

Table 1

List of miRNAs with Significant Differential Expression between Cancer and Normal Samples and Significant Correlation with Patient Survival

	Fold Change	P Value	FDR
miR-483	22.26	2.8E-03	1.8E-05
miR-33b	6.76	5.2E-06	6.4E-08
miR-33a	5.51	1.5E-09	3.0E-11
miR-3065	3.58	9.0E-03	4.8E-05
miR-378c	0.39	4.5E-07	6.9E-09
miR-30e	0.49	5.7E-10	1.3E-11

Correlating miRNA Expression Levels with Clinical Variables

The six dysregulated miRNAs identified in patient survival analysis were analyzed for any relationship between their expression and various clinical variables, including pathological stage, neoplasm grade, and primary tumor size (pathological T stage). The Kruskal-Wallis test was used to determine the significance of correlation (P<.05). Higher miR-33a and miR-33b expression levels were correlated with a greater extent of metastasis, which further suggested their oncogenic role (Figure 3, A and B). Lower miR-30e and miR-378c expression levels were correlated with higher tumor grade, higher pathological stage, and greater primary tumor size, suggesting a tumor-suppressing role for these miRNAs (Figure 3, C and D).

Figure 3

Correlation of candidate miRNAs with clinical varibales. (A-B) Boxplots depicting significant correlation between expression levels of survival-associated miRNAs (A) miR-33a and (B) miR-33b and distant metastasis status. (C-D) Boxplots depicting significant correlation of (C) miR-30e and (D) miR-378c expression levels with neoplasm grade, pathological T stage, and pathological stage.

Multivariate Cox regression analysis revealed that the association of miR-30e downregulation and miR-33a, miR-33b, and miR-483 upregulation with poor patient outcome was largely independent of clinical prognostic factors such as smoking history, gender, clinical stages, and lymphovascular invasion (P<.05; Supp. Table 1).

Associating miRNA Expression Levels with Common Genomic Alterations in Bladder Urothelial Carcinoma

The potential function of the survival-correlated miRNAs could be elucidated by associating their expression levels to changes in the genome that can drive the pathogenesis and progression of bladder urothelial carcinoma, such as CNV or mutational changes.

Somatic CNVs can be associated with the activation or inactivation of miRNAs in cancers, with a large fraction of miRNAs located within genomic sites that undergo frequent copy number gains or losses [15]. We focused on analyzing correlation of survival-related miRNA expressions with CNVs of 39 loci of the genome where high CNV was previously found in bladder cancer [7]. The Wilcoxon rank sum test was applied to determine the significance of correlation (P<.05). CNVs in 32 unique loci correlated with miRNA expression levels, with the most significant correlations caused by the deletion of 2q22.1 and 3p14.2, which correlated with dysregulation of miR-33a, miR-33b, and miR-30e (Figure 4, A-C). The amplification of 3q26.2 correlated with the increased expression of miR-33a, mR-33b, and miR-3065 (Figure 4, B-D). Other significant correlations were shown in Figure 4.

Figure 4

Boxplots depicting significant correlations of survival-associated miRNAs’ expressions with CNV. CNVs of loci 3p14.2, 2q34, and 2q22.1 were significantly associated with the expression of (A) miR-30e, (B) miR-33a, and (C) miR-33b. CNV of locus 3q26.2 was associated with the expression of (B) miR-33a, (C) miR-33b, and (D) miR-3065. Significant associations with CNVs of other loci were seen for (E) miR-378c expression. Samples with expression values greater than three times the interquartile range above the median were considered outliers and not displayed.

A number of mutations have shown to be actively involved in the pathogenesis of bladder cancer, although the underlying molecular mechanism of cancer-related pathways is poorly understood and thought to be more complex than previously conceived. The most common mutations associated with bladder urothelial carcinoma, including mutations of TP53, MLL2, FGFR3, and 23 other genes, were the focus of our analysis [8]. Using the Wilcoxon rank sum test (P<.05), we found that the downregulation of miR-30e and miR-378c, the two miRNAs found to be downregulated in bladder cancer tissues in our analysis, correlated with the presence of TP53 mutation (Figure 5, A and C). TP53 has been identified as a target of the miR-30 family as well as associated to miR-30e in other cancers, but the mechanism of miR-30e interaction with TP53 has yet to be explored [16]. Additionally, the presence of mutations in RHOA was associated with elevated expression of miR-30e and lowered expression of miR-33b, while presence of mutations in FGFR3 was associated with elevated expressions of miR-30e and miR-378c (Figure 5).

Figure 5

Mutational correlations using Wilcoxon rank sum. Boxplots depicting significant correlations of the expressions of survival-correlated miRNAs (A) miR-30e, (B) miR-378c, and (C) miR-33b with mutation status of commonly mutated genes in bladder urothelial carcinoma, including FGFR3, TP53, and RHOA.

Identification of Dysregulated miRNAs Implicated in Smoking and Correlation with Clinical Variables

Since tobacco smoking is the most widely established risk factor in bladder cancer, we attempted to correlate miRNA dysregulation with patient smoking status to explore the potential mechanisms of bladder cancer pathogenesis due to smoking. Differential expression analysis was performed on smokers’ cancer tissues versus lifelong nonsmokers’ cancer tissues (Figure 6A). The six survival-correlated miRNAs identified above were not implicated with smoking in the differential expression analysis results. Twenty-nine miRNAs were identified to be differentially expressed between smokers and nonsmokers (P<.05 and FDR<0.05). Twenty-one miRNAs with fold change greater than two were selected for correlations with clinical variables (Table 2). miR-431 and miR-323a appeared to be upregulated in cancer tissues of smokers versus that of nonsmokers, and their increased expressions correlated with increased primary tumor size, suggesting an oncogenic role for these miRNAs (Figure 6B). Increased miR-431 expression also correlated with increased tumor invasion of surrounding lymph and vascular tissues (Figure 6B).

Figure 6

Candidate miRNAs associated with tobacco smoking. (A) Heatmap of miRNA expression (in counts per million) of cancer samples from nonsmoking patients compared to that of cancer samples from smoking patients. In order from low expression to high expression, the heatmap colors are green, dark green, black, dark red, and red. (B) Boxplots depicting correlations between the expressions of smoking-associated miRNAs miR-323a and miR-431 and pathological T stage and lymphovascular invasion.

Table 2

List of miRNAs with Significant Differential Expression between Cancer Samples from Smokers and Nonsmokers

	Fold Change	P Value	FDR
miR-527	6.50	3.4E-05	4.2E-06
miR-518a-2	5.66	1.9E-04	1.9E-05
miR-518c	5.60	1.7E-04	1.8E-05
miR-518e	4.79	1.4E-03	1.1E-04
miR-525	4.45	3.9E-03	2.6E-04
miR-517a	4.10	8.9E-03	4.8E-04
miR-517b	4.08	9.5E-03	4.8E-04
miR-516a-1	3.99	4.2E-04	3.8E-05
miR-520a	3.97	9.5E-03	4.8E-04
miR-194-2	0.27	3.7E-29	3.7E-29
miR-526b	3.67	6.8E-03	4.0E-04
miR-194-1	0.28	1.7E-27	8.4E-28
miR-516a-2	3.60	1.8E-03	1.3E-04
miR-518b	3.55	2.0E-02	8.5E-04
miR-522	3.34	2.9E-02	1.1E-03
miR-519a-2	3.33	1.7E-02	7.7E-04
miR-217	3.28	4.2E-06	7.0E-07
miR-215	0.32	3.7E-06	7.0E-07
miR-431	2.72	2.6E-07	6.5E-08
miR-192	0.39	3.2E-13	1.1E-13
miR-323a	2.51	1.2E-05	1.8E-06

Computational Correlations for Comprehensive Analysis of All Known Bladder Urothelial Carcinoma Genomic Alterations’ Association with miRNA Expression

The REVEALER algorithm, developed by Kim et al., was used to computational correlate all genomic alterations present in bladder cancer samples to miRNA expression [8]. This exhaustive process aims to identify genomic alternations that are functionally responsible for the range of miRNA expression values in the patient samples, as opposed to associations with the general potential driving processes of bladder cancer that we performed in a previous analysis. Utilizing the concepts of mutual information and information coefficient from the information theory, REVEALER calculates a conditional information coefficient (CIC) for each genomic alteration event. When the absolute value of the CIC is around 0.30 or above, the correlation between the event and miRNA expression is reasonably established, based on simulated test data sets [8]. The most significant correlations are displayed in Figure 7.

Figure 7

Association of genomic alterations with miRNA expression using REVEALER. The gradient bar visualizes the expression distribution of miRNAs. The dark red extreme represents the highest expression, while the dark blue extreme represents the lowest expression. Genomic alterations are analyzed for correlation with low miRNA expression in (A) downregulated miRNAs, while they are analyzed for correlation with high miRNA expression in (B) upregulated miRNAs. Thus, the CIC is negative for correlations with miR-30e and miR-378c expressions. Significant correlation can be visually examined by the increased density of dark blue vertical bars to the extreme left of each row compared to the extreme right. Each vertical bar represents the presence of the stated genomic feature in a specific patient. However, care must be taken when asserting correlations based on visual inspection because of the nonlinear distribution of expression values in each plot. Significant correlation was determined with CIC of around 0.30 or higher. A red dividing bar in between genomic features signifies a change in iteration. Results from iterations yielding no significant correlations are not displayed.

Discussion

Noncoding RNAs, including miRNAs, are known to be involved in cancer pathogenesis pathways. The mechanisms of action of miRNA are relatively well studied and characterized relative to other noncoding RNAs such as long noncoding RNAs, making miRNAs promising candidates to serve as potential disease markers [5]. New biomarkers for diagnostic purposes and for prognostic prediction are urgently needed for bladder cancer since, currently, there are no biomarkers that have adequate sensitivity or specificity, and no new treatment method has been approved for more than 30 years. However, since the pathogenesis and molecular mechanisms of bladder cancer have been extensively documented, investigating regulatory molecules such as noncoding RNAs and their association with known mechanisms in bladder cancer initiation or development may yield insights into their clinical significance and function.

We identified 236 miRNAs to be differentially expressed in a patient cohort of 409 bladder urothelial carcinoma patients (FDR<0.05). Out of these candidates, six miRNAs, miR-483, miR-33b, miR-33a, miR-3065, miR-378c and miR-30e, correlated with patient survival through univariate Cox regression analysis (P<.05). miR-378c and miR-30e were downregulated in cancer patients, while the other four miRNAs were upregulated. To the best of our knowledge, none of these miRNAs were previously found to be implicated in bladder cancer. Multivariate Cox regression analysis was performed on pathological staging, smoking history, lymphovascular invasion, and gender for these six miRNAs. The correlation between expression and poor prognosis was independent of the above variables for miR-30e, miR-33a, miR-33b, and miR-483.

The survival-correlated miRNAs were associated with somatic CNVs in 32 loci. CNV of locus 2q22.1, the site of LRP1B and frequently deleted in bladder cancer [7], was strongly implicated with the upregulation of miR-33a, miR-33b, and miR-30e (P<.003, Wilcoxon rank sum). LRP1B, an LDL receptor-related protein, is well known as a tumor suppressor across multiple cancers. LRP1B was found to be a target of miR-30e in broiler chickens, but no association between LRP1B and miR-33a, miR-33b, or miR-30e has been find in humans [17]. CNV of locus 3p14.2, containing the gene FHIT and also frequently deleted in bladder cancer [7], was associated with the downregulation of miR-30e and upregulation of miR-33a and miR-33b in our analysis. FHIT is a known tumor suppressor that was found to mediate the increased expression of miR-30c, from the same family as miR-30e, in lung cancer to suppress genes that promote metastasis [18].

The expression levels of the six survival-related miRNAs were associated with presence of mutations in the genes TP53, RHOA, and FGFR3. Forty-nine percent of bladder urothelial carcinoma patients have a mutation in TP53 [7], a well-studied tumor suppressor gene that is known to have a complex interaction with miRNAs and is both a regulator of miRNA expression and target of miRNAs. The presence of mutation in TP53 correlated with the downregulation of miR-30e and miR-378c, both of which are expected to have tumor suppressive roles in bladder cancer from our differential expression analysis data. Studies in colon carcinoma have suggested that miR-30e increases the expression of the tumor suppressor protein p21 that is downstream of TP53 [19]. The TP53-p21 pathway was also found to be implicated with miR-378c [20]. The presence of mutation in RHOA correlated with the upregulation of miR-30e and the downregulation of miR-33b. RHOA appears to have an oncogenic role in bladder urothelial carcinoma, and the expression of miR-33b was reported to activate RHOA in melanoma cells [21], [22]. FGFR3 is also frequently reported to be mutated in bladder cancer, with 12% of a sample of over 400 patients possessing a mutation of this gene [7]. Mutation in FGFR3 was correlated with the upregulation of miR-30e, miR-378c, and miR-483 in our analysis, but no studies have linked these miRNAs to FGFR3. Collectively, these findings suggest the possibility of using single or a combination of miRNAs to detect genomic alterations in bladder cancer, underscoring the importance of further investigating the underlying mechanisms of genes within these genomic regions in regulating miRNAs to better understand their roles in bladder urothelial carcinoma carcinogenesis.

Because smoking is a major risk factor of bladder cancer, we correlated miRNA expression to smoking history to explore the effects of smoking on miRNA dysregulation and how this relationship can potentially contribute to oncogenesis. miR-431 and miR-323a were significantly upregulated in smoking patients’ tumor tissues compared to nonsmoking patients’. Their upregulation also correlated with increased primary tumor size (pathological T stage). The implication of smoking in the dysregulation of miRNAs has been studied in lung cancer but not bladder cancer, to our best knowledge [23]. miR-323a was reported to be upregulated in smoker versus nonsmoker lung cancer patient cohorts, consistent with our results in bladder cancer [24]. miR-431 was reported to be dysregulated in the lung cells of mice exposed to smoke, although it was identified as having a tumor suppressing role [25]. As in lung cancer, miRNAs may play a role in the carcinogenesis of bladder cancer as well through their interaction with carcinogens in cigarette smoke.

We next used the REVEALER program to perform a more comprehensive correlation of miRNA expression with genomic alterations beyond those that are most common and previously determined to be associated with bladder cancer. The landscape of CNVs and mutations revealed to have significant association with miRNA expression by REVEALER is very different from the associations observed from correlations using the Wilcoxon rank sum test. With the exception of FGFR3 mutation, no mutations neared the threshold of 0.30 for significant correlation, and many fewer CNVs exhibiting significant correlation with each miRNA expression were observed with REVEALER compared to that observed with the Wilcoxon test. These results suggested that REVEALER may be more selective and statistically robust than the Wilcoxon test. However, differences in the statistical approach discourage the assumption that REVEALER is always more conservative in its correlations. For example, FGFR3 mutation is less strongly associated with the dysregulation of miR-483 than the dysregulation of miR-30e and miR-378c according to the Wilcoxon P value, but FGFR3 mutation is more significantly associated with miR-483 dysregulation according to the REVEALER conditional information coefficient, while association of FGFR3 mutation with miR-30e and miR-378c dysregulation was not near the significant threshold of 0.30. This outcome suggested that both methods may be of value when determining associations of miRNA expression with genomic features.

A unique feature of REVEALER is its selection of the genomic alteration with the strongest correlations as a seed for subsequent iterations and then penalizing the CIC for genomic alterations present in patients who also possess the seed alteration, with the reasoning that the dysregulation of miRNA in those patients would be better accounted for by the seed [8]. This approach can lead to the discovery of genomic alterations that may account for a portion of dysregulations but were not recognized as significant by traditional approaches such as the Wilcoxon test because of the researcher’s interest in stronger correlations exhibited by other genomic alterations. More importantly, this innovative feature of REVEALER allows for the stratification of patients based on their genomic abnormalities using miRNAs as a biomarker.

The heterogenous nature of tumors across different patients necessitates the need for identification of specific patient cohorts that respond differently to treatment but prohibits the easy detection of such clinically relevant tumor subtypes. Mutations and somatic CNVs are genomic alterations critical to the mechanism of cancer development [26], [27], and we hope the dissection of their heterogeneity within cancers may lead to informative stratifications of patients. Great advancements have been made in developing methods to determine which CNVs and mutations are functionally relevant in the initiation and progression of cancers [28], [29], but the use of biomarkers to stratify patients based on these critical CNVs and mutations was less well studied. Most studies on bladder cancer biomarkers focused on identification of the tumor itself instead of tumor stratification. To illustrate, by testing for a panel of three miRNAs (miR-15b, miR-135b, and miR-1224-3p) isolated from urine samples, the sensitivity for diagnosing bladder cancer was found to be 94%, but the specificity was only 51% [30]. Because of the high heterogeneity in clinical characteristics of individual diagnoses, adequate specificity in diagnosing bladder cancer using miRNAs as biomarkers remains a challenge. Stratification of tumors during diagnosis, using personalized diagnosis methods that link the presence of specific miRNAs to specific mechanisms of bladder cancer pathogenesis rather than using miRNA presence as a general marker of bladder cancer presence, may improve diagnostic accuracy. Our study demonstrates a workflow that allows for the specific identification of potentially valuable miRNAs for use as biomarkers and the specific CNVs and mutations they may be able to detect.

The following are the supplementary data related to this article.

Supplementary Table 1

Multivariate Cox Regression Analysis Revealed that miRNA Expression Correlation with Survival Was Independent of Clinical Factors including Gender, Smoking History, Presence of Distant Metastasis (Pathological M Stage), Presence of Regional Lymph Node Invasion (Pathological N Stage), Presence of Primary Tumor (Pathological T Stage), and Lymphovascular Invasion for miR-30e, miR-33a, miR-33b, and miR-483.

Author Contributions

W. M. O. conceived and designed the experiments; W. T. L. and V. N. performed the experiments and analyzed the data; W. T. L. wrote the paper; H. Z. and J. W. R. reviewed the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.