The TCGA bladder cancer data (blca_tcga_pub_2017), containing genomic alternations, mRNA expression profiles, and clinicopathological features, were downloaded from the cBioPortal (http://www.cbioportal.org/) as a training group of 413 muscle-invasive BC patients. For validation, three independent datasets, including GSE13507 (165 patients: 62 MIBC and 103 NMIBC), GSE31684 (93 patients: 79 MIBC and 14 NMIBC), and GSE32548 (131 patients, the numbers of MIBC and NMIBC were undescribed), were derived from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The clinicopathological data of patients are presented in Table 1 .

Nuclear genes involved in mitochondrial disease were obtained from the MITOMAP: A Human Mitochondrial Genome Database (http://www.mitomap.org, last update: January 15^th, 2021), and all of these nuclear MTRGs (including 33 structural nuclear MTRGs and 114 non-structural nuclear MTRGs) were integrated as a nuclear mitochondria-related gene set.

The RNA-sequencing data of BC (404 samples) and normal people (28 samples) were collected from the GEPIA (http://gepia2.cancer-pku.cn/) to screen out differentially expressed MTRGs [|log2-fold change (FC)| >1, Q-value<0.01].

Nuclear mitochondria-related gene data were derived using cBioPortal to explore genomic alterations (in-frame indels, missense mutation, splice mutation, truncating mutation, amplification, deep deletion, and copy number alteration) among BC patients from the TCGA cohort.

The TCGA bladder cancer dataset was analyzed to determine whether the nuclear MTRGs alterations correlated with survival of BC patients via univariate Cox proportional hazards regression analysis, which was conducted in R studio (v. 3.4.3, https://rstudio.com/). Then, least absolute shrinkage and selection operator (LASSO) Cox regression with 10-fold cross-validation was conducted by using the “glmnet” package in R studio (21). Meanwhile, multivariate Cox regression analysis was further utilized to identify prognostic nuclear MTRGs and to construct a nuclear MTRGs signature.

The nuclear MTRGs score was calculated by the following formula: nuclear MTRGs score = gene A expression × γ_A + gene B expression × γ_B + gene C expression × γ_C +… + gene Z expression × γ_Z, where γ_Z represents the coefficient for each nuclear MTRG in the multivariate Cox regression model. The median nuclear MTRGs score served as a cutoff value to divide the patients into two groups, high and low nuclear MTRGs score groups, respectively.

The Kaplan-Meier (KM) curves were drawn by using the “survival” package in R studio. And the area under ROC curve (AUC) was calculated to evaluate the prognostic ability of the defined nuclear MTRGs signature via using the “timeROC” package in R studio. In addition, univariate and multivariate Cox regression analyses were implemented to identify the prognostic values for the signature and clinicopathological features. The nomogram and calibration plots were built by using the “rms” package in R studio.

Enrichment analysis was conducted by using the “clusterProfiler” package (22) in R studio to identify significant key genes and/or common pathways involved in tumor progression and metastasis of BC. We mainly focused on Hallmark gene set enrichment analysis (http://www.gsea-msigdb.org/gsea/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/) pathway enrichment analysis, which were visualized by using the “ggplot2” package in R studio. The threshold was defined at a p-value <0.05.

A deconvolution algorithm [TIMER, http://timer.cistrome.org/, (23)] was employed to survey the infiltrating immune cells in the tumor microenvironment (TME) between high and low nuclear MTRGs score groups, by using the transcriptomic data from the TCGA cohort. Meanwhile, the stromal score, immune score, and ESTIMATE score were calculated to further clarify the infiltrated lymphocytes (24).

A work flowchart is exhibited in Figure 1 . Besides, we performed statistical analyses of Chi-square, Fisher test, and Wilcoxon rank test in R studio. A log-rank test was used to estimate the Kaplan-Meier curves of survival analysis between the high and low nuclear MTRGs score groups. A (adjust) p-value <0.05 was considered statistically significant.

Comparison of the mRNA expression data between BC patients and normal control samples from the TCGA cohort was done, displaying that 4% (6/147) of MTRGs was differentially expressed [ Figure S1 , |log₂(FC)| >1 and Q-value<0.01]. Then we further explored the genomic alterations in these BC patients. Nearly 97% (395/408) of BC patients had altered nuclear MTRGs expression levels ( Figure S2A ). Meanwhile, the MTRGs mutation in BC patients from the TCGA cohort was analyzed, showing that 53% (218/412) of BC patients were identified with at least one MTRGs mutation ( Figure S2B ). Furthermore, copy number alterations occurred in nearly 70% (287/408) of BC patients from the TCGA cohort ( Figure S2C ).

We conducted univariate COX proportional hazards regression analysis with nuclear mitochondria-related genes, and 17 of 146 genes were found to be significantly correlated with the survival of BC patients ( Figure 2A ). The overexpression of nine nuclear MTRGs, including ATAD3, GARS, IARS2, MRPS16, NDUFS1, SUCLA2, ATPAF2, DARS2, and FRDA, significantly related to a worse prognosis of BC patients (p<0.05, Figure 2A ). On the other hand, the overexpression of eight nuclear MTRGs, including TRMU, NDUFA1, NDUFA2, COX14, COX7B, SPG7, DGUOK, and COA5, were significantly correlated with improved prognosis of BC patients (p<0.05, Figure 2A ).

The LASSO regression analysis via 10-fold cross-validation demonstrated that 16 out of 17 identified nuclear MTRGs related to patients’ survival ( Figure 2B ), and these identified genes were further analyzed with the multivariate Cox regression analysis, and their coefficients are exhibited in Table 2 . Moreover, BC patients from the TCGA cohort with a high nuclear MTRGs score (N=202) exhibited the poorer outcomes (median OS: 92.90 vs 20.20 months, p<0.0001, Figure 2C ), and the AUC for OS at 1-, 3-, 5-year was 0.71, 0.72, and 0.72, respectively ( Figure 2D ). Moreover, by comparing the AUC of our nuclear MTRGs signature with single identified MTRG for OS at 1-, 3-, 5-year, the results indicated that the integrated signature had better prediction efficacy ( Figure S3 ).

Furthermore, the signature scoring system was verified in the three independent datasets from GEO database, including GSE13507, GSE31684, and GSE32548. It was shown that BC patients with a high nuclear MTRGs score had significantly worse survival in all these independent validation cohorts (a median OS in GSE13507: 31.52 vs 98.00 months, p<0.05; GSE31684: 32.85 months vs unreached, p<0.05; GSE32548: unreached vs unreached, p<0.05; Figures 3A–C ). In addition, the AUC for OS at 1-, 3-, 5-year are as follows: GSE13507: 0.71, 0.60, 0.60 ( Figure 3D ); GSE31684: 0.67, 0.60, 0.62 ( Figure 3E ); and GSE32548: 0.71, 0.62, 0.65 ( Figure 3F ), respectively.

Meanwhile, as GSE13507 and GSE31684 datasets having comparable data of both NIMBC and MIBC, we further applied the nuclear MTRGs scoring system to figure out whether there existed some significant differences between these two subtypes of the disease. Notably, MIBC samples had a significantly higher nuclear MTRGs score than NMIBC (a median nuclear MTRGs score in GSE13507: 3.09 vs 3.06, GSE31684: 2.99 vs 2.92, p<0.05, Figure 4 ).

In the TCGA bladder cancer dataset, there were two defined histologic subtype groups (the group with papillary features or not). We observed that the BC patients without papillary histological features have significantly higher nuclear MTRGs score than those one with the papillary histological features (p<0.0001, Figure 5A ). While the BC patients without papillary-related features had the worse survival (median overall survival: 28.22 vs 44.28 months, p<0.05, Figure 5B ). In addition, the BC patients with high nuclear MTRGs score had the poorer survival (non-papillary: 19.02 vs 64.75 months, p<0.0001; papillary: 32.00 months vs unreached, p<0.001, Figures 5C, D ) no matter in the non-papillary or papillary group.

In the present study, there were 298 male BC patients, and 81/298 of these patients presented concomitant prostate cancer (PCa). By the statistical analysis, there was no significant difference of the nuclear MTRGs score between the BC patients who presented with concomitant prostate cancer and those who did not (p=0.087, Figure S4A ). Moreover, we observed no significant difference of the overall survival between the BC patients who presented with concomitant prostate cancer and those who did not (p=0.57, Figure S4B ). However, the BC patients with high nuclear MTRGs score had the poorer survival (not concomitant PCa: 17.97 months vs unreached, p<0.0001; concomitant PCa: 25.56 months vs unreached, p<0.01, Figures S4C, D ) no matter the BC patients presented with concomitant prostate cancer or not.

We then investigated the clinical features of BC patients in the high and low nuclear MTRGs score groups. Patients in the high nuclear MTRGs score group have older age at diagnosis (median age at diagnosis: 70 vs 66 years old, p<0.05, Figure 6A ). Remarkably, more BC patients with a high nuclear MTRGs score had advanced clinical stages (0 vs 1.01% in stage I, 24.26 vs 39.51% in stage II, 39.11 vs 30.50% in stage III, 36.63 vs 29.01% in stage IV, p<0.05, Figure 6B ), compared to those with low nuclear MTRGs score. According to the analysis of histological grading, more high histological grading diseases were presented in the high MTRGs score group (99.51 vs 90.45% in high histological grading, 0.49 vs 9.55% in low histological grading, p<0.05, Figure 6C ). Moreover, TNM stages were further analyzed, indicating that patients with a high MTRGs score were significantly related to aggressive tumor and metastasis ( Figures 6D–F ). We found a significantly decreased number of patients at T0–T2 but significantly more patients with disease at T3–T4 stage, which were presented in the high nuclear MTRGs score group (25.01 vs 40.44% at T0–T2, 74.99 vs 50.56% at T3–T4, p<0.05, Figure 6D ). However, there was no statistically significant difference between the two groups in lymph node metastasis (77.41 vs 77.29% at N0–N1, 22.59 vs 22.71% at N2–N3, p>0.05, Figure 6E ). By statistical analysis, only eight and three BC patients with a high or low nuclear MTRGs score were at M1 stage, respectively. A weakly significant difference was observed in the number of patients with long-distant metastasis in the high nuclear MTRGs score group (9.64 vs 2.46% at M1, p=0.05, Figure 6F ).

As is known, gender plays an important role in mitochondrial dysfunction (25, 26). However, we did not obtain evidence for significant differences in the clinical feature of gender in BC (p=0.43, Figure 7A ). Subsequently, we explored whether there existed any statistically significant differences of the identified nuclear MTRGs score between the male and female group, whereas it was found that there was no statistically significant difference in the nuclear MTRGs scores between the male and female group, either (p=0.51, Figure 7B ). Moreover, no any statistically significant differences were observed in the overall survival between the male and female group (p=0.40, Figure 7C ). In addition, this result was further validated in another three independent cohort, GSE13507, GSE31684, and GSE32548, from the GEO (p>0.05, Figures 7D–F ). But it was noteworthy that the male BC patients with high nuclear MTRGs score had the worse overall survival (median overall survival: 19.38 months vs unreached, p<0.0001, Figure S5A ). It was also verified in another three independent cohorts (GSE13507: 98.00 vs 134.97 months, p<0.05; GSE31684: 51.52 months vs unreached, p<0.05; GSE31684: 51.52 months vs unreached, p<0.01, Figures S5B–D ). For the female BC group, in the TCGA bladder cancer cohort, the female BC patients with high nuclear MTRGs score had the worse overall survival as well (median overall survival: 16.16 vs 64.75 months, p<0.0001, Figure S5E ). However, it was observed that there were no any statistically significant differences of overall survival in the cohorts of GSE13507, GSE31684, and GSE31684 (p>0.05, Figures S5F–H ).

The profiles of the top 20 most frequently mutated genes in the high and low nuclear MTRGs score group are manifested in Figure 8 . It could be seen that the prevalence of genes was distinct between the high and low nuclear MTRGs score groups. The high nuclear MTRGs group demonstrated a relatively higher prevalence of TP53 (59%, Figure 8A ), compared to that (37%) in the low nuclear MTRGs group ( Figure 8B ). However, by subsequent statistical analysis, there is no significant difference in the prevalence of TP53 as well as any other genes between the high and low nuclear MTRGs score groups (adjusted p-value>0.05, Table S1 ).

Overall, the multivariate Cox regression analysis further showed that the established nuclear MTRG-based signature was more effective to predict the prognosis of BC (p<0.001, Figure 9A ). Then a nomogram was constructed by combining the nuclear MTRGs signature with clinicopathological parameters, including age at diagnosis, clinical stage, T stage, M stage, and histological grading ( Figure 9B ). The concordance index (C-index) of the nomogram was 0.706 (95% CI=0.648–0.764). According to the nomogram, every evaluated patient would have a nomogram score that was associated with the prognosis of BC patients. Additionally, the AUC values of 1-, 3-, 5-year OS for the nomogram were 0.76, 0.75, and 0.75, respectively, showing our model had good and stable predicting ability ( Figure 9C ). And the calibration plots displayed the agreement between the predicted OS and actual OS ( Figure 9D ).

To elucidate functional differences among patients from the TCGA cohort, we further investigated the high and low nuclear MTRGs score groups. By performing the hallmark gene set enrichment analysis, it was found that high nuclear MTRGs score group were highly enriched in E2F targets, G2M checkpoint, Myc targets V1, epithelial-mesenchymal transition, mTORC1 signaling, mitotic spindle, Myc targets V2, etc (p<0.05, Figure 10A ). Besides, KEGG pathway enrichment analysis revealed that the high nuclear MTRGs score group had a significant abundance of cell cycle, DNA replication, mismatch repair, etc. (p<0.05, Figure 10B ).

The stromal score, immune score, and ESTIMATE score represented the infiltration of stromal/immune cells. The stromal score, immune score, and ESTIMATE score were all significantly higher in the high nuclear MTRGs score group (p<0.01, Figure 11A ). In addition, an estimated abundance of infiltrated immune cells via TIMER analysis revealed that CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells were significantly more enriched in the high nuclear MTRGs score group (p<0.0001, Figure 11B ), whereas a significantly higher abundance of B cells was shown in the low nuclear MTRGs score group (p<0.01, Figure 11B ).