Integrating machine learning models with multi-omics analysis to decipher the prognostic significance of mitotic catastrophe heterogeneity in bladder cancer

Abstract

Background

Mitotic catastrophe is well-known as a major pathway of endogenous tumor death, but the prognostic significance of its heterogeneity regarding bladder cancer (BLCA) remains unclear.

Methods

Our study focused on digging deeper into the TCGA and GEO databases. Through differential expression analysis as well as Weighted Gene Co-expression Network Analysis (WGCNA), we identified dysregulated mitotic catastrophe-associated genes, followed by univariate cox regression as well as ten machine learning algorithms to construct robust prognostic models. Based on prognostic stratification, we revealed intergroup differences by enrichment analysis, immune infiltration assessment, and genomic variant analysis. Subsequently by multivariate cox regression as well as survshap(t) model we screened core prognostic gene and identified it by Mendelian randomization. Integration of qRT-PCR, immunohistochemistry, and single-cell analysis explored the core gene expression landscape. In addition, we explored the ceRNA axis containing upstream non-coding RNAs after detailed analysis of pathway activation, immunoregulation, and methylation functions of the core genes. Finally, we performed drug screening and molecular docking experiments based on the core gene in the DSigDB database.

Results

Our efforts culminated in the establishment of an accurate prognostic model containing 16 genes based on Coxboost as well as the Random Survival Forest (RSF) algorithm. Detailed analysis from multiple perspectives revealed a strong link between model scores and many key indicators: pathway activation, immune infiltration landscape, genomic variant landscape, and personalized treatment. Subsequently ANLN was identified as the core of the model, and prognostic analysis revealed that it portends a poor prognosis, further corroborated by Mendelian randomization analysis. Interestingly, ANLN expression was significantly upregulated in cancer cells and specifically clustered in epithelial cells and provided multiple pathways to mediate cell division. In addition, ANLN regulated immune infiltration patterns and was also inseparable from overall methylation levels. Further analysis revealed potential regulation of the MIR4435-2HG, hsa-miR-15a-5p, ANLN axis and highlighted a range of potential therapeutic agents including Phytoestrogens.

Conclusion

The model we developed was a powerful predictive tool for BLCA prognosis and revealed the impact of mitotic catastrophe heterogeneity on BLCA in multiple dimensions, which then guided clinical decision-making. Furthermore, we highlighted the potential of ANLN as a BLCA target.

Introduction

Against the backdrop of a global cancer rate of 1 in 5 people under the age of 75 [1], bladder cancer contributes about 380,000 newly diagnosed cases and 150,000 deaths annually [2]. Bladder cancer is primarily mediated by abnormal proliferation of the uroepithelium and may invade the submucosal muscular layer [3]. Although nearly 80% of patients are initially diagnosed as non-muscle invasive and show a five-year survival rate of up to 85% [4], metastatic spread mediated by its extreme aggressiveness and muscle invasion contribute to a steep decline in survival rates to as low as 6% [5]. Previous studies have explored numerous BLCA-related drivers [6], but the development of high-resolution prognostic stratification models as well as heterogeneity studies remain necessary to advance the management of BLCA and require urgent attention to address.

Mitosis is well known to the public as the main process of eukaryotic cell proliferation in organisms, with cycles that include DNA replication, nuclear division, and more [7]. Malignant tumor cells, on the other hand, exhibit marked cell cycle dysregulation, presented as an over-activated mitotic cycle [8]. The mitotic catastrophe (MC) then acts as a fundamentally endogenous tumor suppressor pathway, propelling the cell to the grave of death when failure of mitosis is perceived [9]. In addition, MC has been reported to be mediated by complex cascade pathways and precedes apoptosis, necrosis, and other cytostatic pathways [10]. Research on the significance of MC for pan-cancer is in full swing. Chen et al. found that integrating paclitaxel and mulberry water extract synergistically promotes MC in bladder cancer [11], and Lima et al. showed a landscape of MC in pancreatic cancer after receiving AD80 [12], which undoubtedly highlights the potential of MC as a chemotherapeutic target. In addition, MC has been reported to be inextricably linked to better survival in colon cancer [9], and prostate cancer [13]. We found that previous studies have constructed reliable prognostic models based on MC-related genes in hepatocellular carcinoma [14], colon cancer [15], and breast cancer [16], and highlighted the impact of their heterogeneity on the immune microenvironment. But so far no stratified model of bladder cancer based on mitotic catastrophe has been developed.

In this study, we used bioinformatics methods as a powerful weapon to construct a robust prognostic model for BLCA patients based on mitotic catastrophe-related genes for the first time, which filled the gap of previous studies. Based on risk stratification, we explored the factors mediating the prognostic differences between high- and low-risk patients with BLCA, such as functional activation heterogeneity, differences in immune infiltration, and so on. Subsequently, we identified the core gene in our model, characterizing its key role identity in BLCA progression. In conclusion, our study identified many factors affecting the progression of BLCA and discovered new reliable target, which provide the basis and valuable insights for subsequent in-depth mechanistic experimental analysis as well as clinical management.

We adopted a whole-to-local analysis approach, first exploring the prognostic differences and influencing factors between model risk strata, and then focusing on the most influential core genes in the model and identifying their identity in BLCA progression. The main part of this study is presented in Fig. 1, the left half of the picture is the part of overall analysis, while the right side is the panel of core gene analysis.

Fig. 1

Key components of the research

Materials and methods

Data acquisition

TCGA-BLCA transcriptomic data as well as clinical phenotypes, mutation data were downloaded from UCSC Xena (https://xena.ucsc.edu), including 400 cancer samples as well as 18 normal samples. Cancer patient-specific clinical survival data is maintained in Supplementary Table 1. Three training sets (GSE13507, GSE32894, GSE32684) as well as single-cell RNA-seq (GSE135337, GSE130001, GSE145281) derived from Gene Expression Omnibus (GEO), whereas “IMvigor210 CoreBiologies” R package was utilized to obtain immunotherapy cohort information. Specifically, the GSE13507 (platform: GPL6102, country: South Korea) cohort contained 165 samples of patients with primary BLCA, the GSE32894 (platform: GPL6947, country: Sweden) cohort contained 224 patients with uroepithelial carcinoma with complete survival information, and the GSE31684 (platform: GPL570, country: USA) cohort included information on 93 patients with BLCA. Specific clinical survival information of these GEO cohorts is provided in Supplementary Table 2–4. In addition, the Imvigor210 cohort consisted of 348 BLCA patients treated with anti-PD-L1 immunotherapy. GSE135337 (platform: 10X Genomics, Illumina-NovaSeq 6000, country: China) contains single-cell sequencing data from seven BLCA samples of different tumor stages, with specific sample parameters listed in Supplementary Table 5. GSE130001 (platform: 10X Genomics, Illumina-HiSeq 2500, country: USA) contains single-cell sequencing data from two muscle-invasive bladder cancer samples, while GSE145281 (platform: 10X Genomics, Illumina-HiSeq 4000, country: USA) contains single-cell sequencing data from 10 BLCA samples treated with anti-PD-L1 immunotherapy. All three samples were used to explore the expression patterns of the core genes in subsequent analyses in different cell types of BLCA. In particular, we analyzed core gene-mediated functional regulation in GSE135337 and GSE130001. Furthermore, mitotic catastrophe-related genes (MCRGs) were consistent with previous literature [15, 16].

Identification of differentially expressed genes (DEGs)

Throughout the study, the “limma” R package was utilized for difference binning of bulk RNA-seq, with the thresholds being an absolute value of log2 FC greater than 1 and an adjusted P-value less than 0.05. The analysis of DEGs of BLCA versus normal bladder tissue was based on the TCGA-BLCA cohort, whereas the DEGs used for the enrichment analysis were obtained from the differential expression analysis between the high and low risk groups of the TCGA-BLCA cohort.

Identification of core modules for DEGs

In order to condense the parts of bladder cancer and normal tissue DEGs that are most relevant to the tumor phenotype, we performed weighted gene co-expression network analysis (WGCNA) [17] based on TCGA-BLCA cohort. We constructed expression networks with maximally scale-free distributions based on optimal soft thresholds and extracted modules of co-expressed genes using the dynamic tree-cutting algorithm, and the modules were analyzed in terms of Pearson's correlation with the phenotype linkage. Referring to previous literature [18], we extracted the most positively and most negatively correlated gene modules with Tumor phenotypes to ensure that the core gene modules were important in promoting and suppressing Tumor phenotypes, and then combined them for subsequent analyses.

Optimal prognostic modeling

After taking the intersection of the WGCNA results with the MCRGs to obtain the fraction of significant dysregulation, univariate cox regression was used to identify the fraction of them that was significantly associated with prognosis. Based on these genes, we integrated 10 machine learning algorithms for sorting and modeling [19]. Details about the specific algorithms involved in modeling are described in Supplementary Methods 1.1. We used TCGA-BLCA as the training set for modeling, while three GEO cohorts, GSE13507, GSE32894, and GSE31684 were used as independent validation. The contents of the three GEO data were not merged with the aim of exploring the generalization ability of the model in different data environments. The modeling approach with the largest average C-index was identified as the best algorithm in the four datasets including the training set. High- and low-risk groups of each cohort were divided according to the median score of the training set, and K-M survival curves as well as ROC curves were used as the assessment method of modeling effectiveness. Subsequently, the independence of the model score as a prognostic factor for BLCA patients was corroborated by double cox regression.

Enrichment analysis

Based on the DEGs between high and low risk groups, the “clusterProfiler” package was utilized to perform enrichment analysis based on GO and KEGG gene sets. In addition, we calculated the Hallmark pathway scores for each sample using the “GSVA” package, and presented the more significant differences between the high- and low-risk groups in the form of heatmaps. Finally, we performed Gene Set Enrichment Analysis (GSEA) tailored to GO and KEGG datasets to show the functional enrichment landscapes of different risk groups from another perspective, which lays the foundation for a deeper understanding of the potential patterns of mitotic catastrophe heterogeneity.

Immunity landscape profiling

We evaluated the abundance of 22 immune cells in each sample by the well-known Cibersort algorithm [20] and presented the risk score and its correlation by the “ggcor” package. Subsequently, utilizing the ssGSEA method, we determined the immunocompetence of the samples based on a set of 29 immune function marker genes. In addition, we demonstrated intergroup differences in the expression of HLA-related genes and the correlation between risk score and the expression of immune checkpoint genes. Finally, we explore the immunotherapy guidance of risk stratification in the Imvigor210 cohort.

Malignancy and mutation analysis

By ssGSEA method, we evaluated the hypoxia score, epithelial-mesenchymal transition (EMT) score and angiogenesis score of each sample based on Hallmark-derived gene sets. The “RCircos” package was utilized to map the localization of the genes on the chromosomes. The mutation data of the TCGA samples were processed by the “maftools” package to calculate the TMB of each sample, and the overall mutation landscape is presented in a beautiful waterfall plot.

Precise extraction of core gene

In order to increase the reliability of core gene extraction, we integrated the cohort data by de-batching using the “sva” package and plotted the effect before and after processing by principal components analysis. We believe that a crucial feature of core genes is that they are not interfered by other genes, so multivariate cox regression was utilized for initial screening. Subsequently, we creatively constructed the SurvSHAP(t) model to analyze the time-based prognostic significance of variables [21, 22]. This model is based on Shapley Additive exPlanations (SHAP), which explains the importance of local features in black-box survival models by aggregating the effects of time-dependent variables [23, 24], and is provided by the “SurvSHAP” package. The most profoundly important genes in the model were selected as the core of the study for subsequent in-depth analysis.

Mendelian randomization analysis between core gene and BLCA

To explore the causal association of the core gene with BLCA, “TwoSampleMR” package was utilized for MR analysis. The eQTL data for the core genes were obtained from the IEU OpenGWAS database, and the BLCA ending data were provided by the FinnGen database. SNPs included in the analysis met a strong correlation with the exposure factor (P < 5e-04), excluded linkage disequilibrium (thresholds of r² < 0.001, kb = 10,000) and had an F-value greater than 10. Inverse variance weighting method was considered the primary decision-making tool for determining causal associations.

qRT-PCR and immunohistochemistry

Details of cell culture are detailed in Supplementary Methods 1.2. Our previous study described specific methods for total RNA extraction and qRT-PCR [25]. The following primers were used in this study:

ANLN forward, 5'-ATCTTGCTGCAACTATTTGCTCC-3'and.

ANLN reverse, 5'-TCCTGCTTAACACTGCTGCTA-3';

Actin forward, 5'-GAAGATCAAGATCATTGCTCCTC-3'and.

Actin reverse, 5'-ATCCACATCTGCTGGAAGG-3’.

Actin was used as an internal control. The procedures were performed three times to ensure accuracy and precision. The relative expression levels were calculated using the 2^−ΔΔCT method.

Immunohistochemical images for the gene were sourced from The Human Protein Atlas, comparing both cancerous and non-cancerous bladder tissues. The bladder tissue number in the database is T-74000. Specifically, we have selected images of staining for antibody CAB062547. Normal bladder tissue was taken from patient number 1938, while sections of high-grade bladder cancer tissue originated from patient number 3079.

Single-cell analysis

Analysis of single-cell profiles of GSE145281 and GSE130001 was based on TISCH (http://tisch.compbio.cn/home/). As for the processing of GSE135337, we relied on the Seurat V4 package under R software. After reading seven BLCA samples, constructing the analyzed objects by CreateSeuratObject function and integrating them, we converted ENSEMBL ID to Gene Symbol. After strict quality control with the criteria of the number of characterized genes between 500 and 6000, and the proportion of mitochondrial genes was less than 10%, we normalized, log-transformed and scaled the samples utilizing the NormalizeData, FindVariableFeatures and ScaleData functions. Subsequently, principal components were determined by PCA and batch effects were removed by RunHarmony. Finally, dimensionality reduction clustering is presented by UMAP with a resolution parameter of 0.5. Annotation markers for cell population types were collected from previous literature [26]. In addition, the extraction of DEGs at the single-cell level was based on the FindAllMarkers function with the threshold set to log2 FC absolute value greater than 0.25 and adjusted P-value less than 0.05. In our analysis, we explored the expression patterns of core gene in epithelial cells by GSE135337 and GSE130001. Instead of combining the two, we first analyzed in GSE135337 and then in GSE130001, in order to realize the “discovery-validation” process and fully explore the reliability of the expression pattern of the target gene.

PPI network construction, ICD score calculation and methylation analysis

The core gene-based protein–protein interaction (PPI) network was obtained from GeneMANIA (https://genemania.org/). Immunogenic cell death scores were obtained by the ssGSEA algorithm of the GSVA package. ICD-related genes were obtained from previous literature [27]. In addition, the methylation-related genes used in the study were derived from previous literature [28]. The MethSurv website (https://biit.cs.ut.ee/methsurv/) was utilized to explore the association of methylation of core gene with BLCA prognosis.

ceRNA regulatory axis construction

We searched four databases (Supplementary Methods 1.3) for miRNAs targeting core gene and took the intersection as the most likely upstream regulatory miRNAs. miRNA binding site prediction to core genes was by RNAhybrid [29]. lncRNAs that interact with miRNAs were predicted from the ENCORI database. The spearman method was adopted for all correlation analyses.

Drug prediction and molecular docking

Drugs with core genes as targets were predicted from the DSigDB database in Enrichr (https://maayanlab.cloud/Enrichr/). Subsequently, the structure files of the drug molecules were downloaded from the PubChem database and the structures of the core gene proteins were obtained from the Protein Data Bank. After removal of water molecules, addition of hydrogen atoms and calculation of charges, molecular docking was carried out at protein eutectic ligands by AutoDock 4.2 as well as AutoDock tools 1.5.7, followed by visualization by PyMOL 3.0.5.

Statistical analysis

Data analysis and graphical plotting for the entire study was done using R software of version 4.3.3. P-value less than 0.05 was considered significant while adjusted P-value was obtained by Benjamini-Hochberg (BH) method. Wilcoxon test was used for comparison of differences between the two groups, if non-specific emphasized, correlation test was taken by Spearman method.

Results

Identification of prognosis-associated MCRGs

Based on the expression profiles of the TCGA-BLCA cohort, the limma package was utilized for differential expression analysis, resulting in the extraction of 5370 DEGs. Subsequently, in order to identify the parts that are more important for tumor characterization, we took a WGCNA approach. A soft threshold of 7 was chosen (Fig. S1A) and we distinguished the differential genes into 6 modules (Fig. S1B-C). Among them, the blue and black modules were of interest to us for presenting the most significant positive and negative correlations with tumors (Fig. S1D-E). We then integrated the genes of the two modules and took the intersection with the MCRGs to obtain 206 dysregulated MCRGs (Fig. S1F). In order to visualize the differential distribution of the MCRGs, delicate volcano and heat maps were used to demonstrate (Fig. S1G-H). Finally univariate cox was utilized to identify the components that were strongly associated with prognosis (Fig. S2). Overall, close to 30% (246/900, 27.33%) of mitotic catastrophe-associated genes were deregulated in BLCA, and of the 246 genes that were deregulated, more than eighty percent (206/246, 83.74%) were found to be closely associated with the tumor phenotype, suggesting that mitotic catastrophe-associated processes in the progression of BLCA lesions may have some important biological significance. Finally, we screened the fraction of them associated with BLCA prognosis by univariate cox regression for subsequent in-depth analysis.

Construction of machine learning prognostic models

Through a combination of machine learning algorithms, we ultimately identified 98 models for prognostic assessment. The CoxBoost joint RSF modeling method was identified as the best one for its highest C-index (Fig. 2A), encompassing 16 prognostic genes: FBN2, MAP1A, IGF1, JAK2, SBSN, TXNIP, LRP1, LDLR, NGFR, GEN1, ANLN, PSME2, FANCD2, CHEK2, ELN, GMNN. Based on the median risk score, the TCGA cohort was categorized into high and low risk subgroups (Fig. 2B). Interestingly, scatter plot showed a clustering of patients who died and a decrease in overall patient survival time as the score increased (Fig. 2C). Subsequent KM curve analysis highlighted significant differences in OS, DFS, and PFS between high- and low-risk subgroups (Fig. 2D-F), demonstrating the model's excellent stratification ability. In addition, the area under the time-dependent ROC was higher than 0.95 (Fig. 2G), which was a near-perfect predictor of patient survival years. Despite the possible suspicion of overfitting, our model also showed remarkable prognostic stratification in the three validation sets. Dividing the validation set according to the median score of the training set, we found equally significant differences in OS between high and low risk groups (Fig. 2H-J), which was encouraging. In GSE13507, the AUCs of 1-, 3-, and 5-year ROCs were 0.708, 0.722, and 0.689, respectively, whereas they were all higher than 0.75 or even 0.8 in GSE32894 cohort and near 0.6 in GSE32684 (Fig. 2K-M). In conclusion, the mitotic catastrophe-related prognostic model we constructed had excellent risk stratification capabilities.

Fig. 2

Construction of prognostic model. (A) C-index heatmap of machine learning algorithms (B) Scatterplot of median stratification (C) Risk score-survival time scatterplot (D-F) KM curves of OS, DFS, PFS in TCGA cohort (G) ROC curves in TCGA cohort (H-J) KM curves of OS in GEO cohorts (K-M) ROC curves in GEO cohorts

Independent prognostic verification

Utilizing pie charts we organized and presented the distribution of clinical characteristics of the TCGA cohort, and unsurprisingly, the proportion of progressive clinical features was higher in the high-risk subgroup (Fig. 3A). Similarly, risk scores increased as the clinical grade progressed, although interestingly, female patients demonstrated higher risk scores (Fig. 3B). We subsequently confirmed by KM curves that the model demonstrated significant stratification ability across groups with different clinical characteristics (Fig. 3C), and were further illustrated by cox analyses that risk classification was an independent prognostic factor in each cohort (Fig. 3D).

Fig. 3

Independent prognosis for risk stratification. (A) Pie charts of the distribution of different clinical characteristics in the high and low risk groups (B) Comparison of risk scores in the context of various clinical characteristics (C) KM survival curves in the context of various clinical characteristics (D) Forest plot of cox regression for risk stratification versus clinical characteristics in each cohort

Functional analysis of differences between groups

After the DEGs were extracted between high and low risk groups based on established thresholds set in the Methods section, we launched an in-depth exploration of functional activation heterogeneity, starting from enrichment analysis. GO-based enrichment analysis revealed that extracellular matrix-associated functions were significantly probed, focusing on the organization as well as the synthesis of the extracellular matrix (Fig. 4A-B). The functional description corresponding to the GO ID in Fig. 4B, we present in Table 1. In addition, we specifically analyzed the differences of above functions between high and low risk groups. Scatter plots presented a positive trend between the above functions and Riskscores, and boxplots further emphasized that these functions were scored higher in the high-risk group. Moreover, the expression of squamous BLCA markers from previous literature [57] showed a similar trend (Fig. S3). Subsequent KEGG enrichment suggests dysregulation of many metabolic pathways, as well as pro-oncogenic pathways such as PI3 K-AKT and Wnt (Fig. 4C). Hallmark-based GSVA then more clearly mapped the upregulated pro-cancer pathways in the high-risk group (Fig. 4D). Finally, in-depth GSEA further corroborated that the function of extracellular matrix was active in the high-risk group, and the activation of Focal adhesion could not be ignored (Fig. 4E–F). In this section, we explored differences in functional landscapes between patients in the mitotic catastrophe-related high- and low-risk groups, revealing possible factors mediating a poorer prognosis and providing direction for an in-depth exploration of disease progression in BLCA.

Fig. 4

Comprehensive enrichment analysis. (A) Bubble diagram of GO enrichment (B) Circle diagram of GO enrichment (C) Bubble diagram of KEGG enrichment (D) Hallmark-based GSVA heatmap (E–F) GSEA in the context of GO (E) and KEGG (F)

Table 1 Specific descriptions corresponding to GO IDs in Fig. 4 B

Immune microenvironment profiling

It is well known that the complexity of the immune microenvironment plays an important role in tumorigenesis and progression, and we thus meticulously analyzed the heterogeneity of immune infiltration patterns among patients based on the expression profiles of the TCGA-BLCA cohort. Initially, we demonstrated the extent of infiltration of 22 immune cells in all BLCA patients by bar plots (Fig. 5A). Comparison of abundance highlighted higher infiltration of CD8 T cells, memory B cells, and follicular helper T cells in the low-risk group and showed a significant negative correlation with the risk score, whereas the opposite was true for mast cells (Fig. 5B-C). Moreover, it is interesting to note that we found higher infiltration of neutrophils in the high-risk group, and the positive correlation between risk score and M2 macrophages is not negligible. The specific correlation coefficients and test p-values are presented in Fig. S4A. Immunofunctional analyses also highlighted upregulated mast cell pathways in the high-risk group (Fig. 5D). In addition, we found only four differentially expressed HLA molecules, all of which were upregulated in the high-risk group (Fig. 5E). These included three coherent nonclassical HLA molecules, HLA-E, F, and G, all of which were located within the class I region and, together with the classical class I antigens, constituted a complete list of active class I genes in humans. Interestingly, the correlation between risk score and the expression of immune checkpoints was complex, with TNFSF9, for example, showing a positive correlation and CD40, for example, the opposite (Fig. 5F). Combining the above heterogeneity of the immune microenvironment between the two risk groups, we further explored the significance of risk stratification as a guide in the Imvigor 210 cohort receiving immunotherapy. Unsurprisingly, overall survival was worse in the high-risk group (Fig. 5G) and patients without a therapeutic response demonstrated higher risk scores (Fig. 5H), reflecting the clinical guidance potential of our model. In this section, we revealed the differences in immune infiltration between high and low risk groups associated with mitotic catastrophe, shed light on factors mediating prognostic heterogeneity in terms of immune cellular dynamics, and emphasized the efficacy of risk scores as predictors of immunotherapy.

Fig. 5

Immunoinfiltration profiling. (A) Histogram of the percentage of immunocytes (B) Box plots comparing immune cell abundance between high and low risk groups (C) Correlation of risk score with immune cells (D) Box plots of differences in immune function between high and low risk groups (E) Differential expression of HLA-related genes (F) Correlation of risk score with immune checkpoint gene expression (G) KM curves for Imvigor 210 cohort (H) Comparison of Risk Score for Immunotherapy Response Subgroups

Malignancy and mutation analysis

To further explore the contributing factors to the risk increasing process, we ventured into the association of risk scores with three notoriously malignant indicators: hypoxia, EMT, and angiogenesis. Indisputably, we found that the three major metrics were upregulated in the high-risk group (Fig. 6A) and showed a positive correlation with the risk score (Fig. 6B), suggesting the importance of these biological pathways in the evolution of risk. Interestingly, the opposite was true for TMB, which was instead higher in the low-risk group (Fig. 6C) and significantly negatively correlated with the risk score (Fig. 6D). Survival analysis showed a better prognosis in the high TMB group (Fig. 6E), especially in patients with low risk accompanied by high TMB (Fig. 6F). The unique significance of the mutation burden caught our attention and prompted us to analyze the genetic changes in depth. We first demonstrated the chromosomal localization of 16 model genes (Fig. 6G), and CNV statistics showed that all but ELN were altered to varying degrees, with TXNIP being observed to be the most amplified (Fig. 6H). Interestingly, LRP1 exhibited the highest frequency of mutations among the model genes at about 7% (Fig. 6I). We then collated the overall mutation landscape, with missense mutations occurring far more frequently than others and C to T base transitions being the most frequent (Fig. S4B-C). And meticulously comparing risk subgroups, we found that ARID1 A had a significantly higher mutation rate in the high-risk group, and the opposite was true for KDM6 A and MUC16 (Fig. 6J-K). The mutation rate ranking of pathways did not differ much between high and low risk subgroups, except that it was more frequent in the low risk group (Fig. S4D-E). In conclusion, we found that risk scores showed a positive correlation with several malignant phenotypes and that there was somatic mutation heterogeneity between patients in the high and low risk groups.

Fig. 6

Malignancy analysis and mutational landscapes. (A) Differences in malignancy scores between high and low risk groups (B) Correlation of risk scores with malignancy scores (C) Differences in TMB between high and low risk groups (D) Correlation of risk scores with TMB (E) KM survival curves based on the median TMB (F) TMB-Risk stratified survival curves (G) Chromosomal localization of model genes (H) CNV of model genes (I) Mutation waterfall plot for model genes (J-K) Mutation waterfall plots for high (J) and low (K) risk groups

Characterization of the core model gene

Initially, we integrated the training set as well as the three validation sets by de-batching to expand the reliability and validity (Fig. S5A-B), and then screened them for deeper prognostic factors that were not affected by interactions through multivariate cox (Fig. S5C). Brilliantly, based on the results of the multivariate cox screening, we constructed a survshap(t) explanatory model to determine the extent to which each factor contributes to the outcome, ultimately highlighting ANLN as an excellent influencer (Fig. 7A-B). We then focused on ANLN, which was upregulated in expression in cancer samples at the overall sample as well as paired sample level (Fig. 7C-D). The AUC for diagnostic ROC was as high as 0.849 (Fig. 7E), highlighting the value of ANLN as a BLCA characterization factor. Interestingly, we found ANLN expression rising with progressive disease stage in the TCGA cohort (Fig. 7F-G). In addition, ANLN expression was also significantly higher in deceased samples than in surviving samples (Fig. 7H). Finally, we affirmed the potential of ANLN as a prognostic factor by survival analysis, with significant differences in OS, DSS, and PFS between the high and low expression groups in the TCGA cohort (Fig. 7I), and focused on the same OS stratification in the three GEO cohorts (Fig. 7J). In this panel, we identified ANLN as a core gene for the model, and its high expression predicted a poorer prognosis for BLCA patients.

Fig. 7

Screening for model core gene. (A) SurvSHAP (t) model feature importance ranking (B) Aggregated SurvSHAP (t) values (C) Differences in ANLN expression between tumor and normal samples (D) Differential ANLN expression in TCGA paired samples (E) Diagnostic ROC curves of ANLN (F-G) Differences in ANLN expression in patients with different Stages (F) and T-stages (G) (H) Differences in ANLN expression between dead and surviving patients (I) OS, DSS, PFS curves for ANLN stratification in TCGA cohort (J) OS curves for ANLN stratification in each GEO cohort

Mendelian randomization analysis between core gene and BLCA

Three SNPs of ANLN were screened for inclusion in the MR analysis, and all four methods, with the exception of the MR Egger method, demonstrated a positive correlation between BLCA and SNP effect on ANLN (Fig. 8A). The forest plot accurately portrayed the significance of exposure factors and outcomes, with ANLN being a significant risk factor in the Inverse Variance Weighted method (Fig. 8B). In addition, funnel plots identified the absence of bias in SNPs (Fig. 8C) and affirmed the risk factor status of ANLN for BLCA by leave-one-out analysis (Fig. 8D). On the other hand, inverse MR analysis outlined that the causal relationship of BLCA to ANLN was not statistically significant, thus firming our judgment of ANLN as an oncogene (Fig. 8E–H).

Fig. 8

MR analysis. (A) Impact of ANLN-related SNPs on BLCA outcomes (B) Forest plot of ANLN-related SNPs on BLCA outcomes (C) funnel plot (D) leave-one-out sensitivity test (E–H) Reverse Mendelian randomization with BLCA-associated SNPs as exposure and ANLN as outcome

Expression validation

Immunohistochemical profiles derived from the HPA database confirmed the high expression of ANLN in BLCA at the protein level (Fig. 9A), and in addition, we performed a more direct validation at the mRNA level by qRT-PCR. Unsurprisingly, the experimental results of cell lines as well as tissues were in agreement with the above (Fig. 9B-C).

Fig. 9

Exploration of expression modes. (A) Immunohistochemical profiles of ANLN (B) qRT-PCR of ANLN in cell lines (C) qRT-PCR of ANLN in tissues

Single-cell analysis

Further analysis at single-cell resolution provided us with a chic view of the ANLN expression landscape. Initially, we found that ANLN expression was not probed in the GSE145281 dataset composed entirely of immune cells (Fig S6 A-C). We then performed in-depth profiling in GSE135337 containing 7 BLCA samples. Clustering yielded 19 clusters (Fig. 10A), and based on the marker provided by the original authors, we defined the cells into five categories (Fig. 10B), and interestingly, epithelial cells accounted for the vast majority of the cells, which the original authors determined to be all of malignant phenotype by inferCNV. Subsequently, we found that ANLN expression showed distributional specificity, being significantly enriched in the cell cluster corresponding to cluster 5 (Fig. 10C). This cluster of epithelial cells was isolated from the bulk of the population, suggesting that it may be characterized by a different cellular grouping. We therefore subjected cluster 5 to differential expression analysis with the remaining epithelial cell clusters to explore intercluster heterogeneity, and ended up with 98 DEGs, of which only 5 were down-regulated in expression. We visualized the intercluster expression differences of these 98 DEGs by a heatmap (Fig. S6D). The heatmap shows that in addition to ANLN, recognized proliferation markers such as MKI67, CDK1, and PCNA are upregulated in cluster 5. Interestingly, the enrichment results in the context of Disease Ontology as well as DisGeNET showed that DEGs are involved in tumor progression, especially in urological tumors (Fig. 10D-E). GO and KEGG analyses consistently highlighted the involvement of DEGs in a number of important processes as well as structures of cytokinesis, such as nuclear division, DNA replication, cell cycle, and so on (Fig. 10F-G). This suggested to us that Cluster 5, which exhibits high ANLN expression, has a strong splitting activity, therefore it was identified by us its as a proliferative subgroup. And ANLN may play an indispensable role in mediating such divisive tendencies.

Fig. 10

Comprehensive single-cell analysis. (A) Cluster distribution of GSE135337 (B) Cell type annotation of GSE135337 (C) Expression landscape of ANLN in GSE135337 (D) Disease Ontology enrichment (E) DisGeNET enrichment (F) GO enrichment (G) KEGG enrichment (H) Clusters of GSE130001 (I) Cell type of GSE130001 (J) Expression landscape of ANLN in GSE130001 (K) Expression of ANLN in different clusters (L) Intercellular communication of C1 and C3 (M) GSEA between clusters (N) Display of E2 F_Targets score and G2M_checkpoint score

We then focused our attention on another BLCA dataset, GSE130001, to further analyze the identity of ANLN. After collating the cell cluster subgroups of the GSE130001 dataset (Fig. 10H-I), we found that ANLN expression was particularly significant in epithelial cells (Fig. 10J). Notably, cluster 5 epithelial cells deviated from the predominant epithelial clusters, suggesting some possible heterogeneity. We extracted DEGs of cluster 5 with the remaining clusters in GSE130001. KEGG enrichment analysis showed that DEGs were involved in cellular senescence as well as the P53 pathway (Fig. S6E). We visualized the P53 pathway score more visually and found that it was indeed significantly upregulated in cluster 5 (Fig. S6F, Fig. 10M). Also upregulated were TNF-α signaling via NF-κB and TGF-β signaling pathway, pathways involved in senescence-associateed secretory phenotypes (Fig. S6F, Fig. 10M). In addition, the important cellular senescence-related marker CDKN1 A (P21) also had the highest expression in cluster 5 (Fig. S6G-H). We therefore identified cluster 5 as a senescent subpopulation in epithelial cells. Subsequently, we found that epithelial cell clusters 1 and 3 exhibited the most vigorous ANLN expression (Fig. 10K). Intercellular communication analysis showed that both cell clusters communicated closely with cell clusters 5 and 12 (Fig. 10L). Interestingly, both C1 and C3 showed frequent MDK-SDC4/2 receptor-ligand pair cross-talk against C12, and they showed docking activity at MDK-SDC4/1 as well as MDK-(ITGA6 + ITGB1) against C5 (Fig. S7A-B). We then performed further cell cluster GSEA and interestingly, C1 and C3 exhibited strong E2 F_Targets, G2M_Checkpoint, and Mitotic_Spindle activation. C12, which communicates closely with these two cell clusters, exhibited activation of immune-related pathways such as Interferon_Alpha/Gamma_Response, TNFA_Signaling_Via_NFKB, whereas C5 exhibited activation of cancer-related pathways, such as P53, KRAS, etc., in addition to the above pathways (Fig. 10M). We then visualized the pathways upregulated in the ANLN-enriched region more intuitively (Fig. 10N).

Comprehensive functional analysis of ANLN

In-depth GSEA-based analysis revealed differences in pathway activation between high and low ANLN expression groups (Fig. 11A), providing insights into the functional mediation of ANLN. Similar to the above findings, we probed that the high ANLN group had significant activation of function at all temporal phases of mitosis. Interestingly, the low ANLN group, on the other hand, showed a more intense inactivation of metabolic pathways involving arachidonic acid, retinol, drugs, etc. The PPI network highlighted the physical interactions of ANLN with proteins such as MEN1, ZWILCH, and UBAP2L (Fig. 11B).

Fig. 11

Comprehensive functional analysis of ANLN. (A) GSEA between high and low expression groups (B) PPI network around ANLN (C) Differences in immune cell abundance between high and low expression groups (D) Correlation of Tregs abundance with ANLN (E) Differences in immune function between high and low expression groups (F) Correlation of partial immune function with ANLN (G) Survival curves stratified by ANLN in Imvigor 210 cohort (H) Differences in ANLN expression between immunotherapy effect groups (I) Differences in TMB between ANLN groups (J) Differences in ICD score between ANLN groups (K) Differences in methylation-related gene expression between high and low expression groups (L) Methylation-related gene expression correlates with ANLN

We then explored different degrees of immune activation between high and low expression groups to determine the identity of ANLN in the tumor microenvironment. Tregs were enriched in the low ANLN group, whereas the opposite was true for neutrophils (Fig. 11C). Furthermore, Tregs showed a significant negative correlation with ANLN. (Fig. 11D). Immunofunction, on the other hand, was quite different between the two groups, emphasizing that most of it was upregulated in the high ANLN group (Fig. 11E). Notably, inflammatory response showed a positive correlation with ANLN, coinciding with the Hallmark-based GSEA, highlighting its indispensable influence in the high expression group. In addition, interestingly, MHC class 1-related effects as well as CD8 T cell functions were also upregulated (Fig. 11F). Given that the function of MHC class 1 activated as well as toxic T cells is often considered critical for immunotherapy, we explored the potential of ANLN as a predictor of immunotherapy.The KM curves highlighted a better prognosis in the high expression group (Fig. 11G), which was corroborated by the higher expression of ANLN in the treatment-responsive group (Fig. 11H). Given the variation in immunotherapy outcomes, we further analyzed the association of several immunotherapy influences with ANLN expression. Interestingly, both TMB and immunogenic cell death (ICD) score were higher in the high ANLN expression group (Fig. 11I-J).

Methylation plays a complex role in tumor progression, and the methylation picture associated with ANLN is unclear. Interestingly, we found that the ANLN high-expression group showed a striking consistency with significant upregulation of methylation-related gene expression (Fig. 11K). And the correlation between the majority of methylation-related genes and ANLN was greater than 0.3, even as high as 0.65 (Fig. 11L), reflecting that ANLN may play an important facilitator role in mediating methylation. In addition, we analyzed the effect of methylation of ANLN itself on the prognosis of BLCA, and found that none of the three loci had a significant role (Fig. S8A), so we excluded this direction for follow-up research.

In addition, we were inspired by the squamous tendency exhibited by the high-risk group and therefore ventured to explore the correlation between ANLN and markers such as keratinization. Interestingly, the correlation between genes such as the KRT6 family and ANLN was greater than 0.4, suggesting that ANLN plays an important role in mediating keratinization (Fig. S8B-C). This was also corroborated by the up-regulated squamous-related scores in the high ANLN group (Fig. S8D-E).

In this panel, we reveal the identity of ANLN in BLCA patients through a multidimensional analysis, highlighting its involvement in mitotic processes and its association with better immunotherapeutic responses in the immune microenvironment. In addition, ANLN also significantly regulated methylation as well as squamous carcinogenesis trends.

ceRNA regulatory system

We used ceRNA as an entry point to delve into the potential upstream regulatory mechanisms of ANLN. Through the intersection of four database targets (Fig. 12A), we extracted the most probable upstream miRNAs of ANLN, among which hsa-miR-15a-5p received our attention due to its significant negative correlation with ANLN (Fig. 12B) and obvious down-regulation in cancer tissues (Fig. 12C). Interestingly, the KM curve showed that patients with high expression of hsa-miR-15a-5p had a better prognosis (Fig. 12D), and its minimal free energy for binding to ANLN was as high as −27.0 kcal/mol (Fig. 12E), which indicated that it had a high potential to inhibit the oncogenic factor ANLN. We then collected lncRNAs closely linked to hsa-miR-15a-5p in the ENCORI database and targeted MIR4435-2HG by correlation analysis, which was significantly negatively correlated with hsa-miR-15a-5p (Fig. 12F), whereas the positive correlation coefficient with ANLN was as high as 0.439 (Fig. 12G). Interestingly, its expression in tumor samples was very high (Fig. 12H) and was associated with poorer prognosis (Fig. 12I). Previous studies have shown that MIR4435-2HG can mediate the progression of head and neck cancer by promoting EMT, and we therefore ventured to explore whether a similar possibility exists in BLCA. Unsurprisingly, it exhibits a very strong positive correlation with EMT scores, with a coefficient of 0.663 (Fig. 12J). Specifically, it was negatively correlated with E-cadherin and positively correlated with three other EMT markers (Fig. 12K), including N-cadherin, suggesting the promotion of EMT is also a potential mechanism for its regulation in BLCA. Finally, we present the above conclusions in an overall model diagram (Fig. 12L), illuminating the direction of subsequent experiments on the mechanism.

Fig. 12

ceRNA network. (A) Identification of potential upstream miRNAs (B) Correlation between hsa-mir-15a-5p and ANLN (C) Differential expression of hsa-mir-15a-5p (D) Stratified survival curves of hsa-mir-15a-5p (E) Binding site of hsa-mir-15a-5p with ANLN (F) Correlation of MIR4435-2HG with hsa-mir-15a-5p (G) Correlation of MIR4435-2HG with ANLN (H) Differential expression of MIR4435-2HG (I) Stratified survival curves of MIR4435-2HG (J) Correlation of MIR4435-2HG with EMT score (K) Correlation of MIR4435-2HG with EMT markers (L) Model diagram of the regulatory network

Drug prediction and molecular docking

In the DSigDB database, we selected five drugs with possible effects with ANLN based on the adjusted P-value (Fig. 13A), and subsequently verified the credibility of drug-target binding by molecular docking experiments (Fig. 13B-F). Surprisingly, Phytoestrogens and Scriptaid showed remarkable affinities with ANLN as high as −9.7 and − 8.2 kcal/mol, respectively, whereas the binding potential of MeIQx as well as Diuron to ANLN is not negligible. The development of new drugs targeting ANLN can help to promote the breakthrough of clinical precision therapy, and our work points out the direction for further research in the future.

Fig. 13

Drug prediction and molecular docking. (A) Potential drugs (B-F) Molecular docking site demonstration

Discussion

Over the past years, the management of malignant invasive bladder cancer, especially the muscle-invasive type, has become more and more refined based on continuous breakthroughs in molecular characterization [30]. Recent studies have constructed a number of prognostic frameworks for bladder cancer, such as those based on glycolytic response [31], drug ADME-related genes [32], nicotine metabolism [33], and more. The prognostic potential of mitotic catastrophe-associated genes has been validated in breast cancer, hepatocellular carcinoma as described in the introduction, but the prognostic significance of mitotic catastrophe heterogeneity in bladder cancer is unclear. Our study is dedicated to filling this gap, providing insights into the evolution of BLCA progression, offering new options for the identification of a precise prognosis for BLCA as well as improving clinical management.

Our study identified prognostically relevant MCRGs in dysregulated genes by univariate cox regression as well as the CoxBoost algorithm, and subsequently constructed prognostic stratification models by the RSF algorithm. The 16 genes included in the model construction were FBN2, MAP1A, IGF1, JAK2, SBSN, TXNIP, LRP1, LDLR, NGFR, GEN1, ANLN, PSME2, FANCD2, CHEK2, ELN, GMNN. FBN2 encodes an elastin-related substance whose aberrant methylation plays an important role in the malignant invasion of lung cancer [34] and may serve as a prognostic factor in bladder cancer [35]. The microtubule protein MAP1A has been more extensively studied as a neurological disease target, such as Parkinson's disease [36]. It has also been reported as an autophagy-associated prognostic factor in BLCA [37], which is associated with microenvironmental immune infiltration and builds a ceRNA regulatory network [38]. Overexpression of Insulin-like growth factor-1 (IGF-1) promotes stemness in esophageal squamous carcinoma cells [39] and mediates the progression of hepatocellular carcinoma [40], while JAK2 is notorious for constituting a well-known pro-carcinogenic pathway with STAT3 [41]. SBSN enhances SRC-STAT3 signaling to stimulate BLCA metastasis by potentiating EGFR kinase activity [42], whereas TXNIP plays the role of an anti-cancer factor by inhibiting ERK kinase [43]. LRP1 mediates the DLL4-Notch pathway and promotes M2 macrophage differentiation thereby inhibiting immunotherapy efficacy [44], and, interestingly, another model gene, NGFR, mediates immunosuppression by inhibiting T cell sensitivity [45]. The role of LDLR and GEN1 in bladder cancer is unclear, but are important malignant factors in breast cancer [46, 47]. High expression of PSME2 predicts a higher degree of immune infiltration with a focus on M1 macrophages, which is widely recognized in pan-cancer [48, 49]. FANCD2 is highly expressed in tumor cells as a pro-oncogene and is associated with a poorer prognosis [50], but downregulation of its expression dramatically adds sensitivity to many chemotherapeutic agents [51,52,53]. While ELN is included in the management of bladder cancer prognosis as an inflammation-related feature [54], GMNN is complicit in the malignant progression of adrenocortical carcinoma [55]. Anillin was subsequently fully analyzed as the core of the model, which we will detail later. In summary, the 16 genes involved in the model construction highlight the prognostic potential from all aspects and the robustness of our model.

GO-based gene enrichment highlighted the upregulation of extracellular matrix (ECM) shaping-related functions, as evidenced by GSVA scores as well as GSEA. The extracellular matrix, as a major structure of the tumor microenvironment in solid tumors, encompasses a large number of glycoproteins, receptors, and even cytokines that provide fertile ground for cancer cell proliferation as well as proliferation [56]. Focusing on bladder cancer, we found that the ECM involved in the formation of the bladder structure at various layers showed different degrees of alterations, such as the up-regulation of the expression of fetuin, and thus a variety of enzymes related to matrix shaping of the ECM were considered as reliable biomarkers for BLCA [57]. On the other hand, dysregulation of ECM rigidity in BLCA affects smooth muscle function thus resulting in altered bladder function in patients [58]. Another interesting finding was that the high-risk group had a higher propensity for squamous carcinomatosis, as evidenced by activation of functions such as epidermis development and significant upregulation of keratinization-related marker genes. Previous molecular typing studies for BLCA have demonstrated the similarity of squamous carcinomatosis characterizing BLCA cells to undifferentiated basal-like cells, with a stronger tendency to proliferate, and the worst prognosis of all the typologies [59, 60]. The high proliferation rate accompanying squamous cell differentiation in the high-risk group may trigger more intense mitotic catastrophe-associated effects, and our study provides some novel insights into the progression of BLCA.

The profiling of the immune microenvironment presents a more subtle picture of the BLCA microcosm. In our analysis, CD8 T cells and follicular helper T cells were heavily infiltrated in the low-risk group and negatively correlated with risk scores. CD8 T cells, the mainstay of tumor toxicity killing, could predict a better prognosis and correlated with a stronger response to immunotherapy in bladder cancer [61]. Interestingly, follicular helper T-cells work both ways, promoting B-cell-mediated antibody responses and interacting with CD8 T-cells to enhance the immune effect, creating a synergistic effect on solid tumor killing [62]. There were previous studies showing that patients with BLCA lymph node metastases had more follicular helper T cells in lymph node drainage [63], so from another perspective, its high infiltration may also be a bystander in progressive BLCA. Thus its exact immune microenvironmental identity deserved further in-depth exploration. Furthermore we found dramatically that M2 macrophage abundance was significantly positively correlated with risk score. In addition to the well-known fact that M2 macrophages mediate immunosuppression as well as angiogenesis leading to a poorer prognosis in BLCA [64], there was evidence linking high M2 macrophage infiltration to basal-type BLCA [65], which certainly confirmed, on the other hand, that patients in the high-risk group may exhibit a more pronounced basal-like feature. Interestingly, the significant upregulation of mast cells in the high-risk group is also of interest. Mast cells have a complex identity in the microenvironment, where they were induced by chemokines to recruit immune cells after migrating to the vicinity of the tumor, but more evidence suggested that they mediate angiogenesis to promote metastasis and can reduce anti-tumor immunity in BLCA [66]. In addition, previous literature highlighted a higher density of mast cells in high-grade BLCA tissue [67]. And the urinary levels of some mast cell-related mediators, such as N-methylhistamine, also changed when treated with BCG [68]. In conclusion, mast cells are also a potential factor in the innovation of prognostic biomarkers for BLCA patients and in the study of BLCA progression. Another concern is the non-classical HLA molecules whose expression is upregulated in the high-risk group, mainly HLA-E, F and G. They have been reported to be important players in tumor immune escape as well as features such as drug resistance [69], which provides another perspective to explain the poorer prognosis. The complexity of the immune microenvironment leads to different efficacy of immunotherapy, and ultimately we found that the low-risk group has stronger immune response, which provides protocol support for precise clinical decision-making.

The analysis about mutations opens up new perspectives, and we find it interesting that the low-risk group exhibits a higher mutation burden. Previous studies have shown that the tubulointerstitial subtype BLCA has a higher genomic mutation rate, especially for KDM6 A [70], which coincides with our characterization of the low-risk group, suggesting that the low-risk group may consist of the tubulointerstitial subtype BLCA, which is significantly differentiated from the basal-like high-risk group. On the other hand, antigenic changes due to high mutations also contributed to the enhancement of the immune response [71], resulting in a rise in immunotherapy efficacy.

Freedom from reciprocal interference is in our view an important feature of the key core genes, so we performed a multivariate cox analysis screen after integrating the data from each cohort. Subsequently and wonderfully, we incorporated SurvSHAP(t) into the core gene screen to rank the time-dependent prognostic importance of genes. Although SHAP interpretation has grown considerably in the past, its use in bioinformatics analysis from a transcriptomic perspective is still relatively lacking, and our work further advances the integration of new methods in machine learning and provides a solution for deeper learning in genomics. ANLN, was then identified as the core of this model.

Anillin (ANLN) encodes an actin-binding protein that coordinates the interaction of actin with the cytoskeleton and plays an important role in mitosis and cell motility [72]. In addition, it has been reported to mediate malignant processes such as stemness upregulation and intercellular adhesion in breast cancer [73]. This is a notorious gene in pan-cancer and has also been recognized as a potential biomarker for bladder cancer [74]. In our study, Mendelian randomization analysis came first to elucidate the causality of ANLN as a risk factor for BLCA. IHC as well as qRT-PCR provided definitive evidence of upregulation in cancer cells.

Expression profiles under single-cell sequencing provide additional perspectives for in-depth studies. The first thing that deserves attention is the heterogeneity among cell clusters, which is manifested as the deviation of some clusters among the same cell type in the umap plot. In our study, cluster 5 in GSE135337 was identified as a proliferative subpopulation by its high expression of proliferation markers such as MKI67, CDK1, PCNA [75] and significant activation of mitosis-related functions. Previous studies have also shown that the proliferative subpopulation is segregated from other clusters in a single-cell perspective [76, 77], suggesting a distinct gene expression profile. In the case of GSE130001, the cluster 5 epithelial cell population was also segregated from the main epithelial cell population. We found that in addition to the well-known senescence pathway P53 [78] was upregulated in this cluster, TNF-α signaling via NF-κB and TGF-β signaling pathway were also significantly upregulated. These two pathways are closely related to the senescence-related secretion phenotype, and their up-regulation also implies an increase in the degree of cellular senescence [78]. In addition, CDKN1A (P21), a key senescence marker [78], also showed the highest expression among clusters in cluster 5. We therefore identified it as a senescent subpopulation. As also reflected in previous studies, clusters showing a clear senescence trend were separated from other clusters in the umap plot [79].

Interestingly, we found that ANLN was specifically expressed in some epithelial cells of BLCA and its high expression was found to be linked to more vigorous mitosis from a single-cell perspective, focusing on G2M checkpoint, Mitotic Spindle, and other links. We found that DEGs between clusters highly expressing ANLN and the other clusters were significantly enriched to many cancers, including urologic cancers. Moreover, these DEGs were also enriched to processes in all phases of mitosis. Given that cancers are characterized by exuberant cell division, it is understandable why these DEGs are enriched in various cancer phenotypes. And in addition, proteins that closely interact with ANLN also focused on cytokinesis functions, such as ZWILCH [80]. In our analysis, the MDK-SDC4 pathway, which was found to mediate major intercellular communication, was shown to play a key role in regulating cancer cell proliferation and migration [81]. Furthermore, activation of this pathway after chemotherapy portends activation of extracellular matrix remodeling functions, which triggers malignant manifestations, such as distant spread [82]. In our single-cell analysis, the number of ANLN-positive cells was low. We suggest that this is related to the biological context of ANLN, which plays an important function in mitosis, and its up-regulated expression in epithelial cells is usually associated with exuberant mitotic activity. Thus, ANLN is not a housekeeping gene that is significantly expressed in pan-cellularity. And since cells in the divisional phase are, after all, a minority in single-cell sequencing, it makes sense that the number of ANLN-positive cells is low. Given that malignancy is marked by uncontrolled mitosis in cancer cells, we believe that ANLN has the potential to serve as a reliable biomarker for BLCA. This expression property of ANLN deserves more direct biological validation and deeper investigation. We plan to construct a cohort of BLCA patients in our hospital to verify the histologic differences in ANLN expression in high- and low-risk patients, including through immunohistochemical staining. In addition, we plan to utilize spatial transcriptome sequencing and other means to explore the specific expression localization of ANLN, and to deeply investigate the significance of ANLN expression patterns in the progression of BLCA. It is believed that these subsequent studies will better address the current limitations.

The identity of ANLN in the immune microenvironment was also meticulously dissected, and it is noteworthy that inflammation-related functions were upregulated with its expression, which was also reflected in the GSEA analysis. Inflammatory response is not only one of the triggers of BLCA, but also participates in mediating the proliferation and proliferation of BLCA cells [83], and the specific pro-inflammatory mechanism of ANLN as well as the targeted therapeutic research will be one of the main focuses of our subsequent studies. On the other hand, our study highlighted that ANLN could be used as an observable indicator of better prognosis for immunotherapy in BLCA patients. Interestingly, this ran counter to the risk factor status of ANLN in conventional BLCA cohorts. We analyzed this in depth and attributed the possible reasons for this phenomenon to the following four aspects. First, the reduced abundance of Tregs and active CD8 T-cell-associated function in the ANLN high-expression group are consistent with the recognized characteristics of “hot tumors”. There is no doubt that hot tumors are more responsive to immunotherapy [84]. Second, we observed a higher tumor mutation burden in the ANLN high-expression group, which is consistent with our previous discussion that more neoantigens resulting from mutations promote immune response effects [71]. Subsequently, given that high ANLN expression predicted more vigorous mitosis, we followed this as a breakthrough for in-depth analysis. Previous literature suggests that the Warburg effect is also stronger in tumor microenvironments with more divisive activity, and that this metabolic reprogramming is accompanied by a high production of components such as reactive oxygen species [85]. Reactive oxygen species are an important factor in triggering immunogenic cell death [86]. Consistent with this main line, our analysis revealed a higher immunogenic cell death score in the ANLN high expression group. And the occurrence of immunogenic cell death can stimulate T-cell immunity and enhance the cancer cell killing effect of immunotherapy [87]. Finally, we found that ANLN was significantly and positively correlated with the expression of many major BLCA basal-type markers as well as with squamous carcinomatosis scores, suggesting that BLCA cells with robust ANLN expression may have a more pronounced basal-type profile. In contrast, previous studies have emphasized that basal-type BLCA are more responsive to immunotherapy [70]. Taken together, these four points give us reason to strongly believe that high ANLN expression predicts better immunotherapy outcomes, and thus it makes sense that ANLN serves as a better prognostic marker in the Imvigor210 cohort treated with anti-PD-L1 therapy.

Dysregulation of m6 A modification plays an important role in the progression of BLCA [88], and we found that ANLN showed a strong positive correlation with the expression of the 20 m6 A genes acquired, which provides an idea for epigenetic correlation studies of BLCA. In addition, the inextricable positive association between ANLN and squamous features, and the possibility that ANLN may mediate malignant evolution by advancing keratinization of BLCA cells, provide novel insights into the biological functions of ANLN. Deeper, we proposed a regulatory network of MIR4435-2HG, hsa-miR-15a-5p, and ANLN axis and highlighted the possibility that MIR4435-2HG mediates EMT, revealing possible upstream and downstream regulatory mechanisms. Ultimately, we screened a series of highly binding targeted drugs, including Phytoestrogens, based on the important target of ANLN, advancing the progress of precision therapy and making excellent contributions to the development of new drugs for BLCA.

While the results of our research are encouraging, it is important to recognize the limitations that exist. The first is related to the reliability and validity of the prognostic model; our construction process relied on the content of publicly available datasets, which may have led to biases in sample selection and ethnographic characteristics. Secondly the lack of treatment-specific information in the dataset leads to limitations in the significance of therapeutic decisions regarding the heterogeneity of mitotic catastrophes, especially immunotherapy and chemotherapy. In addition regarding the function of the core gene ANLN, we have examined it on a large scale mainly through bioinformatics methods, which may be at variance with the actual situation. Moreover, in the single-cell analysis part of the core genes, we did not integrate multiple datasets for the validation of a sufficiently large number of cells, which may have some impact on the reliability of the results. Taking these limitations into account, we plan to incorporate more prospective cohorts to assess model efficacy in the future, and to obtain more BLCA tissue samples for single-cell sequencing so as to obtain more data to further affirm the reliability of the results. Furthermore, we will try to validate the functional landscape of ANLN in detail through specific and fulfilling cytological experiments, and animal experiments.

Conclusion

Our study is the first to construct a BLCA prognostic model based on mitotic catastrophe-associated genes, and its risk-stratification ability was validated in multiple datasets. Integrating multidimensional analyses of functional activation, immune infiltration, and somatic mutations, we provided an in-depth discussion of the heterogeneity affecting prognosis. Subsequently, the core gene ANLN, identified by SuvrSHAP(t), was recognized by us as a promising BLCA target. Finally, we screened potential drugs based on the core gene and demonstrated binding efficacy with molecular docking. In conclusion, our study provides valuable insights by pointing the way for subsequent experiments to deeply analyze the malignant progression of BLCA.