Based on 10 cuproptosis molecules (18) ( Figure 1A ), we sought to establish a cuproptosis-associated gene signature for BLCA prognostication using the following public databases ( Figure 1B ). (i) The TCGA cohort of BLCA (legacy) that includes 407 tumor samples and 19 bladder nontumorous tissues (31). Patient information, pathology/histology, transcriptome, mutation, and copy number variation (CNV) data were downloaded from https://gdc.cancer.gov/. Aneuploid score was from reference (32). Tumor mutation burden (TMB) was calculated using Rpackage TCGAmutations. (ii) GSE13507 (33, 34), GSE154261 (35), and GSE176307 (36) BC cohorts. The data in these cohorts were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). (iii) E-MTAB-4321 cohort. The data were downloaded from https://www.ebi.ac.uk. (iv) IMvigor210 cohort. The data were from IMvigor210CoreBiologies (37, 38). For RNA sequencing data above, gene expression levels were measured using Transcripts Per Kilobase Million (TPM) and log2(x + 1) transformed. For array results (determined by 4×44K v2 microarray kit), we determined transcript abundances using probe-set values; and when multiple probes targeted the same mRNAs, the probes with largest mean values were chosen and then standardized using “Limma” package (39). During RNA sequencing and array processing, we conducted a two-step filtering. First, those genes with undetectable expression in >75% of samples were discarded. Second, we further got rid of genes with expression median absolute deviation (MAD) ≤0.01 and at the bottom 25%. The present study did not contain experimental analyses directly from human samples and animals, and thus needed no ethics permission.

To identify cuproptosis-associated genes, we first conducted single sample gene set enrichment analysis (ssGSEA) to obtain the normalized enrichment score (NES) according to expression levels of 10 cuproptosis genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS and CDKN2A) in the TCGA BLCA cohort ( Figure 1A ) using Rpackage “GSVA” with kcdf=Caussian and method = ssgsea. Weighted gene co-expression network analysis (WGCNA) and Pearson’s correlation were then performed to make sample clustering (tree) followed by the construction of a unsigned scale-free network, adjacency matrix and the topological overlap matrix, which eventually formed different modules ( Figure 2A ). Briefly, Pearson’s correlation was used to make sample clustering trees and none of the TCGA samples were outliners (with height >220) ( Figure 2A left). To define the optimal soft threshold value, we set up 1:20 as a power value, and when the scale independence reached 0.9 while mean connectivity was <100, the soft threshold value 12 was obtained ( Figure 2A middle). Based on this soft threshold setting, we constructed a unsigned scale-free network, adjacency matrix and topological overlap matrix through which the number of genes in each module was defined (maxBlockSize = 6000 and minModuleSize = 50). The function “WGCNA::blockwiseModules” was employed to assign genes into appropriate modules ( Figure 2A right). The correlation between each module and cuproptosis ssGSEA-NES together with clinical variables (stage and grade) was then evaluated, and by setting correlation R >0.30, we acquired the yellow (R = 0.34 with 559 genes) and turquoise modules (R = 0.32 with 2525 genes) ( Figure 2A right). Further filtering out genes at bottom Rs in these two modules (MM correlation R > 0.5, and GS correlation R > 0.2, P < 0.05) reduced the gene numbers to 392 and 1 586 in yellow and turquoise module, respectively. We then made COX and LASSO regression analyses to determine effects of these genes on patient progression-free survival (PFS) and expression differences between tumors and NTs ( Figures 2B, C ). Finally, 11 genes were obtained as the cuproptosis-associated 11 gene signature, which we named CuAGS-11. These 11 genes include C18orf54, NEIL3, ANLN, AHCY, PSMG1, TTC5, SRPRB, XPOT, ZC3HAV1L, SLC25A15 and P3H4.

Based on expression levels of 11 genes above, the CuAGS-11 score in each sample was calculated using the following formula:

Score = Σ βi × RNAi, where βi is the coefficient of the i-th gene in multivariable Cox regression analysis, and RNAi is RNA abundance of gene i. The obtained score values were further standardized using the scale function. Patients were classified into the high- and low-risk groups using the median score as a cut-off point. Differences in survival (OS, PFS and RFS) and BCG or ICI treatment efficacy between the CuAGS-11 high- and low-risk groups were then compared.

Time-dependent ROC curves and area under curves (AUCs) were used to estimate the accuracy of identified survival predictors (CuAGS-11 model and stage) in BLCA patients and made using Rpackage “timeROC”. We performed Cox regression analysis to evaluate the effect of the CuAGS-11 score and clinical parameters on survival in the TCGA and GSE13507 BLCA cohorts and established a predictive nomogram by using independent survival predictors in both cohorts to predict 1-, 3-, and 5-year survival (OS and/or PFS). The model-based predictive survival time against the observed one was plotted using calibration curves. R package “regplot” was used to make nomograms and to assess their predicative ability.

Based on the median CuAGS-11 score values in BLCA cohorts, we categorized tumors into low- and high-risk groups. Reference gene signatures required for Kyoto Encyclopedia of Genes and Genomes (KEGG) and Hallmark analysis were downloaded from https://www.gsea-msigdb.org/gsea/index.jsp (h.all.v2022.1.Hs.symbols.gmt’ and ‘c2.cp.kegg.v2022.1.Hs.symbols.gmt’). Differences in KEGG and hallmark pathways between two risk groups were evaluated using GSEA (version 4.2.1). Adjusted P value <0.05 and FDR <0.25 were regarded as significantly over- or under-represented pathways.

Proliferation statuses were estimated using expression levels of Ki-67 mRNA and cell cycle scores, respectively. Cell cycle, stemness and EMT signature scores were calculated based on ssGSEA or as the median z-score of signature gene panels for each sample. These signatures are as follow: Cell Cycle: CDK2, CDK4, CDK6, BUB1B, CCNE1, POLQ, AURKA, KI-67 and CCNB2 (40). Stemness score was assessed according to the mRNA expression-based stemness developed by Malta et al. (41). EMT score was calculated based on the dbEMT signature from http://dbemt.bioinfo-minzhao.org/ (42).

TIDE score is evaluated according to myeloid-derived suppressor cell (MDSC), macrophage M2, T cell Dysfunction and Exclusion (43). TIDE score for BLCA cohorts treated with Atezolizumab was calculated online at http://tide.dfci.harvard.edu/. mRNA abundance was standardized with use of the all sample average expression as the normalization control prior to TIDE score calculations.

All statistical analyses in the present study were conducted by using R package version 4.0.5. We performed Wilcox and K-W sum tests to determine differences between two groups and among multi groups, respectively. Spearman’s Rank-Order Correlation coefficient was used to assess correlation coefficient R between two variables. Survival analyses were carried out by using log-rank test, and Kaplan–Meier survival curves for visualization of OS, PFS and RFS were done using “Survival” and “Survminer” packages. Univariable and multivariable Cox regression analyses were employed to measure effects (HR and 95% CI) of quantitative predictive parameters on OS, PFS or RFS. P values < 0.05 were considered statistically significant. FDR correction was also performed to measure a statistical significance (< 0.25) when needed.

We first evaluated 10 cuproptosis molecules as survival predictors but failed to set up a satisfied model in the TCGA BLCA cohort (data not shown) (31). By analyzing cuproptosis-correlated genes as described in Methods, we developed the cuproptosis-associated 11 gene signature (CuAGS-11). These 11 genes include C18orf54, NEIL3, ANLN, AHCY, PSMG1, TTC5, SRPRB, XPOT, ZC3HAV1L, SLC25A15 and P3H4. The expression of these 11 genes was significantly higher in BLCA tumors than in their NT counterparts ( Figure 2D ). Survival analyses unraveled the significant association of PFS with each of 11 factors when patients were categorized into high and low groups using a median expression value as the cutoff ( Figure 2E ). We then calculated CuAGS-11 score in each tumor, and divided patients into high- and low-risk groups using the median CuAGS-11 score value as a cut-off point. The CuAGS-11 score-based grouping of the TCGA BLCA tumors was significantly associated with patient age, gender, grade, stage and metastasis ( Table S1 ).

BLCAs are in general stratified into luminal and basal subtypes according to their featured gene expression signatures (44–48). The luminal subtype is overrepresented by urothelium differentiation markers and transcription factors, while the basal one is poorly differentiated (5, 44, 48). To examine a potential association between molecular and CuAGS-11 subtypes, we analyzed 233 BLCA tumors well characterized for their differentiation subtypes in the TCGA cohort. As shown in Figure 3A , the basal subtype was significantly enriched in the CuAGS-11 high-risk group (high- vs low-risk: P = 5.193E-07).

Because the basal BLCA subtype is enriched with cycling and stem- and/or mesenchymal-like cells (5, 44), we further determined proliferation, stemness and EMT markers in those tumors. For proliferation analyses, Ki-67 was first used as the specific biomarker, and the CuAGS-11 high-risk tumors expressed significantly higher levels of Ki-67 mRNA (high- vs low-risk: P = 1.35E-19) ( Figure 3B ). Then, cell cycle scores based on ssGSEA were evaluated and similar results were obtained (high- vs low-risk: P = 2.21E-29) ( Figure 3B ). BLCA stem cell and EMT phenotype analyses showed significantly higher stemness and EMT scores in CuAGS-11 high-risk tumors (high- vs low-risk: P = 9.16E-07 and 2.95E-06 for stemness and EMT, respectively) ( Figures 3C, D ). Consistent with these findings, the GSEA hallmark analysis revealed overrepresentation of G2M checkpoint, mitotic spindle, E2F targets, glycolysis, PIK3-AKT-MTOR signaling and among others in the CuAGS-11 high-risk tumors ( Figure 3E ). GSEA KEGG analysis showed that cell cycle, MTOR, ERBB2, basal transcription factor, TP53 and other pathways were highly enriched in the CuAGS-11 high-risk tumors ( Figure 3F ).

We then probed whether there were differences in genomic alterations between the CuAGS-11 high- and low-risk tumors. First, global genomic aberrations including aneuploidy and TMB were evaluated. Aneuploidy scores were significantly higher in CuAGS-11 high- than low-risk tumors (P = 0.048) ( Figure 3G ), while there was no difference in TMB ( Figure 3H ). Alterations of individual genes were then compared, and we observed a significantly higher frequency of TP53 gene aberrations in the CuAGS-11 high-risk tumors (high- vs low-risk: 55% vs 34%, P<0.05 and FDR<0.05) ( Figures 3I, J ).

Having established the CuAGS-11 score model in BLCAs, we then evaluated whether it had impacts on patient survival (OS and PFS) from the TCGA BLCA cohort as a discovery one ( Table S1 ). These patients were categorized into high- and low-risk groups using the median CuAGS-11 score value as a cut-off point. Patients in the high-risk group had significantly shorter OS (P = 0.0001) and PFS (P = 0.001), as assessed by the Kaplan-Meier analysis ( Figure 4A ). Univariable COX regression OS analyses of age, gender, stage, grade, and the CuAGS-11 model showed that age (>60 yrs), advanced stages and CuAGS-11 (high-risk) were associated with significantly shorter OS ( Figure 4B ), while all these three variables remained highly significant in the multivariable COX regression analysis ( Figure 4C ). We observed a similar association between patient age (>60 yrs), advanced stages, and CuAGS-11 score (high-risk) and shorter PFS in the univariable COX regression analysis ( Figure 4D ). Multivariable COX analyses unraveled that only stage and CuAGS-11 risk score were independent predictive variables for patient PFS ( Figure 4E ).

We further determined the impact of the CuAGS-11 score on patient survival in the GSE13507 BLCA cohort to validate the findings obtained from the TCGA BLCA patients above. The GSE13507 cohort included 165 patients with BLCA (33, 34), and there were only OS data available. The patient characteristics were summarized in Table S2 . As expected, patients in the CuAGS-11 high-risk group had significantly shorter OS (P = 0.0002) ( Figure 4F ). In univariable COX regression analyses, OS was significantly associated with age, grade, stage and CuAGS-11 score ( Figure 4G ), whereas age, stage and CuAGS-11 score served as independent prognostic factors, according to multivariable Cox regression analysis ( Figure 4H ).

We then conducted time-dependent ROC and AUC analyses to evaluate the predictive ability of the CuAGS-11 model in the TCGA and GSE13507 BLCA cohorts. For TCGA patients, AUCs for 1, 3 and 5 year PFS by the CuAGS-11 model were 0.669, 0.634 and 0.674, respectively ( Figure 5A left panel). Like CuAGS-11 model, the stage was also an independent prognostic factor for OS and/or PFS in both cohorts, and moreover, the BLCA stage was a well-established predictor for long-term survival. Thus, we made a comparison of 5-year PFS prediction between CuAGS-11 model and stage. A slightly bigger AUC was observed for the CuAGS-11 model ( Figure 5A middle panel). We further combined CuAGS-11 with stage together and resulting AUCs increased substantially in predicting all PFS time points ( Figure 5A right panel). We then conducted the same analyses for OS in both TCGA and GSE13507 BLCA cohorts, and largely similar results were obtained ( Figures 5B, C ). Accordingly, we further established prognostic nomograms composed of CuAGS-11 score and stage. In the TCGA cohorts, the CuAGS-11/stage nomogram almost precisely predicted the possibility of both PFS and OS at 1, 3 and 5 years ( Figures 5D, E ). A highly accurate prediction of OS by the CuAGS-11/stage nomogram was observed in the GSE13507 cohort, too ( Figure 5F ).

The clinic-pathological variables have been mainly used to evaluate response to BCG therapy (8). We sought to determine whether the CuAGS-11 score could serve as such a predictive biomarker. The GSE154261 BLCA cohort (35) analysis of 73 BCG-treated patients showed that the recurrence and non-recurrence rates were 58.3% and 41.7% in the CuAGS-11 high-risk group, while 27.0% and 73.0% in the low-risk one, respectively (P = 0.009) ( Figure 6A top) ( Table S3 ). In the low-risk group, all the patients displayed stable disease status, whereas 45% of the high-risk patients underwent progression (P = 0.027) ( Figure 6A bottom). Consistently, the recurrence-free survival (RFS) and PFS were both significantly shorter in the CuAGS-11 high-risk group (P = 0.003 and 0.002, respectively) ( Figure 6B top and bottom). We further analyzed the E-MTAB-4321 cohort of 88 T1 BLCA patients who received BCG therapy (49) ( Table S4 ). None of 44 CuAGS-11 low-risk patients had disease progression, while 4 in the high-risk group exhibited progressive disease, although the difference was not statistically significant (P = 0.116) ( Figure 6C ). Nevertheless, PFS was significantly shorter in the CuAGS-11 high-risk patients (P = 0.042) ( Figure 6C ). Finally, 37 patients with MIBC in the TCGA cohort (31) were treated with BCG, and the recurrence rates were 72.2% and 36.8% for the CuAGS-11 high- and low-risk groups, respectively (P = 0.049) ( Figure 6D ) ( Table S5 ). OS was significantly shorter in the CuAGS-11 high-risk group (P = 0.040), whereas PFS was also shorter in this group but did not reach a statistical level (P = 0.092) ( Figure 6E ).

ICI therapy has been applied for MIBCs with good efficacy in subsets of patients (3). We further assessed whether the CuAGS-11 score was able to predict patient response to ICIs. For this purpose, IMvigor210 and GSE176307 cohorts were analyzed as the discovery and validation sets, respectively. A total of 348 patients received Atezolizumab (anti-PD-L1 antibody) therapy and 298 of them had response information available in IMvigor210 cohorts (37, 38) ( Table S6 ). These 298 patients were divided into the CuAGS-11 high- and low-risk groups based on the median score value. Patient responses to Atezolizumab were categorized into complete remission (CR), partial remission (PR), stable disease (SD) and progressive disease (PD), respectively (38). CR, PR, SD and PD were 2.7%, 8.0%, 23.5% and 65.8% in the high-risk group, whereas 14.1%, 20.8%, 18.8% and 46.3% in the low-risk group, respectively (P = 7.018E-06) ( Figure 7A ). OS information was available in this cohort and Kaplan-Meier analysis showed a dramatically shortened OS in the CuAGS-11 high-risk group (P = 9.20E-09) ( Figure 7A ). The GSE176307 cohort of 34 BLCA patients treated with Atezolizumab (36) was then analyzed for validation ( Table S7 ). All 17 patients in the high-risk group underwent disease progression, while more than half of patients in the low-risk group acquired CR (29.4%) or PR (23.5%) and only 35.3% of them had BC progression (high- vs low-risk, P = 8.65E-05) ( Figure 7B ). In accordance with their response rates, significantly longer OS and PFS were observed in the CuAGS-11 low-risk group (P = 8.12E-06 and 4.04E-03 for OS and PFS, respectively) ( Figure 7B ).

Given these observations, we further probed potential differences in infiltrated stroma and immune cells between the CuAGS-11 high- and low-risk tumors. To this end, we compared their TIDE scores. TIDE has been shown to predict ICI responses and determine mechanisms underlying tumor immune evasion (43). In the IMvigor210 cohort, the CuAGS-11 high-risk tumors displayed significantly higher TIDE score than did the low-risk tumors (P = 1.95E-05), and more specifically, robustly higher T cell exclusion score was observed in CuAGS-11 high-risk tumors (high- vs low-risk: P = 1.96E-05) ( Figure 7C ). The GSE176307 cohort analysis showed similar score differences in TIDE (P = 0.008) and T cell exclusion (P = 0.008) between high- and low-risk groups ( Figure 7D ).