We obtained gene expression profiles and clinical data of MIBC patients from The Cancer Genome Atlas (TCGA) database, including 391 tumor samples, after excluding non-MIBC cases and those with incomplete data on age, recurrence-free survival (RFS), and overall survival (OS). Additional gene expression profiles and clinical data were retrieved from the Gene Expression Omnibus (GEO) database, focusing exclusively on MIBC patients. RNA-seq data were normalized to enhance comparability with microarray data. Three GEO datasets (GSE13507, n = 61; GSE31684, n = 78; GSE48075, n = 73) were merged to create a meta-GEO MIBC cohort (n = 212) for validating the prognostic value of CXCR6. For patients receiving anti-PD-L1 immunotherapy (IMvigor210), gene expression profiles and clinical data were obtained from the IMvigor210CoreBiologies repository.

CXCR6 expression was analyzed using the UALCAN platform (https://ualcan.path.uab.edu) and GEPIA2 (http://gepia2.cancer-pku.cn/). These tools were utilized to compare CXCR6 expression between bladder cancer and normal tissues, as well as across different cancer stages.

Pan-cancer survival analysis for CXCR6 was conducted using the TIMER2.0 platform (http://timer.cistrome.org/). Additionally, Kaplan–Meier curves were used to assess survival times. To standardize the cutoff value for survival analysis in both the TCGA and meta-GEO cohorts, CXCR6 expression levels were normalized (mean = 0, SD = 1). The ‘‘minimum P value’’ method was applied to identify the optimal cutoff for CXCR6 expression, ensuring the best stratification of patients based on OS.

GSEA was performed to investigate the potential mechanisms of CXCR6 in MIBC by identifying statistically significant gene sets between high and low CXCR6 expression groups. GSEA software v3.0 was used for the analysis, utilizing gene sets from the Molecular Signatures Database (MSigDB).

For immune microenvironment analysis, MIBC tissues were stratified into two groups based on CXCR6 expression levels. The CIBERSORT algorithm was applied to the TCGA dataset to assess the relationship between CXCR6 expression and 22 tumor-infiltrating immune cell subsets. Additionally, the correlation between CXCR6 expression and immune cell infiltration was further analyzed using TIMER2.0.

The scRNA-seq data from nine bladder cancer samples were downloaded from the GEO database under GSE222315 (16). Batch correction was performed using Harmony, and quality control along with clustering was conducted using Seurat, following standard preprocessing protocols (17). Principal component analysis (PCA) was used to reduce the dimensionality of the gene expression matrix, extracting the most essential features. Inter-cell distances were computed with the FindNeighbors function, which informed the construction of a shared-nearest neighbor graph. This was followed by clustering and cell grouping using the Louvain method. For data visualization, the PCA-reduced data was projected using UMAP (Uniform Manifold Approximation and Projection). To identify feature genes, the FindAllMarkers function was applied, and significant genes were determined using the Wilcoxon rank-sum test. Marker genes identified from the feature set were then utilized to assign cell type identities, referencing the CellMarker database (18) and relevant literature (19–24). Intercellular interactions were analyzed with CellChat, a tool for studying ligand-receptor signaling in specific pathways (25). To investigate the pathways and functions associated with CXCR6+ T/NK cells and CXCL16+ myeloid cells, gene enrichment analysis was conducted for each subcluster. Additionally, the scRNA-seq databases used to validate CXCR6 and CXCL16 expression were obtained from TISCH2 (http://tisch.comp-genomics.org).

This study included 131 MIBC patients (pT2-T4, N0-N3, M0) who underwent radical cystectomy at Fudan University Shanghai Cancer Center (FUSCC) between 2008 and 2012. All patients provided informed consent and had not received neoadjuvant chemotherapy. Tumor tissue microarrays were constructed, and single-staining IHC was performed as previously described (26). The IHC score for CXCR6 was calculated by multiplying the staining intensity (0 = negative, 1 = weak, 2 = moderate, 3 = strong) by the percentage of positive cells (0-100%). Anti-CXCR6 antibody (1:300, Proteintech, #36380) was used for staining, with a median IHC score of 180 serving as the cutoff. CXCR6 expression levels below 180 were classified as low expression. Two pathologists independently reviewed the slides, and all tissue microarray samples were stained concurrently.

Pearson’s χ² test was used to compare categorical variables, while the t-test was employed for continuous variables. For nonparametric data, the Mann-Whitney U test (Wilcoxon Rank Sum Test) was applied to compare outcomes between two independent groups. Survival curves were analyzed using the Kaplan-Meier method and compared with log-rank tests. Univariate and multivariate analyses were performed using the Cox regression model, with variables significant in the univariate analysis included in the multivariate model. A prognostic nomogram and calibration plots were created using R software (v3.0.2) with the ‘rms’ package. All statistical analyses were carried out using SPSS v21.0 (IBM Corp), and a two-sided P-value <0.05 was considered statistically significant.

First, we investigated the expression of CXCR6 in bladder cancer. Using the GEPIA2 and UALCAN online platforms, we found that CXCR6 expression was significantly lower in bladder cancer tumor tissues compared to adjacent normal tissues ( Figures 1A, B ). Moreover, CXCR6 expression decreased with advancing tumor stage ( Figure 1C ).

Next, we explored the prognostic significance of CXCR6 in MIBC. A forest plot of overall survival (OS) prognostic analysis across different cancer types using TISCH2 revealed that CXCR6’s predictive value varied among cancers ( Figure 1D ). Specifically, in bladder cancer, low CXCR6 expression was associated with worse patient outcomes. In particular, we assessed the prognostic value of CXCR6 in MIBC using the TCGA dataset. Applying the “minimum P value” approach to define the cutoff, 264 patients (67.5%) were classified as having low CXCR6 expression, while 127 patients (32.5%) had high expression. The characteristics of patients in the high and low CXCR6 expression groups were summarized in Supplementary Table 1 . Low CXCR6 expression was significantly associated with higher tumor grade. Kaplan–Meier analysis further demonstrated that low CXCR6 expression was significantly correlated with worse OS and RFS in the TCGA cohort ( Figures 1E, F ). This finding was validated in the meta-GEO database, where patients with high CXCR6 expression showed a trend toward better OS compared to those with low expression ( Figure 1G ). Additionally, IHC analysis in a 131-patient FUSCC cohort further confirmed that low CXCR6 expression was associated with poorer OS ( Figures 1H, I ). Cox regression analyses were performed to evaluate whether CXCR6 expression was an independent predictor of survival outcomes. After adjusting for confounding factors such as age, gender, tumor weight, stage, and necrosis, multivariate analysis showed that CXCR6 expression was an independent prognostic factor for both OS and RFS ( Supplementary Table 1 ).

To provide a practical tool for prognosis, we developed a nomogram to predict 3- and 5-year OS following surgery ( Figure 1J ). The model incorporated age at surgery, tumor stage, and CXCR6 expression as key predictors, with higher total scores indicating worse prognosis. Calibration plots showed that the nomogram accurately predicted 3- and 5-year survival, closely aligning with the ideal model ( Figures 1K, L ). In summary, high CXCR6 expression is a favorable prognostic marker in MIBC patients.

To further explore the molecular mechanisms underlying survival outcomes, we conducted a biological process and pathway analysis of MIBC from the TCGA dataset using GSEA. The results revealed that immune activation-related pathways were significantly enriched in tumors with high CXCR6 expression, including adaptive immune response, interferon gamma response, cell activation, T cell activation, leukocyte adhesion, and lymphocyte activation ( Figure 2A ; Supplementary Table 2 ). These findings are consistent with previous studies that emphasize the critical role of CXCR6 in regulating tumor immunity.

To further investigate the correlation between CXCR6 expression and immune cell infiltration, we utilized CIBERSORT to analyze CXCR6 expression in 391 MIBC samples and assess the proportions of 22 immune cell subtypes. The analysis revealed that several immune cell types, including naïve B cells, plasma cells, naïve CD4+ T cells, activated CD4+ memory T cells, activated NK cells, monocytes, macrophages (M0, M1, M2), and activated dendritic cells (DCs), were significantly less abundant in the low CXCR6 expression group compared to the high expression group ( Figure 2B ). To further validate the correlation between CXCR6 expression and immune cell infiltration, we used TIMER2.0. Consistent with previous findings, CXCR6 expression was significantly positively correlated with CD8+ T cells, CD4+ T cells, neutrophils, and dendritic cells (DC) ( Figure 2C ). These results suggest that CXCR6 plays a key role in promoting immune activation in MIBC.

To explore the role of CXCR6 in cellular composition and molecular profiles, we analyzed scRNA-seq data from 9 bladder cancer samples (GSE222315) (16). After integrating the transcriptional data, we identified six major cell types based on known cell-type-specific markers: epithelial cells, T/NK cells, fibroblasts, myeloid cells, endothelial cells, and B cells ( Figures 3A–C ). Notably, T/NK cells exhibited the highest levels of CXCR6 expression among these clusters ( Figures 3D, E ). Similar results were obtained from the public single-cell bladder cancer dataset (GSE14952) analyzed using TISCH2, further confirming the pivotal role of this receptor in these immune subsets ( Supplementary Figure 1A ).

Additionally, cell communication analysis using receptor-ligand interaction networks revealed that T/NK cells were connected with various cell types, including myeloid cells ( Figures 3F, G ). A closer examination indicated that T/NK-myeloid interaction was primarily mediated by the CXCL16-CXCR6 axis ( Figures 3H, I ). Notably, myeloid cells displayed the highest levels of CXCL16 expression ( Supplementary Figure 1B ), suggesting that CXCR6+ T/NK cells are actively recruited and influenced by CXCL16+ myeloid cells. This was further supported by similar findings from the TISCH2 analysis of the GSE145281 dataset ( Supplementary Figure 1C ). These findings suggest that CXCR6 is predominantly expressed in T/NK cells and plays a critical role in promoting T/NK-myeloid interaction through the CXCL16-CXCR6 axis.

Next, we further divided T/NK cells into four subclusters: T effector, T regulatory (T reg), T proliferative, and NK cells ( Figure 4A ; Supplementary Figure 2A ). We found that CXCR6 was predominantly expressed in T effector cells and NK cells ( Figures 4B, C ). Similarly, we subdivided myeloid cells into four subclusters: monocytes, macrophages, and DC1 and DC2 ( Figure 4D ; Supplementary Figure 2B ). CXCL16 was mainly expressed in macrophages and DC1 cells ( Figures 4E, F ).

Subsequently, we performed pathway enrichment analysis on CXCR6+ T cell, CXCR6+ NK cells, CXCL16+ macrophages, and CXCL16+ DC1 cells to further elucidate the functional roles of these cell subclusters. The results showed that CXCR6+ T cells were significantly enriched in pathways related to antigen receptor-mediated signaling, positive regulation of T cell activation, T cell receptor signaling, regulation of response to cytokine stimulus, and T cell-mediated cytotoxicity ( Figure 4G ). For CXCR6+ NK cells, the analysis highlighted the enrichment of pathways involved in natural killer cell-mediated cytotoxicity, T cell receptor signaling, and antigen processing and presentation ( Figure 4H ). In CXCL16+ macrophages, pathways related to cytokine-cytokine receptor interaction, chemokine signaling, positive regulation of phagocytosis, and regulation of macrophage activation were significantly enriched ( Figure 4I ). Finally, CXCL16+ DC1 cells were enriched in pathways related to antigen processing and presentation, regulation of T cell activation, positive regulation of T cell activation, cell activation involved in immune response, and dendritic cell antigen processing and presentation ( Figure 4J ). Taken together, these findings suggest that the CXCL16-CXCR6 axis mediates important immune interactions between T/NK cells and myeloid cells, particularly macrophages and DC1 cells, which are essential for orchestrating a robust anti-tumor immune response (27–30).

Given the earlier findings that CXCR6 promoted anti-tumor immunity in MIBC, we investigated its potential role in predicting immunotherapy response. First, we conducted a correlation analysis between CXCR6 and previously reported immunotherapy biomarkers using GEPIA2. The results revealed a significant positive correlation between CXCR6 expression and CD8A ( Figure 5A ), effector T cell signature ( Figure 5B ), and CD274 (PD-L1) ( Figure 5C ) in bladder cancer, suggesting that CXCR6 may hold potential as a predictive biomarker for immunotherapy response.

To further investigate the association between CXCR6 and immunotherapy outcomes, we analyzed the IMvigor210 cohort, a well-established phase II clinical trial evaluating the efficacy of anti-PD-L1 (atezolizumab) immunotherapy in advanced MIBC patients (31, 32) ( Figure 5D ). Our analysis revealed that patients with high CXCR6 expression demonstrated greater disease control ( Figure 5E ) and higher complete response rates ( Figure 5F ) compared to those with low CXCR6 expression. Additionally, lower CXCR6 expression was associated with poorer treatment responses ( Figure 5G ), whereas higher CXCR6 expression correlated with improved overall survival (OS) following immunotherapy ( Figure 5H ). The immune phenotypes of tumors in the IMvigor210 cohort have been detected, so we investigated the associations of CXCR6 expression among different immune phenotypes. We found that high CXCR6 expression was significantly associated with the inflamed immune phenotype ( Figure 5I ), which is known to respond more favorably to immunotherapy, reinforcing CXCR6 as a potential biomarker for predicting immunotherapy efficacy.