Three independent cohorts consisting of MIBC patients treated with ICB therapy and the Cancer Genome Atlas (TCGA) cohort were included for analysis. The detailed patient characteristics can be seen in Table 1. Normalization and log₂-transformation were conducted in all RNA-seq and microarray data.

The training cohort enrolled 168 MIBC patients three times weekly treated with 1,200 mg atezolizumab from IMvigor210 trial, patients with other cancer types or unconfirmed overall responses were excluded, corresponding data were downloaded by the IMvigor210CoreBiologies R package (18). A total of 76 MIBC patients treated with at least one dose of anti-PD-1 or anti-PD-L1 immunotherapy with response and survival information were recruited as validation I cohort (19). A small series of expression profiles taken from 11 MIBC patients treated with pembrolizumab was considered as validation II cohort used for expression validation. The last validation III cohort contained 386 MIBC patients from TCGA database. Full transcriptome data and characteristics of the patients were assessed from http://www.cbioportal.org/ in July 2022. Patients achieving complete response (CR) or partial response (PR) were regarded as responders, while patients with progressive disease (PD) or stable disease (SD) were defined as non-responders in the training and validation I cohort. Responses of immunotherapy in TCGA cohort were conducted using the Tumor immune dysfunction and exclusion (TIDE) analysis, responders were patients with TIDE score <0, otherwise, non-responders (20).

As illustrated in Figure 1, our study included three phases. In the discovery phase, we measured the performances of 50 cancer hallmarks, identified and validated IFNα response as the leading factor for the favorable prognosis of ICB therapy from the survival, functional and tumor microenvironment perspectives. Secondly, differentially expressed prognostic genes mapping on IFNα response pathways were included for machine learning algorithms to select robust genes with better performance in therapy response prediction. Moreover, the expression levels and predictive abilities of the above candidate genes were verified in external validation cohorts.

In the training cohort, we quantified levels of cancer hallmarks based on transcriptional profiles and hallmark gene sets acquired from the Molecular Signatures Database (MsigDB) through ssGSEA algorithms (21). Subsequently, the prognostic significance of cancer hallmarks was evaluated by univariate Cox analysis in MIBC patients after ICB therapy. We also applied gene set enrichment analysis (GSEA) to compare the enriched pathways between responders and non-responders referring to the hallmark gene sets (22). As a result, IFN-α response (IFNAR) was found with the lowest hazard ratio (HR) value. Subsequently, according to the IFNAR-related score, we divided patients into high-, middle- and low-score groups, the survival differences, responder proportion, PD-L1 protein expression levels on immune cells (IC) and immune phenotype were compared among groups.

The overall infiltration of immune cells, stromal cells and tumor cell purity were inferred by ESTIMATE algorithm (23). Moreover, we evaluated the infiltration of 22 immune cell subpopulations in MIBC biopsies through CIBERSORT (Cell-type Identification by Estimating Relative Subsets of RNA Transcripts), a deconvolution algorithm to characterize immune cell composition using gene expression profiles (24). Immune cell abundance variations in the high- and low- IFNAR-related score groups were detected, correlations between immune cell infiltration and IFNAR were also calculated.

Ninety-five genes that participated in the process of IFNAR were derived from MSigDB. We then screened differentially expressed genes (DEGs) between responders and non-responders through “limma” package when false discovery rate (FDR) was less than 0.05 (25). Meanwhile, prognostic genes were identified by univariate Cox analysis with a threshold of p < 0.05. Prognostic DEGs were selected for further analysis. Subsequently, two machine learning algorithms, least absolute shrinkage and selection operator (LASSO) Cox regression and random survival forest (RSF) analysis were commonly applied to perform gene selection. LASSO used 10-fold cross-validation to estimate the penalty parameters by “glmnet” package to avoid over-fitting. RSF adapts random forests to survival analysis based on ensemble trees. Variable importance (VIMP) evaluates the predictive ability alterations of RSF model when genes are randomly permuted, higher VIMP indicates greater significance, while the average depth of genes among all survival trees was implied by minimal depth, smaller values suggest increased importance (26). They measure the impact of genes from different points of view, eligible genes identically selected by LASSO, VIMP and minimal depth were obtained for further validation.

To predict the therapeutic response to ICB therapy of MIBC patients in TCGA cohort, we applied TIDE algorithm to evaluate diverse mechanisms in tumor immune evasion, comprising of immunosuppressive factors induced cytotoxic T lymphocytes (CTLs) exclusion and dysfunction. Before analysis in TCGA cohort, we conducted TIDE in the training and validation I cohort to test its predictive ability in ICB therapy response. The infiltration of immunosuppressive myeloid suppressor cells (MDSC), cancer-associated fibroblasts (CAFs) and M2 subtypes of tumor-associated macrophages (TAM.M2) were also evaluated (20). Lower TIDE scores imply better clinical efficacy of immune checkpoint inhibitors.

We performed all statistical analyses with R software. The D’Agostino and Pearson omnibus normality tests were initially carried out to assess whether the data fit a normal distribution. When parameters were normally distributed, a two-tailed unpaired t-test, one-way ANOVA with Tukey’s correction and the Pearson correlation would be conducted. Once data did not achieve the assumptions of parametric tests, the Mann–Whitney U test, one-way ANOVA using Kruskal–Wallis with Dunn’s correction and Spearman correlation would be employed. Results met the level of 5% (p < 0.05) were considered statistically significant.

In the training cohort, we quantified the performance of 50 cancer hallmarks, then each HR value was calculated and ranked through univariate Cox analysis (Supplementary Table S1; Supplementary Figure S1). Among hallmarks, IFNAR ranked first with the lowest HR value (Figure 2A, HR = 0.667, p = 0.008), and high IFNAR-related score indicated favorable overall survival (OS) (Figure 2B, p = 0.002). Moreover, GSEA also revealed that IFNAR was significantly higher enriched in patients with responses to ICB therapy (Figure 2C; Supplementary Table S2). Based on the IFNAR-related score, we equally divided patients into three groups, the fraction of patients achieving CR or PR were 37.5%, 21.4%, and 16.1% in the high-, middle- and low-score group, respectively (Figure 2D, p = 0.024). Besides, the IFNAR-related score was higher in the group of patients achieving CR, patients with the higher protein expression level of PD-L1 on IC and patients with inflamed immune phenotype than other groups (Figures 2E–G). Likewise, in the validation I cohort, IFNAR was ordered first among cancer hallmarks (Figure 2H, HR = 0.640, p = 0.006) (Supplementary Table S3), high IFNAR-related score was linked to better OS (Figure 2I, p = 0.041), IFNAR was significantly upregulated in responders (Figure 2J; Supplementary Table S4). Together, our study suggested that the IFNAR pathway was a leading favorable factor with promising predictive value for the prognosis of ICB therapy in MIBC patients.

It has been well established that the tumor microenvironment influences clinical efficacy of immunotherapy. Therefore, we utilized ESTIMATE and CIBEROSRT to calculate the immune cell infiltration in tumor tissues. Both for training and validation I cohorts, the high-IFNAR-related-score subgroup exhibited higher immune and stromal scores, lower tumor purity than the low-score subgroup did (Figure 3A). In particular, immune effector cells, including CD8⁺ T cells, CD4⁺ memory-activated T cells and type 1 proinflammatory macrophage (M1) were increasingly infiltrated, while CD4⁺ memory-resting T cells were decreasingly infiltrated in the high-IFNAR-related-score subgroup in both cohorts (Figure 3B). Furthermore, the IFNAR-related score was positively correlated with above mentioned immune effector cells and negatively associated with CD4⁺ memory-resting T cells (Figure 3C). Collectively, these results illustrated high levels of IFNAR were accompanied by an immune-active tumor microenvironment in MIBC.

We acquired IFNAR-related genes (n = 95) for subsequent analyses. With the threshold of FDR <0.05, twenty-four upregulated DEGs and 5 downregulated DEGs in tumor tissues of responders compared to non-responders were identified (Figure 4A). Meanwhile, seventeen prognostic genes with p < 0.05 were also discovered (Figure 4B). Ten intersected prognostic DEGs were then applied to machine learning algorithms (Figure 4C). Then four genes (CXCL10, LAMP3, TAP1, TRIM5) were selected by LASSO Cox regression analysis (Figure 4D), three genes (CXCL10, LAMP3, IRF1) were acquired based on VIMP and minimal depth through RSF (Figure 4E). Among them, CXCL10 and LAMP3 were co-selected as promising predictive biomarkers for the responses of ICB therapy.

In the training cohort, CXCL10 and LAMP3 were significantly higher expressed in responders (Figure 5A, p = 0.006, p = 0.019, respectively). The same trend was also observed in another two validation cohorts (Figures 5B, C). Besides, both CXCL10 and LAMP3 expression in the tumor tissues indicated favorable OS in the training cohort (Figure 5D, p = 0.006, p = 0.001, respectively) and validation I cohort (Figure 5E, p = 0.029, p = 0.005, respectively). Moreover, we noted the patients who achieved CR or PR were more frequent in CXCL10 ^high LAMP3 ^high subgroups, and the combination of CXCL10 and LAMP3 expression significantly indicated improved OS in both cohorts (Figures 5F, G). Moreover, the expression levels of CXCL10 and LAMP3 were higher in patients with inflamed immune phenotype and higher PD-L1 protein expression levels on IC than in other groups (Supplementary Figure S2). Together, our study suggested that CXCL10 plus LAMP3 could serve as an ideal and stable predictive biomarker for patient response to ICB therapy in clinical settings.

In the training and validation I cohort, lower TIDE scores were seen in patients with CR/PR, and patients in the high-TIDE-score group showed worse overall survival (Supplementary Figure S3), which proved the predictive ability of TIDE in ICB therapy response. To further evaluate the universal applicability of CXCL10 and LAMP3 in predicting the responsiveness to ICB therapy, we evaluated therapy responses of MIBC patients through TIDE algorithm in TCGA cohort. The TIDE score was significantly higher in the low-IFNAR-related-score subgroup, which indicates lower therapeutic sensitivity (Figure 6A, p < 0.001). Furthermore, a strong correlation was found between IFNAR-related scores and TIDE scores (R = −0.32, p = 1.7e-10), the IFNAR-related score was also negatively associated with exclusion (R = −0.14, p = 0.006) and dysfunction (R = −0.41, p = 2.2e–16) (Figure 6B). As for the combination of CXCL10 and LAMP3, the TIDE score was lower and the responder frequency was higher in the CXCL10 ^high LAMP3 ^high subgroup (Figures 6C, D). Intriguingly, the scores of three reported immunosuppressive cells that suppress tumor T-cell infiltration, consisting of CAFs, MDSCs and TAM.M2 were higher in the CXCL10 ^low LAMP3 ^low subgroup (Figure 6E). The above results implied that high expression levels of CXCL10 and LAMP3 may predict a tumor microenvironment that favors immunotherapeutic response and indicate better ICB therapy response.