Abstract
Bladder cancer (BC) or carcinoma (BLCA) is predominantly derived from urothelium and includes non-muscle invasive BC (NMIBC) and muscle invasive BC (MIBC). Bacillus Calmette-Guerin (BCG) has long been applied for NMIBC to effectively reduce disease recurrence or progression, whereas immune checkpoint inhibitors (ICIs) were recently introduced to treat advanced BLCA with good efficacy. For BCG and ICI applications, reliable biomarkers are required to stratify potential responders for better personalized interventions, and ideally, they can replace or reduce invasive examinations such as cystoscopy in monitoring treatment efficacy. Here we developed the cuproptosis-associated 11 gene signature (CuAGS-11) model to accurately predict survival and response to BCG and ICI regimens in BLCA patients. In both discovery and validation cohorts where BLCA patients were divided into high- and low-risk groups based on a median CuAGS-11 score as the cutoff, the high-risk group was associated with significantly shortened overall survival (OS) and progression-free survival (PFS) independently. The survival predictive accuracy was comparable between CuAGS-11 and stage, and their combination-based nomograms showed high consistence between predicted and observed OS/PFS. The analysis of 3 BLCA cohorts treated with BCG unveiled lower response rates and higher frequencies of recurrence or progression coupled with shorter survival in CuAGS-11 high-risk groups. In contrast, almost none of patients underwent progression in low-risk groups. In IMvigor210 cohort of 298 BLCA patients treated with ICI Atezolizumab, complete/partial remissions were 3-fold higher accompanied by significantly longer OS in the CuAGS-11 low- than high-risk groups (P = 7.018E-06). Very similar results were obtained from the validation cohort (P = 8.65E-05). Further analyses of Tumor Immune Dysfunction and Exclusion (TIDE) scores revealed that CuAGS-11 high-risk groups displayed robustly higher T cell exclusion scores in both discovery (P = 1.96E-05) and validation (P = 0.008) cohorts. Collectively, the CuAGS-11 score model is a useful predictor for OS/PFS and BCG/ICI efficacy in BLCA patients. For BCG-treated patients, reduced invasive examinations are suggested for monitoring the CuAGS-11 low-risk patients. The present findings thus provide a framework to improve BLCA patient stratification for personalized interventions and to reduce invasive monitoring inspections.
Introduction
Bladder cancers (BCs) or carcinomas (BLCAs) are the commonest urological malignancy worldwide, and up to 95% of them are originated from the urothelium (–). At diagnosis, the majority of BLCAs (70% – 80%) are non-muscle invasive (NMIBCs) while 20% - 30% present with muscle invasive BCs (MIBCs). NMIBCs display a high frequency of recurrence, but patients in general have a favorable outcome with long-term survival and only a small fraction (15% - 20%) advance into MIBCs (, ). Local and distant disseminations occur frequently in MIBCs, and many patients die from aggressive or metastatic diseases within 5 years (, ). During the last decades, the major clinical interventions of BLCAs largely include surgery plus intravesical Bacillus Calmette-Guerin (BCG) instillation for intermediate/high-risk NMIBCs and neoadjuvant chemotherapy for MIBCs (, –). BCG as a traditional immunotherapeutic strategy has been very successful in NMIBC treatments and this protocol is still recommended by the international guidelines as the standard care to reduce BC recurrence and progression in the present BLCA care (–, ). In the recent years, promoting anti-cancer immunity by using immune checkpoint inhibitors (ICIs) as a novel strategy has been developed for clinical application and totally revolutionized the therapeutic landscape of BLCAs and other cancer types (). By targeting PD-1/PD-L1, CTLA4, or other immune checkpoint proteins, the ICI therapy demonstrates robust efficacy in subsets of BLCA patients (, ). To improve stratification for better immunotherapeutic applications, numerous studies have paid great attention to biomarker identification for response to ICIs (, , –). Attempts to look for BCG treatment predictors are far behind, and the evaluation of BCG response depends mainly on cystoscopy, cytology and/or bladder biopsy nowadays (), although several molecules are shown to serve as potential factors (). In short, searching for reliable prognostic factors for patient survival, recurrence, NMIBC progression to MIBCs and treatment response are critical unmet needs, and patients with high-risk can thus be pinpointed for active surveillance and personalized intervention, thereby reducing BLCA-associated morbidity and mortality.
Copper, as an essential mineral nutrient, has long been appreciated to participate in cancer development and progression, and the copper signaling is actively involved in cancer cell proliferation, survival and metastasis (). More recently, Tsvetkov et al. defined a copper-dependent form of regulated cell death named cuproptosis (). During the cuproptotic process, FDX1, a reductase, and copper together bring on the lipoylation and aggregation of mitochondrial enzymes responsible for the tricarboxylic acid (TCA) cycle, and promote Fe-S cluster protein degradation, which consequently result in proteotoxic stress and eventual cell death (). It is currently unclear whether cuproptosis plays a role in carcinogenesis or copper-mediated oncogenic function can be targeted by inducing cuproptosis. Nevertheless, recent clinical investigations have shown that cuproptosis-related factors serve as predictors for outcomes and treatment response in several cancer types (–). By studying patients with clear cell renal cell carcinoma (ccRCC), we showed that the cuproptosis-associated 13 gene signature (CuAGS-13) was a robust predictor for outcomes and response to ICI and targeted therapies in ccRCC (). The association between cuproptosis and BLCA has also been explored using 10 cuproptosis factors, or cuproptosis-related genes and long non-coding RNAs (–). These different models consistently showed their prognostic values in outcome prediction of BLCA patients, and some of them also revealed the significant impact of cuproptosis-related factors on invasion, drug resistance, tumor microenvironments and immune cell infiltrations. However, it remains elusive whether the cuproptosis-based models are capable of predicting response to immunotherapy and survival benefits, the key clinical-related issues. Most of the above studies only examined the relationship between cuproptosis-based models and immune cell infiltrations in BLCA tumors (, , , ), while only in the report by Li et al, 34 BLCA patients treated with Atezolizumab (anti-PD-L1 antibody) were analyzed for their response rate using a 14 cuproptosis-related gene-containing signature (). Unfortunately, the authors ignored the IMvigor210 cohort with 348 BLCA patients receiving Atezolizumab. The results obtained from a small cohort of 34 patients, together with lack of survival analysis, are far from the conclusive proof of the cuproptosis-related effect on ICI therapy. In addition, BCG as an immunomodulator has long been successful in NMIBC treatment as described above (, ), but it is currently unclear whether cuproptosis or its associated gene signature can predict BCG response. In the present study, we analyzed several BLCA cohorts and developed the BLCA-specific cuproptosis-associated 11 gene signature (CuAGS-11) as a useful predictor for outcomes and therapeutic efficacy of BCG and ICIs.
Materials and methods
Study workflow, data collection and processing of BLCA tumors and bladder nontumorous tissues (NTs)
Based on 10 cuproptosis molecules () (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 (). Patient information, pathology/histology, transcriptome, mutation, and copy number variation (CNV) data were downloaded from https://gdc.cancer.gov/. Aneuploid score was from reference (). Tumor mutation burden (TMB) was calculated using Rpackage TCGAmutations. (ii) GSE13507 (, ), GSE154261 (), and GSE176307 () 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 (, ). 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 (). 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.
Figure 1
Identification of cuproptosis-associated genes whose expression correlates with 10 cuproptosis factors
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.
Figure 2
Construction of the CuAGS-11 risk score
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 receiver operating characteristic (ROC) curves and construction of a survival predictive nomogram
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.
Gene set enrichment analysis (GSEA)
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.
Analyses for proliferation, cancer stemness, and epithelial–mesenchymal transition (EMT) scores
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 (). Stemness score was assessed according to the mRNA expression-based stemness developed by Malta et al. (). EMT score was calculated based on the dbEMT signature from http://dbemt.bioinfo-minzhao.org/ ().
Tumor immune dysfunction and exclusion (TIDE) score analysis
TIDE score is evaluated according to myeloid-derived suppressor cell (MDSC), macrophage M2, T cell Dysfunction and Exclusion (). 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.
Statistical analysis
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.
Results
The CuAGS-11 model establishment based on the TCGA cohort of BLCA analysis
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) (). 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).
Enrichments of BLCA basal subtype and aggressive features in the CuAGS-11 high-risk tumors
BLCAs are in general stratified into luminal and basal subtypes according to their featured gene expression signatures (–). The luminal subtype is overrepresented by urothelium differentiation markers and transcription factors, while the basal one is poorly differentiated (, , ). 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).
Figure 3
Because the basal BLCA subtype is enriched with cycling and stem- and/or mesenchymal-like cells (, ), 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).
Genomic alterations and their association with the CuAGS-11 model
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).
The CuAGS-11 score for survival prediction in BLCAs
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).
Figure 4
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 (, ), 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).
Figure 5
The CuAGS-11 model as a predictor for response to BCG treatment
The clinic-pathological variables have been mainly used to evaluate response to BCG therapy (). We sought to determine whether the CuAGS-11 score could serve as such a predictive biomarker. The GSE154261 BLCA cohort () 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 () (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 () 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).
Figure 6
The CuAGS-11 model as a predictor for response to Atezolizumab therapy
ICI therapy has been applied for MIBCs with good efficacy in subsets of patients (). 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 (, ) (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 (). 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 () 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).
Figure 7
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 (). 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).
Discussion
BLCAs derived from urothelial cells are heterogenous and include NMIBCs and MIBCs with different outcomes and different clinical interventions (, ). Reliable biomarkers are highly required to stratify patient risk and to tailor treatment regimens (). Based on the cuproptosis-associated gene analysis, we develop the CuAGS-11 model for such a purpose. Our results presented herein demonstrate that CuAGS-11 serves as a useful predictor for BLCA patient survival and response to immunotherapies including BCG and ICIs.
Copper has long been appreciated to participate in oncogenesis (). Cuproptosis, a newly identified form of regulated cell death (RCD), is copper-dependent cell death caused by FDX1-mediated mitochondrial protein lipoylation (). It is currently unclear whether cuproptosis, like apoptosis or other types of RCD, plays any parts in carcinogenesis. We recently developed a CuAGS-13 model to accurately predict ccRCC outcomes and patient response to targeted and ICI therapies (), however, this same model failed to do so in BLCAs (data not shown). Intriguingly, the direct application of 10 cuproptosis factors did not show any prognostic values in either BLCAs or ccRCCs. Moreover, the genes in the CuAGS-11 model are totally different from those in the CuAGS-13 () or from any other cuproptosis-related models. Thus, cuproptosis-based models are context-dependent and their prognostic powers very significantly. In the present study, the CuAGS-11 model is constructed based on 10 cuproptosis factors, but the scores are strongly correlated with many key aggressive characteristics in BLCAs in positive manners. Thus, the CuAGS-11 model represents a classifier integrated with pathological/molecular and many other BLCA features. Despite such broad connections, the CuAGS-11 score predicts BLCA patient survival independently of any other variables.
BCG, a live attenuated strain of Mycobacterium bovis, was first applied for patients with NMIBC in 1976 (), and since then, intravesical BCG therapy has been the most effective treatment strategy to prevent recurrence for intermediate- and high-risk NMIBCs (). If without adjuvant therapies, recurrence occurs in up to 70% of patients who received tumor resection (). Mechanisms underlying BCG against BLCAs remain incompletely understood, but accumulated studies indicate that BCG as an immunomodulator stimulates both innate and acquired immune responses, thereby exerting a therapeutic efficacy (, ). Li et al. recently showed that BCG treatment failure was associated with enhanced PD-L1 and FGFBP1 expression (). In addition, BCG may have direct impacts on BLCA cells. Despite tremendous advances in next generation sequencing and other high-throughput technologies, there is still lack of reliable molecular predictors for response to BCG. Clinic-pathological variables combined with cystoscopy have so far served as major approaches to predict and measure potential response to BCG (, ). In clinical practice nowadays, cystoscopy, cytology and/or bladder biopsy are used to determine response at 3 months and 6 months following a BCG induction regimen (, ). Cystoscope examination is invasive and costly. In the present investigations, we showed a high accuracy of the CuAGS-11 model in prognosticating BCG responders. Based on our analyses of 3 BLCA cohorts treated with BCG, almost all patients with disease progression were observed in the CuAGS-11 high-risk group. For T1 patients, no recurrence occurred in the CuAGS-11 low-risk group. These proof-of-concept results suggest that the invasive examination may not be required or at least reduced for CuAGS-11 low-risk patients, which is worthy of further studies.
In the recent years, the ICI therapy has been applied for advanced cancers including MIBCs (). As only subsets of MIBCs respond to ICIs, distinguishing potential responders from non-responders should greatly contribute to personalized application of ICIs and the development of accurate predictive biomarkers is thus critical. Here we found that the CuAGS-11 model was similarly useful for stratifying ICI responders in 2 cohort patients treated with Atezolizumab. In the IMvigor210 cohort, the CR/PR rate was > 3-fold higher in the low- than high-risk group. The GSE176307 cohort analysis showed that all 17 patients underwent progression in the CuAGS-11 high-risk group, whereas 9 of 17 (54.7%) low-risk patients obtained CR or PR. Consistently, patients in the low-risk group had significantly longer PFS and/or OS. TMB has been shown as a reliable predictor for response to ICI therapy in BLCAs (), but the CuAGS-11 signature was not related to TMB. TIDE analyses of these two cohorts consistently showed that CuAGS-11 high-risk tumors were characterized by T cell exclusion, which suggests that the CuAGS-11 model can help identify BC tumors with an immunologically cold phenotype. Poor response to ICIs is thus conceivable in patients with CuAGS-11 high tumors. It is currently unclear what is a mechanistic link between CuAGS-11 score and T cell exclusion in BLCAs and whether cuproptosis is involved in immune cell fate decision or T cell-mediated tumor destruction. Elucidating these important issues should contribute to improvement of BLCA immunotherapy.
In conclusion, we constructed the CuAGS-11 score model for prediction of survival and response to BCG and ICI therapies in BLCAs. This model, although derived from cuproptosis-associated factors, is a classifier integrated with molecular and other features of BLCA. Importantly, for patients receiving BCG, recurrence occurs predominantly in the CuAGS-11 high-risk group, and no disease progression was observed in the low-risk patients. Thus, it may be unnecessary to monitor the CuAGS-11 low-risk patients using invasive examinations routinely. The present findings further show that the CuAGS-11 model is helpful to identify BLCA tumors with an immunologically cold phenotype and to distinguish between ICI responders and non-responders. Taken together, the CuAGS-11 score model may significantly improve BLCA patient stratification for tailored patient interventions, reducing BLCA-associated morbidity and mortality. It is worth validating these observations in clinical practices.
Statements
Data availability statement
Source data downloaded from public databases are provided with this paper. Any additional information required to reanalyze the data reported in this paper is available from the corresponding authors upon reasonable request.
Author contributions
HY, TL, YF, and DX designed the study. HY and YX performed bioinformatic analysis. YF and DX supervised the study. HY, TL, YF, and DX wrote and revised the manuscript. All authors reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by grants from Scientific Research Foundation of Qilu Hospital of Shandong University (Qingdao) (No. QDKY2019QN17), National Natural Science Foundation of China (No. 82103557, and 81972475), Shandong Provincial Natural Science Foundation (No. ZR2020QH245), the Swedish Cancer Society (No. 19 0018 Pj), Swedish Research Council (2018–02993), the Cancer Society in Stockholm (201393), and Karolinska Institutet (2022–01889).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1126247/full#supplementary-material
References
1
CambierSSylvesterRJColletteLGonteroPBrausiMAvan AndelGet al. EORTC nomograms and risk groups for predicting recurrence, progression, and disease-specific and overall survival in non-muscle-invasive stage Ta-T1 urothelial bladder cancer patients treated with 1-3 years of maintenance bacillus calmette-guerin. Eur Urol (2016) 69(1):60–9. doi: 10.1016/j.eururo.2015.06.045
2
SiegelRLMillerKDJemalA. Cancer statistics, 2020. CA Cancer J Clin (2020) 70(1):7–30. doi: 10.3322/caac.21590
3
PatelVGOhWKGalskyMD. Treatment of muscle-invasive and advanced bladder cancer in 2020. CA Cancer J Clin (2020) 70(5):404–23. doi: 10.3322/caac.21631
4
ChoiWCzerniakBOchoaASuXSiefker-RadtkeADinneyCet al. Intrinsic basal and luminal subtypes of muscle-invasive bladder cancer. Nat Rev Urology (2014) 11(7):400–10. doi: 10.1038/nrurol.2014.129
5
KnowlesMAHurstCD. Molecular biology of bladder cancer: new insights into pathogenesis and clinical diversity. Nat Rev Cancer (2015) 15(1):25–41. doi: 10.1038/nrc3817
6
JiangSRedelman-SidiG. BCG In bladder cancer immunotherapy. Cancers (Basel) (2022) 14(13):3073. doi: 10.3390/cancers14133073
7
ReikeMJIngersollMAMullerDCZuiverloonTCMStrandgaardTKamatAMet al. Biology of BCG response in non-muscle invasive bladder cancer - 2021 IBCN updates part III. Urol Oncol (2022), S1078-1439(22)00340-4. doi: 10.1016/j.urolonc.2022.09.012
8
KamatAMFlaigTWGrossmanHBKonetyBLammDO’DonnellMAet al. Expert consensus document: consensus statement on best practice management regarding the use of intravesical immunotherapy with BCG for bladder cancer. Nat Rev Urology (2015) 12(4):225–35. doi: 10.1038/nrurol.2015.58
9
ChenHYangWXueXLiYJinZJiZ. Neoadjuvant immunotherapy and chemoimmunotherapy for stage II-III muscle invasive bladder cancer. Front Immunol (2022) 13:986359. doi: 10.3389/fimmu.2022.986359
10
MoralesAEidingerDBruceAW. Intracavitary bacillus calmette-guerin in the treatment of superficial bladder tumors. J Urol (1976) 116(2):180–3. doi: 10.1016/S0022-5347(17)58737-6
11
SamsteinRMLeeCHShoushtariANHellmannMDShenRJanjigianYYet al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet (2019) 51(2):202–6. doi: 10.1038/s41588-018-0312-8
12
ManciniMRighettoMNoessnerE. Checkpoint inhibition in bladder cancer: clinical expectations, current evidence, and proposal of future strategies based on a tumor-specific immunobiological approach. Cancers (Basel) (2021) 13(23):6016. doi: 10.3390/cancers13236016
13
YiXZhengXXuHLiJZhangTGePet al. IGFBP7 and the tumor immune landscape: a novel target for immunotherapy in bladder cancer. Front Immunol (2022) 13:898493. doi: 10.3389/fimmu.2022.898493
14
de KouchkovskyIZhangLPhilipEJWrightFKimDMNatesanDet al. TERT promoter mutations and other prognostic factors in patients with advanced urothelial carcinoma treated with an immune checkpoint inhibitor. J Immunother Cancer (2021) 9(5):e002127. doi: 10.1136/jitc-2020-002127corr1
15
LiuTLiSXiaCXuD. TERT promoter mutations and methylation for telomerase activation in urothelial carcinomas: new mechanistic insights and clinical significance. Front Immunol (2022) 13:1071390. doi: 10.3389/fimmu.2022.1071390
16
WangZSoWZLohKYLimYKMahendranRWuQHet al. Predictive biomarkers of response to bacillus calmette-guerin immunotherapy and bacillus calmette-guerin failure for non-muscle invasive bladder cancer. Int J Urol (2022) 29(8):807–15. doi: 10.1111/iju.14921
17
GeEJBushAICasiniACobinePACrossJRDeNicolaGMet al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat Rev Cancer (2022) 22(2):102–13. doi: 10.1038/s41568-021-00417-2
18
TsvetkovPCoySPetrovaBDreishpoonMVermaAAbdusamadMet al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science (2022) 375(6586):1254–61. doi: 10.1126/science.abf0529
19
YuanHQinXWangJYangQFanYXuD. The cuproptosis-associated 13 gene signature as a robust predictor for outcome and response to immune- and targeted-therapies in clear cell renal cell carcinoma. Front Immunol (2022) 13:971142. doi: 10.3389/fimmu.2022.971142
20
ZhuBWangSWangRWangX. Identification of molecular subtypes and a six-gene risk model related to cuproptosis for triple negative breast cancer. Front Genet (2022) 13:1022236. doi: 10.3389/fgene.2022.1022236
21
ZhangZZengXWuYLiuYZhangXSongZ. Cuproptosis-related risk score predicts prognosis and characterizes the tumor microenvironment in hepatocellular carcinoma. Front Immunol (2022) 13:925618. doi: 10.3389/fimmu.2022.925618
22
YangWWuHTongLWangYGuoQXuLet al. A cuproptosis-related genes signature associated with prognosis and immune cell infiltration in osteosarcoma. Front Oncol (2022) 12:1015094. doi: 10.3389/fonc.2022.1015094
23
WuMBaoJLeiYTaoSLinQChenLet al. Comprehensive analysis of the cuproptosis-related model to predict prognosis and indicate tumor immune infiltration in lung adenocarcinoma. Front Oncol (2022) 12:935672. doi: 10.3389/fonc.2022.935672
24
WangGXiaoRZhaoSSunLGuoJLiWet al. Cuproptosis regulator-mediated patterns associated with immune in fi ltration features and construction of cuproptosis-related signatures to guide immunotherapy. Front Immunol (2022) 13:945516. doi: 10.3389/fimmu.2022.945516
25
SongSZhangMXiePWangSWangY. Comprehensive analysis of cuproptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front Immunol (2022) 13:978909. doi: 10.3389/fimmu.2022.978909
26
ZhangSYuSDuanHXiaWWangCShenH. Revealing prognostic and tumor microenvironment characteristics of cuproptosis in bladder cancer by genomic analysis. Front Genet (2022) 13:997573. doi: 10.3389/fgene.2022.997573
27
SongQZhouRShuFFuW. Cuproptosis scoring system to predict the clinical outcome and immune response in bladder cancer. Front Immunol (2022) 13:958368. doi: 10.3389/fimmu.2022.958368
28
LiHZuXHuJXiaoZCaiZGaoNet al. Cuproptosis depicts tumor microenvironment phenotypes and predicts precision immunotherapy and prognosis in bladder carcinoma. Front Immunol (2022) 13:964393. doi: 10.3389/fimmu.2022.964393
29
ZhangYLiXLiXZhaoYZhouTJiangXet al. Comprehensive analysis of cuproptosis-related long noncoding RNA immune infiltration and prediction of prognosis in patients with bladder cancer. Front Genet (2022) 13:990326. doi: 10.3389/fgene.2022.990326
30
BaiYZhangQLiuFQuanJ. A novel cuproptosis-related lncRNA signature predicts the prognosis and immune landscape in bladder cancer. Front Immunol (2022) 13:1027449. doi: 10.3389/fimmu.2022.1027449
31
RobertsonAGKimJAl-AhmadieHBellmuntJGuoGCherniackADet al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell (2017) 171(3):540–56.e25. doi: 10.1016/j.cell.2017.09.007
32
ThorssonVGibbsDLBrownSDWolfDBortoneDSOu YangTHet al. The immune landscape of cancer. Immunity (2018) 48(4):812–30 e14. doi: 10.1016/j.immuni.2018.03.023
33
KimWJKimEJKimSKKimYJHaYSJeongPet al. Predictive value of progression-related gene classifier in primary non-muscle invasive bladder cancer. Mol Cancer (2010) 9:3. doi: 10.1186/1476-4598-9-3
34
LeeJSLeemSHLeeSYKimSCParkESKimSBet al. Expression signature of E2F1 and its associated genes predict superficial to invasive progression of bladder tumors. J Clin Oncol (2010) 28(16):2660–7. doi: 10.1200/JCO.2009.25.0977
35
RobertsonAGGroeneveldCSJordanBLinXMcLaughlinKADasAet al. Identification of differential tumor subtypes of T1 bladder cancer. Eur Urol (2020) 78(4):533–7. doi: 10.1016/j.eururo.2020.06.048
36
RoseTLWeirWHMayhewGMShibataYEulittPUronisJMet al. Fibroblast growth factor receptor 3 alterations and response to immune checkpoint inhibition in metastatic urothelial cancer: a real world experience. Br J Cancer (2021) 125(9):1251–60. doi: 10.1038/s41416-021-01488-6
37
MariathasanSTurleySJNicklesDCastiglioniAYuenKWangYet al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501
38
BalarAVGalskyMDRosenbergJEPowlesTPetrylakDPBellmuntJet al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet (2017) 389(10064):67–76. doi: 10.1016/S0140-6736(16)32455-2
39
RitchieMEPhipsonBWuDHuYLawCWShiWet al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007
40
MotzerRJBanchereauRHamidiHPowlesTMcDermottDAtkinsMBet al. Molecular subsets in renal cancer determine outcome to checkpoint and angiogenesis blockade. Cancer Cell (2020) 38(6):803–17.e4. doi: 10.1016/j.ccell.2020.10.011
41
MaltaTMSokolovAGentlesAJBurzykowskiTPoissonLWeinsteinJNet al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell (2018) 173(2):338–54.e15. doi: 10.1016/j.cell.2018.03.034
42
ZhaoMLiuYZhengCQuH. dbEMT 2.0: an updated database for epithelial-mesenchymal transition genes with experimentally verified information and precalculated regulation information for cancer metastasis. J Genet Genomics (2019) 46(12):595–7. doi: 10.1016/j.jgg.2019.11.010
43
JiangPGuSPanDFuJSahuAHuXet al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
44
MoQNikolosFChenFTramelZLeeYCHayashiKet al. Prognostic power of a tumor differentiation gene signature for bladder urothelial carcinomas. J Natl Cancer Inst (2018) 110(5):448–59. doi: 10.1093/jnci/djx243
45
SjodahlGErikssonPLiedbergFHoglundM. Molecular classification of urothelial carcinoma: global mRNA classification versus tumour-cell phenotype classification. J pathology (2017) 242(1):113–25. doi: 10.1002/path.4886
46
DamrauerJSHoadleyKAChismDDFanCTiganelliCJWobkerSEet al. Intrinsic subtypes of high-grade bladder cancer reflect the hallmarks of breast cancer biology. Proc Natl Acad Sci U S A (2014) 111(8):3110–5. doi: 10.1073/pnas.1318376111
47
ChoiWPortenSKimSWillisDPlimackERHoffman-CensitsJet al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell (2014) 25(2):152–65. doi: 10.1016/j.ccr.2014.01.009
48
GuoYYuanXLiKDaiMZhangLWuYet al. GABPA is a master regulator of luminal identity and restrains aggressive diseases in bladder cancer. Cell Death differentiation (2020) 27(6):1862–77. doi: 10.1038/s41418-019-0466-7
49
HedegaardJLamyPNordentoftIAlgabaFHoyerSUlhoiBPet al. Comprehensive transcriptional analysis of early-stage urothelial carcinoma. Cancer Cell (2016) 30(1):27–42. doi: 10.1016/j.ccell.2016.05.004
50
CastellanoESambaCEstesoGSimpsonLVendrameEGarcia-CuestaEMet al. CyTOF analysis identifies unusual immune cells in urine of BCG-treated bladder cancer patients. Front Immunol (2022) 13:970931. doi: 10.3389/fimmu.2022.970931
51
LiFZhangHWangYYaoZXieKMoQet al. FGFBP1 as a potential biomarker predicting bacillus calmette-guerin response in bladder cancer. Front Immunol (2022) 13:954836. doi: 10.3389/fimmu.2022.954836
Summary
Keywords
BCG, cuproptosis, ICI, immunotherapy, prognosis, bladder carcinoma
Citation
Yuan H, Xiu Y, Liu T, Fan Y and Xu D (2023) The cuproptosis-associated 11 gene signature as a predictor for outcomes and response to Bacillus Calmette-Guerin and immune checkpoint inhibitor therapies in bladder carcinoma. Front. Immunol. 14:1126247. doi: 10.3389/fimmu.2023.1126247
Received
17 December 2022
Accepted
21 April 2023
Published
03 May 2023
Volume
14 - 2023
Edited by
Jianzhong Ai, Sichuan University, China
Reviewed by
Neveen Said, Wake Forest Baptist Medical Center, United States; Jamshid Hadjati, Tehran University of Medical Sciences, Iran
Updates
Copyright
© 2023 Yuan, Xiu, Liu, Fan and Xu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yidong Fan, fanyd@sdu.edu.cn; Tiantian Liu, liu.tiantian@sdu.edu.cn; Dawei Xu, Dawei.Xu@ki.se
†ORCID: Dawei Xu, orcid.org/0000-0003-3141-4524
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.