3.1. Construction and Validation of the MRGs Prognostic Model
To comprehensively explore the role of mitophagy in BCa, we initially collected 824 mitophagy-related genes (MRGs) (with relevance scores >1) from the Genecards database (Supplementary Table 1). Subsequently, to identify MRGs related to prognosis in BCa, Kaplan-Meier survival analysis was performed in TCGA and GSE13507 cohorts, resulting in 2,165 genes associated with prognosis. We selected the intersection of 824 MRGs and 2,165 prognosis-related genes, ultimately choosing 115 MRGs (Fig. 1A). Further univariate Cox regression analysis identified 26 candidate MRGs with prognostic significance (Fig. 1B). Multiple machine learning algorithms were employed to identify model genes.
In the end, a combination of Lasso and SuperPC algorithms was selected for the final model, which exhibited the most stable performance with the smallest 5-year AUC variation across all three cohorts (Fig. 1C-E). Correlation analysis revealed the relationship between BCa prognosis and 16 model genes (Supplementary Table 2), along with their co-expression patterns (Fig. 1F). Apart from CDK5RAP3, TRIM27, and SCO2, all other MRGs were identified as risk factors for BCa prognosis (p<0.05). This analysis provided insights into the interactions among these genes.
Afterwards, heatmaps showed that low-risk group had higher expression of protective genes, while the high-risk group had higher levels of risk genes (Fig. 2A). In the TCGA cohort, the model effectively predicted patient survival, with mean AUC values for 1-year, 3-year, and 5-year survival exceeding 0.7 (0.69, 0.73, and 0.74, respectively) (Fig. 2B, C). Similarly, KM analysis and AUC supported the reliability of the model in the GSE13507 cohort, indicating its sensitivity and specificity in predicting prognosis (Fig. 2D, E). To further validate the independence of our model, univariate as well as multivariate Cox analyses were performed in the TCGA cohort (Fig. 2F, G). These two Cox regression analyses confirmed that stage and age, along with risk score were significant prognostic factors for OS in BCa patients.
Based on the prognostic model established previously, we then developed a nomogram to provide a quantitative method for predicting BCa prognosis (Fig. 2H). Calibration plots demonstrated a strong consistency between predicted and actual 1-, 3-, and 5-year OS rates (Fig. 2I). This model offers clinical utility for personalized treatment planning and decision-making.
3.2. Mutation Analysis of Prognostic Model Genes
Notably, the 16 model genes exhibited relatively frequent mutations, with a mutation rate of 16.91% in BCa patients (Fig. 3A). In terms of copy number variations (CNVs), SCD, MAP1A, CSNK2A2, ANXA2, GFPT2, and SAR1A showed more copy number gains than losses, while the remaining genes exhibited more losses than gains (Fig. 3B). The chromosomal locations of CNVs for all 16 MRGs were mapped (Fig. 3C). We then analyzed the differential expression of these 16 genes between BCa and normal tissues in the TCGA database. The results revealed significant differential expression for all genes except ANXA2, CSNK2A2, PRDX6, and MAP1A (Fig. 3D). Abnormal gene mutations, CNVs, and mRNA dysregulation suggest the disruption of mitophagy in BCa. Moreover, considering the important role of epigenetics modifications in gene expression regulation, we assessed the levels of DNA methylation in MRG promoter sequences (Fig. 3E). Compared to non-tumor tissues, the DNA methylation levels surrounding the promoters of all model genes were lower in tumor tissues, suggesting that epigenetic activation might contribute to the abnormal expression of MRGs.
Gene set enrichment analysis (GSEA) was then used to explore the potential biological functions between the groups with high or low risk. Pathways such as myogenesis, epithelial-mesenchymal transition (EMT), hypoxia, angiogenesis, and Hedgehog signaling were significantly enriched in the groups with high risk (Supplementary Fig. 1A). Conversely, pathways including interferon-α and interferon-γ responses were enriched in the groups with low risk significantly (Supplementary Fig. 1B). These findings indicate activation of stromal pathways in the groups with high risk, while the lower risk groups exhibited heightened immune activity, explaining the differences between two groups in survival value. Additionally, differential mutation analysis revealed a higher mutation rate of the key cell cycle regulator TP53 in the groups with high risk compared to the low (53% vs. 43%) (Supplementary Fig. 1C, D). Most TP53 mutations in the group with high risk were missense mutations, which impair its ability to inhibit the cell cycle and promote oncogenic functions. This likely contributes to the increased cell cycle activity observed in the high-risk group.
3.3. Correlation Analysis of Prognostic Model with Immune Landscape and Chemotherapy Response
To evaluate immune infiltration in TCGA-BLCA patients, we employed the CIBERSORT algorithm (Fig. 4A). Notably, the low-risk group showed higher infiltration of immune cells such as memory activated dendritic cells, CD8+ T cells, B cells and follicular helper T cells , while the high one exhibited increased levels of M0 and M2 macrophages (Fig. 4B). These differences highlight distinct immune microenvironments between the two groups.
New metrics like the immune phenotype score (IPS), stromal score, and ESTIMATE score were also assessed to predict patients' responsiveness to immunotherapy. We found elevated IPS scores and lower stromal and ESTIMATE scores in the group with low risk compared to the high (Fig. 4C-E), suggesting a better immunotherapy effect of the patients in the group with low risk. To further explore whether the risk model can predict BCa patient response to immunotherapy, we applied it to the IMvigor210 cohort, which received anti-PD-L1 therapy. Consistent with the TCGA and GSE13507 cohorts, the prognostic model effectively differentiated survival outcomes in these patients (Fig. 4F). Interestingly, patients with lower risk scores exhibited an "inflamed" phenotype characterized by increased immune cell infiltration and improved responses to immunotherapy. In contrast, high-risk patients predominantly showed "desert" or "excluded" phenotypes, indicative of poor immunotherapy responses (Fig. 4G). The low-risk group also had a higher objective response rate (complete or partial response) to anti-PD-L1 therapy compared to the high-risk group (Fig. 4H), underscoring the strong association between risk score and responsiveness to anti-PD-L1 therapy in BCa patients.
Utilizing data from the Genomics of Drug Sensitivity in Cancer (GDSC) database, we predicted the chemotherapy responses between the high-risk and low-risk groups in the TCGA-BLCA cohort. The high-risk group had significantly higher IC50 values for cisplatin, gemcitabine, and doxorubicin, all of which are commonly used chemotherapeutic agents in BCa (Fig. 4I-K). This suggests that patients in the low-risk group may experience better outcomes when treated with these drugs, further supporting the role of mitophagy in BCa chemoresistance.
To delve deeper into the cellular localization of the 16 model genes in BCa, single-cell analysis was conducted. Dimensionality reduction, clustering, and cell annotation revealed seven distinct cell types (Supplementary Fig. 2A, B, Supplementary Table 3), with most model genes being predominantly expressed in epithelial cells (Supplementary Fig. 2C). This suggests that aberrant mitophagy likely occurs in tumor epithelial cells. Given this, we assessed whether the expression of these genes could impact epithelial cell function. By utilizing the AUCell algorithm, epithelial cells were divided into high and low expression groups for the 16 MRGs, and differentially expressed genes (DEGs) were identified (Supplementary Table 4). GO enrichment analysis of these DEGs indicated significant alterations in the mitochondrial matrix and protein-containing mitochondrial complexes within BCa epithelial cells (Supplementary Fig. 2D). KEGG enrichment analysis indicated that the 16 model genes altered the levels of mitophagy in epithelial cells. Furthermore, significant changes were observed in drug resistance pathways, the p53 signaling pathway, and immune-related pathways between the high and low 16-MRG groups (Supplementary Fig. 2E, Supplementary Table 5). In summary, our single-cell data analysis indicated that most dysregulated MRGs were predominantly expressed in tumor epithelial cells, suggesting abnormal mitophagy in BCa tumor epithelial cells and further emphasizing the close link between BCa and mitophagy.
3.4. The Potential Role of DARS2 in Bladder Cancer Progression
To identify the most important MRG among the 16 genes for bladder cancer prognosis, we performed univariate Cox regression analysis, identifying 13 genes with HR > 1, indicating that higher expression of these genes is associated with poorer prognosis (Fig. 5A). Next, to further identify core genes closely related to bladder cancer progression, we integrated multiple single-cell datasets of BCa and extracted epithelial cells for Monocle3 pseudotime analysis. The results showed that UBC, MAP1A, DARS2, DCXR, and GFPT2 were progressively upregulated during the malignant transformation from normal epithelial cells to NMIBC-associated epithelial cells and MIBC-associated epithelial cells, suggesting that these five genes are closely related to BCa progression (Fig. 5B-D).
DepMap is a large-scale gene screening project aimed at identifying gene dependencies in cancer cells by integrating data from the Cancer Genome Atlas (TCGA) and the Cancer Cell Line Encyclopedia (CCLE)[26, 27]. In the DepMap project, CRISPR-Cas9 technology was used to knock out genes in 1,100 cancer cell lines to assess their impact on cell proliferation and survival[27]. Based on the aforementioned genes UBC, MAP1A, DARS2, DCXR, and GFPT2, we utilized DepMap to further understand their roles in various cancers, including BCa. Notably, only UBC, DARS2, and DCXR showed significant inhibition of tumor cell proliferation and survival upon knockout (Fig. 5E). Specifically, UBC and DARS2 stood out, with 113 of 333 cell lines showing dependency on UBC and 315 of 1,150 cell lines dependent on DARS2. Finally, to confirm the interaction between these core genes and mitophagy, we used Alphafold3 artificial intelligence to predict the protein-protein interaction potential between UBC, DARS2, and PINK1, a key gene in mitophagy, and selected the gene with the highest interaction potential for subsequent experimental studies (Fig. 5F, G). The results revealed that the interaction potential between DARS2 and PINK1 was significantly higher than that between UBC and PINK1. Therefore, DARS2 was chosen for further in-depth investigation.
To investigate DARS2's role in BCa further, we collected surgical pathology specimens from patients who underwent radical cystectomy at Nanfang Hospital's Department of Urology. Both DARS2 mRNA and protein expression were found to be significantly elevated in tumor tissues compared to adjacent normal tissues (Fig. 5H, I). Immunohistochemical (IHC) staining of paraffin-embedded tissue sections confirmed that DARS2 expression was markedly higher in BCa tissues than in adjacent non-cancerous tissues (Fig. 5J). These findings indicate that DARS2 is highly expressed in bladder cancer tumor tissues, likely contributing to the progression of the disease.
3.5. DARS2 Promotes G1-to-S Phase Transition by Upregulating CDK4 Expression and thereby Inhibiting Cellular Senescence
To gain further insight into the biological function of DARS2 in bladder cancer (BCa), we performed GO and KEGG enrichment analyses using the TCGA_BLCA dataset. The analyses revealed that DARS2 is predominantly enriched in pathways linked to cellular senescence, cell cycle regulation, DNA replication, and DNA damage repair (Fig. 6A). Prior studies have demonstrated a strong connection between cellular senescence, cell cycle regulation, and DNA repair processes, all of which are critical in cancer development and progression[28]. Given this, we explored whether DARS2 promotes BCa progression by modulating cellular senescence. We first assessed DARS2 expression across five BCa cell lines and utilized siRNA to knock down DARS2 expression in the T24 and SW780 cell lines, which exhibited the highest DARS2 expression levels (Fig. 6B-E). Knockdown of DARS2 significantly increased the proportion of cells positive for SA-β-gal staining, suggesting that DARS2 inhibition induces cellular senescence in BCa cells (Fig. 6F).
Since cellular senescence is characterized by irreversible cell cycle arrest and is closely associated with reduced proliferation, we further investigated the effects of DARS2 on cell proliferation and cell cycle progression[18]. Our results demonstrated that DARS2 knockdown inhibited the proliferation of BCa cells, leading to a significant arrest of cells in the G1 phase and a reduced number of cells progressing to the S and G2 phases (Fig. 6G-I). This indicates that DARS2 facilitates the transition from G1 to S and G2 phases, thereby limiting cellular senescence and promoting cell proliferation. Furthermore, we observed that DARS2 enhanced the migration and invasion capabilities of BCa cells in vitro (Supplementary Fig. 3A).
The ability of eukaryotic cells to transition from the G1 phase to the S and G2 phases is primarily controlled by the G1 phase checkpoint, in which CDK4 and CDK6 interact with cyclin D1 to play a critical role[29]. Additionally, P53, a key protein involved in the DNA damage checkpoint, induces P21 transcription, leading to G1 phase arrest[30]. Therefore, we further assessed the expression levels of CDK4, CDK6, P53, and P21. The results showed that DARS2 knockdown reduced CDK4 expression, while the expression levels of CDK6, P53, and P21 remained unchanged (Fig. 6J-L). Our results demonstrated that DARS2 promotes G1-to-S phase transition by upregulating CDK4 expression, thereby inhibiting cellular senescence in BCa cells.
3.6. DARS2 Promotes PINK1-Mediated Mitophagy in Bladder Cancer Cells
Although we initially identified DARS2 from MRGs, its role in mitophagy remains unknown. Therefore, we conducted further analyses and experimental validations. GSEA enrichment analysis of the GSE13507 dataset revealed that DARS2 is significantly enriched in mitophagy-related pathways (Supplementary Fig. 3B). Additionally, DARS2 showed significant correlations with multiple mitophagy-related genes (Supplementary Fig. 3C), with a particular focus on the interaction between DARS2 and PINK1. Our experiments further revealed that knockdown of DARS2 significantly downregulated PINK1 protein expression, although PINK1 mRNA levels remained unchanged (Supplementary Fig. 3D, E). Previous reports indicate that DARS2 encodes aspartyl-tRNA synthetase, which plays a critical role in mitochondrial amino acid synthesis[15]. Hence, abnormal expression of DARS2 may impair mitochondrial protein synthesis. Based on this biological function, we hypothesized that DARS2 might influence PINK1 mitochondrial synthesis. We established stable DARS2 knockdown BCa cells and isolated mitochondrial proteins (Supplementary Fig. 4A, B). As predicted, DARS2 knockdown significantly reduced mitochondrial PINK1 protein levels (Supplementary Fig. 4C). These results demonstrate that DARS2 promotes mitochondrial synthesis of the mitophagy-related protein PINK1.
To investigate whether changes in PINK1 expression after DARS2 regulation affect mitophagy, we first assessed the overall autophagy levels in the cells. By detecting LC3II/I, we found that autophagy levels decreased following DARS2 knockdown (Fig. 7A, B). Abnormal accumulation of ROS leads to irreversible mitochondrial damage, indirectly triggering mitophagy, an essential pathway for mitochondrial quality control[31]. Conversely, normal levels of mitophagy effectively remove damaged mitochondria and stabilize intracellular ROS levels under physiological conditions[32]. Given the close relationship between ROS and mitophagy, we first measured ROS levels before examining mitophagy. The results showed a significant increase in ROS levels in BCa cells following DARS2 knockdown (Fig. 7C). Further experiments revealed that the total number of mitochondria, as labeled by Mitotracker, significantly increased following DARS2 knockdown, while TMRE staining indicated a reduction in the proportion of mitochondria with normal membrane potential (Fig. 7D). This suggests that DARS2 knockdown may cause impaired clearance of damaged mitochondria, leading to the accumulation of depolarized mitochondria. We then used Mitotracker-green to label mitochondria and Lysotracker-red to label lysosomes. Super-resolution and confocal fluorescence microscopy showed a reduction in the colocalization of mitochondria and lysosomes following DARS2 knockdown, indicating decreased mitophagy levels (Fig. 7E, F). We further transfected BCa cells with the mt-Keima plasmid, which localizes to mitochondria, to assess the interaction between mitochondria and lysosomes. The plasmid exhibits green fluorescence in neutral pH mitochondria and red fluorescence in acidic pH lysosomes. Consistent with previous results, DARS2 knockdown reduced the colocalization of mitochondria with lysosomes, indicating decreased mitophagy (Fig. 7G). Finally, transmission electron microscopy revealed a significant decrease in mitophagy following DARS2 knockdown, with more mitochondria displaying vacuolization and abnormal swelling, likely due to impaired clearance of damaged mitochondria (Fig. 7H). Our results demonstrate that DARS2 promotes PINK1-mediated mitophagy in BCa cells.
3.7. DARS2 Inhibits Cellular Senescence and Promotes Tumor Growth by Enhancing PINK1-Mediated Mitophagy
To further elucidate the link between DARS2-regulated mitophagy and cellular senescence, we overexpressed PINK1 and added CCCP in DARS2 knockdown BCa cells. CCCP (carbonyl cyanide m-chlorophenyl hydrazone) is a protonophore that increases the permeability of the mitochondrial inner membrane to H+, causing loss of membrane potential and inducing rapid mitochondrial depolarization and mitophagy[33]. Interestingly, PINK1 overexpression and CCCP treatment not only significantly increased autophagy levels but also partially restored the decreased CDK4 expression caused by DARS2 knockdown (Fig. 8A, B).
Moreover, PINK1 overexpression and CCCP treatment significantly reduced SA-β-gal positivity and increased cell proliferation compared to the DARS2 knockdown group (Fig. 8C, D). More importantly, this treatment alleviated G1 phase arrest caused by DARS2 knockdown (Fig. 8E). These findings suggest that DARS2 promotes CDK4 expression and G1-to-S phase transition by enhancing PINK1-mediated mitophagy, thereby delaying cellular senescence and promoting BCa cell proliferation.
To assess the effects of DARS2 on tumor growth in vivo, we performed subcutaneous tumor xenograft experiments using nude mice (Fig. 9A-C). The results demonstrated that DARS2 knockdown significantly reduced both tumor size and weight, an effect that was partially reversed by PINK1 overexpression (Fig. 9D, E). Immunohistochemical analysis of the tumor tissues confirmed that PINK1 and CDK4 expression levels were markedly reduced following DARS2 knockdown, with partial restoration upon PINK1 overexpression, consistent with our earlier experimental findings (Fig. 9F). These results indicate that DARS2 knockdown exerts an anti-tumor effect by promoting PINK1-dependent mitophagy.
To further explore the clinical significance of the DARS2-PINK1-CDK4 axis, we utilized the ssGSEA algorithm to calculate the overall expression levels of the DARS2-PINK1-CDK4 axis in each sample across multiple bladder cancer cohorts, generating a Prognostic Score (Fig. 9G). Our results showed that the DARS2-PINK1-CDK4 axis is highly expressed in BCa and is significantly positively correlated with higher tumor grade, stage, and distant metastasis across multiple datasets (Fig. 9H-N). Additionally, this axis was significantly associated with poorer overall survival outcomes and suboptimal responses to various drug treatments in patients (Fig. 9O-Q, Supplementary Table 6). In summary, DARS2 promotes PINK1-mediated mitophagy, upregulates CDK4 expression, and facilitates G1-to-S phase transition, thereby delaying cellular senescence and promoting tumor growth. Targeting the DARS2-PINK1-CDK4 axis holds promise for improving the prognosis of BCa patients.