Gene expression data of bladder cancer with clinical variables were obtained from The Cancer Genome Atlas (TCGA) database, and patients with MIBC (T2 to T4) were selected for subsequent analysis. Another independent dataset was downloaded from Gene Expression Omnibus (GEO) for external validation. The ESTIMATE algorithm was used to calculate scores to predict the level of immune and stromal cell infiltration for each patient (Yoshihara et al., 2013). The R package “ESTIMATE” has been widely utilized in cancer-related studies (Liu et al., 2021a; Liu et al., 2021b; Liu et al., 2021c). The optimal cutoff values of immune and stromal scores were determined with maximally selected log-rank using the R package “survminer” (Li et al., 2020), and patients were then divided into immune-high/-low and stromal-high/-low groups. Kaplan–Meier curves were utilized to evaluate the association of immune and stromal scores with survival outcomes.

Differential expression genes (DEGs) between immune-high versus immune-low and stromal-high versus stromal-low subgroups were identified using the “limma” R package with a setting of |fold change| > 2 and a p value < 0.05, visualized with a heatmap and intersected with a Venn plot. Enrichment analyses of gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were performed to reveal the biological process, cellular component, molecular function, and molecular pathways that the intersected DEGs were associated with.

The protein–protein interaction (PPI) network of the DEGs was constructed with the STRING database (https://string-db.org/), with an interaction confidence of 0.99, and the core modules of the network were identified, defined as a collection of genes with no less than three nodes within the network. The prognostic value of the genes present in the core modules of the PPI network was evaluated with Kaplan–Meier curves. Then an independent GEO dataset GSE31684 was used for the external validation of the prognostic genes. A p value < 0.05 indicates that the correlation is significant.

The first 100 most correlative genes of each prognostic gene were identified and then intersected. Subsequent enrichment analyses and network construction were performed. The expression of the prognostic genes was compared between tumor and normal samples with data downloaded from the UALCAN cancer database (Chandrashekar et al., 2017) and the Human Protein Atlas (Uhlén et al., 2015). The connection between the prognostic genes and clinical variables was presented in a Sanguini diagram using the R package “ggalluval.” The percentage abundance of TICs was predicted and displayed using the R package “pheatmap,” and the CIBERSORT algorithm was utilized to compare the distribution of TICs according to the expression level of prognostic genes with the R package “ggplot.” We further evaluated the expression of eight immune checkpoint–relevant genes to reveal a potential role of the prognostic genes in immunotherapy. Correlation between TICs, immune checkpoint–relevant genes expression levels, and TMB with LCK and CD3E was further performed for pan-cancers, including 32 kinds of tumors. Another two independent bladder cancer cohorts receiving anti-PD1/PD-L1 inhibitor immunotherapy (GSE176307 (Rose et al., 2021) and IMvigor210 (Mariathasan et al., 2018)) were recruited to validate the predictive value of LCK and CD3E for immune response (responder: partial response [PR] or complete response [CR]; non-responder: stable disease [SD]; or progressed disease [PD]). The R package “Maxstat,” “Survminer,” “Survival,” and “ggplot2” were used to assess the prognostic significance of LCK and CD3E. The correlation of LCK and CD3E with immune phenotypes (inflamed, excluded, and desert) in the IMvigor210 cohort was also analyzed.

A total of 361 and 78 MIBC patients were identified from TCGA and GSE31684 (Table 1), respectively. The optimal cutoff values were 1157.37 for the immune score and -1106.15 for the stromal score (Figures 1A, B). In total, 92 patients from TCGA were assigned to the immune-high group and 269 to the immune-low group, while 253 patients were classified as stromal-high and 108 stromal-low. Higher immune infiltration was associated with improved overall survival (p = 0.042), while increased stromal infiltration was associated with worse overall survival (p = 0.009) (Figures 1C, D).

In total, 2033 DEGs were identified in stromal-high/low groups and 1843 DEGs in the immune-high/low groups (Figures 1E, F). Then 1234 DEGs were eventually intersected (Figure 2A), and the enrichment analyses indicated a role in cytokine–cytokine receptor interaction and T-cell activation (Figure 2B). Network construction identified a core module of 173 genes (Figure 2C). Survival curves demonstrated that CD3E (TCGA: p = 0.041, GEO: p = 0.022) and LCK (TCGA: p = 0.026, GEO: p = 0.024) were the prognostic genes of MIBC after external validation with GSE31684 (Figures 2D–G).

The top 100 genes that co-express with LCK and CD3E are partially presented in Figures 3A, B. Interestingly, CD3E was the most correlative gene of LCK, with a Spearman coefficient r = 0.86 (p < 0.001) (Figures 3C, D). The intersection of the top 100 co-expressed genes of LCK and CD3E were found to be associated with T-cell activation and differentiation and the T-cell receptor signaling pathway (Figures 3E–G).

The expression of LCK and CD3E was found to be much lower in tumor samples than in normal samples, with a median expression value (transcript per million) of 6.685 vs 3.116 for LCK (Figure 4A) and 17.484 vs 5.58 for CD3E (Figure 4B). Immunohistochemical staining showed a consistent expression trend of LCK (Figures 4C, D) and CD3E (Figures 4E, F) between normal tissue and bladder cancer samples. The association between clinical variables, including age, gender, and pathological stage, with the expression of LCK and CD3E, was also displayed, from which we could observe that there was a tendency of the distribution of high/low LCK/CD3E across different pathologic stages of MIBC (Figures 4G, H). Supplementary Figure S1 more quantitatively shows that MIBC patients in earlier stages (stages II and III vs stage IV) had a slightly higher percentage of high-LCK and high-CD3E expression (high-LCK: 53.11 vs. 44.19%; high-CD3E: 51.55vs.vs 48.06%; Supplementary Figure S1A). Further subgroup analyses consistently indicated that higher expressions of LCK and CD3E were found in MIBC patients diagnosed with earlier T stage (T2&T3 vs. T4: high-LCK: 50.16vs.vs 48.28%; high-CD3E: 50.32 vs 46.55%; Supplementary Figure S1B), and in particular earlier N stage (N0 vs. N1 and N2 vs. N3: high-LCK: 52.04 vs. 46.22vs.vs 28.57%; high-CD3E: 51.58 vs 49.58% vs. 28.57%; Supplementary Figure S1C) and earlier M stage (M0 vs. M1: high-LCK: 45.76 vs. 0%; high-CD3E: 42.94 vs. 12.50%; Supplementary Figure S1D).

The percentage abundance of TICs in MIBC patients is shown in Supplementary Figure S2 and Supplementary Table S1, from which we have a general view of the percentage of different TICs in MIBC samples. The CIBERSORT algorithm revealed the significantly different distribution of TICs between LCK-high/low and CD3E-high/low subgroups; for instance, the infiltration of memory B cells, CD8⁺ T cells, CD4⁺memory activated, and macrophage M1 was higher in the LCK-high and CD3E-high subgroups (Figures 5A, B). Further analysis indicated a higher expression of several key immune checkpoint–relevant genes such as PD-1, PD-L1, PD-L2, and CTLA4 in LCK-high and CD3E-high subgroups (Figures 5C, D). Moreover, this positive correlation between memory B cell, CD8⁺ T cell, CD4⁺ memory activated, and macrophage M1 with LCK and CD3E was also consistent across different cancer types (Figures 6A,B).

Interestingly, LCK and CD3E had a generally positive correlation with those eight immune checkpoint genes for pan-cancers, except for thymoma (Figures 7A, B). Pan-cancer analysis also revealed a positive correlation between LCK/CD3E and TMB in several types of cancers, despite the coefficients being small (Figures 7C, D). Further analyses showed that the correlation of LCK and CD3E with TMB in bladder cancer was weak, although a statistical significance was reached (Figures 7E, F). Last but not least, there was a moderate correlation between CD39, a previously reported gene predictive of response rate to PD-1 inhibitors, with LCK (r = 0.48, p < 0.001) and CD3E (r = 0.52, p < 0.001) (Figures 7G, H).

Figures 8A, B show a greater percentage of high LCK (62.50 vs. 47.22%) and high CD3E (68.75% vs. 54.17%) in MIBC patients defined as immunotherapy responders in GSE176307 dataset, which was consistent with the tendency in the IMvigor210 cohort (LCK: 61.76 vs. 46.52%, CD3E: 55.88 vs. 48.26%, Figures 8C, D). More importantly, the percentage of responders distributed in high-LCK and high-CD3E groups was obviously higher, considering the generally relatively limited response rate of immunotherapy, than the low-expression groups, both in GSE176307 (LCK: 22.73 vs. 13.64%, CD3E: 22.00 vs. 13.16%, Figures 8E,F) and IMvigor210 cohorts (LCK: 28.19vs.vs 17.45%, CD3E: 25.50 vs. 20.13%, immune cells), excluded (immune cells accumulated but not efficiently infiltrated) and desert phenotypes (absence of CD8⁺ T cells) (Figures 8I,J). Similarly, high-LCK and high-CD3E group also had more inflamed phenotypes (LCK: 42.28vs.vs 8.26%, CD3E: 42.74vs.vs 7.50%, Figures 8K, 8L). Supplementary Figure S3 presents a tendency of higher expression of LCK and CD3E in the responder group versus the non-responders (Supplementary Figures S3A–S3D), but statistical significance was not reached. However, a significant trend of LCK/CD3E expression was observed across the immune phenotypes (Supplementary Figure S3E–S3F). Last, the prognostic value of LCK and CD3E was validated by indicating higher expression of LCK (GSE176307: HR 0.44, 95% CI 0.26–0.75, p = 0.003; IMvigor210: HR 0.45, 95% CI 0.26–0.76, p = 0.003, Figures 8M,N) or CD3E (GSE176307: HR 0.64, 95% CI 0.48–0.86, p = 0.003; IMvigor210: HR 0.58, 95% CI 0.39–0.87, p = 0.008, Figures 8O,P) and was associated with an improved overall survival for bladder cancer patients receiving immunotherapy.