We obtained the RNA-seq datasets for CD56^bright NK cells and CD56^dim NK cells along with 3 stromal cell types (i.e., epithelial, endothelial, and fibroblast cells), and 18 other immune cell types (i.e., mast cells, memory B cells, naïve B cells, eosinophil, monocyte, neutrophil, γδ-T cells, immature myeloid dendritic cells, mature myeloid dendritic cells, M1-macrophages, M2-macrophages, naïve CD8⁺ T-cells, helper T-cells, regulatory T-cells, central CD4⁺ memory T-cells, effector CD4⁺ memory T-cells, central CD8⁺ memory T-cells, and effector CD8⁺ memory T-cells) from the Human Bulk Cell-type Catalogue (HBCC) (25) database. To adjust the gene-transcript abundance for the intra and inter-dataset variability, we used a Bayesian multilevel noise modelling approach, Cellsig (25). Cellsig used the cellular differentiation hierarchy to compute the missing transcript-abundances for every node in the hierarchy ( Supplementary Figure S1 ).

We used the adjusted cell type transcriptomes from Cellsig to generate the transcriptional signature matrix ( Supplementary Table S1 ) using CIBERSORTx (26); with G.min = 25, G.max = 30, and q.value = 0.01, while the rest of the parameters were kept default. To assess the specificity of the transcriptional signature, we performed principal component analysis (PCA) on the HBCC dataset with the identified cell type markers.

We obtained the bulk RNA-seq datasets and the matched clinical information of BLCA patients from the cancer genome atlas (TCGA) through the GDC Data Portal (27). We removed any duplicated sample-transcript pairs within the retrieved data, then scale-normalised the transcript-abundances for the sequencing depths by trimmed mean of M values (TMM) method (28). We selected only the transcriptomes of primary tumour (n = 408) and adjacent normal tissue samples (n = 19) for the analysis. To estimate the relative abundances of the tumour-infiltrating immune cells, we applied the cellular deconvolution method on the patient transcriptomes with the CIBERSORT deconvolution algorithm (29) implemented in the tidybulk R-package (30). We used our generated cell type signature matrix as the deconvolution reference, while the other parameters were kept default.

To correlate the cell type signature abundances with the patient outcomes, we carried out survival analysis with Kaplan–Meier (KM) estimates with survminer (31) package, where the patients were split into groups (i.e., high, and low) based on the median abundance of a cell type abundance. The differences between the survivals of the groups were compared with log-rank (Mantel–Cox) test (32). P-values of the KM curves were adjusted following the Benjamini-Hochberg’s (BH) method ( Supplementary Table S2 ).

We conducted differential gene expression analysis to contrast between different patient groups with the DESeq2 (33) package implemented in tidybulk (30). Additionally, the p-values of the contrasts were multiple-test corrected with Benjamini-Hochberg’s (BH) false discovery rate (FDR) method (34). After the contrast, we only selected the genes with |log₂FoldChange| ≥ 1 and adjusted p-value < 0.05 as significantly differentially expressed (DE).

Pathways associated with a target gene set were identified using with functional overrepresentation analysis (ORA) (35) by clusterProfiler (36) package. We performed the analysis against various pathway modules, such as gene ontology biological process (GOBP) (37), KEGG (38), Reactome (39), Wikipathways (40), and BioPlanet pathway database (41). Enriched pathways with p-values < 0.05 were considered significant.

We also obtained an independent single-cell RNAseq dataset (GSE135337) (42) from the Gene Expression Omnibus (GEO) database (43). This dataset recorded the single cell transcriptome from 7 bladder cancer patients with varying tumour stages and grades. The patients were newly diagnosed and did not receive any prior treatments. We analysed this dataset with the Seurat (44) package. While preprocessing, we kept only the high-quality cells that satisfy the criteria of nFeature_RNA > 200, nFeature_RNA < 5000, and percentage of mitochondrial transcripts < 10. Also, we removed the ribosomal transcripts from the count information. All samples were pre-processed and then normalised by the SCTransform method with v2 regularization separately before we merged them together with SCTransform. The top 15000 features per sample were used to identify the integration anchor set for the merging. After the integration of the samples, principal component analysis (PCA) was performed up to the first 15 dimensions. To cluster the cells, first the neighbour cells were searched through the shared nearest-neighbor (SNN) graph approach with default parameters. Then the cells were clustered through the Leiden algorithm in igraph method with a clustering resolution of either 0.2 or 0.6. After the clustering, we identified the cluster-specific markers by the Wilcoxon Rank Sum test with a log-FoldChange threshold of 0.25 ( Supplementary Table S3 ). In order to identify the characteristics of the clusters, SingleR (45) package was used where the cell type marker references were taken from the EncodeBlueprint database. To assess the presence of CD56^bright and CD56^dim NK signatures, the clusters were screened using the AddModuleScore function. To obtain the ligand-receptor interactions and estimate the communications between the clusters, cellchat (46) package was utilised.

Statistical significance of the comparison between the means between two groups in the boxplots were achieved through the non-parametric Wilcoxon signed-rank test (47), while comparison of means between multiple groups was performed by non-parametric ANOVA (Kruskal–Wallis) test (48).

NK cells are thought to play a significant role in cancer immune surveillance through direct cytotoxicity of tumour cells and immunomodulation (49, 50). We hypothesized that subset-specific functionalities of NK cells may be present in the BLCA TME that could regulate anti-tumour immunity. To address this, we first sought out the core transcriptional markers of CD56^bright and CD56^dim NK cells that segregate these phenotypic subsets from each other, and from 21 other stromal and immune cell types, including T lymphocyte subsets ( Figures 1A, B ). In the PCA analysis with the generated transcriptional signatures, a clear contrast between stromal and immune cell subsets was observed ( Figure 1A ). Within immune cells, naïve and memory B cells clustered more closely with the myeloid cell subsets, whereas the NK cell subsets clustered more closely with T lymphocyte subsets. Importantly, the NK cell subsets formed distinct clusters separated from each other and from all other immune cell subsets, including the CD8⁺ T cell subsets ( Figure 1A ). Moreover, the classical gene markers for CD56^dim and CD56^bright NK cells showed preferential subset-specific expression patterns, respectively ( Figure 1B ). These results show that we have identified gene signatures specific for CD56^dim and CD56^bright NK cell subsets that can be distinguished from each other and from other immune and stromal cell-types.

We next estimated the relative abundances of the tumour-infiltrating NK subsets in the TCGA-BLCA patient cohort by cellular deconvolution. Bladder cancers are primarily of epithelial origin (51). Indeed, we found epithelial cells to be the most abundant stromal cells in the BLCA TME followed by fibroblasts ( Figure 1C ). Amongst, immune cell types, signatures of immature dendritic cells, memory B cells, and CD8⁺ central memory T cells were the most abundant in BLCA ( Figure 1D ). Interestingly, we observed the CD56^bright NK cell signature to be the major infiltrating NK subset in BLCA patients ( Figures 1E, F ). Furthermore, patients from different histopathological grades and stages also harboured similar NK subset infiltration profiles ( Supplementary Figures S2A–E ). Intriguingly, patients with more advanced tumour-stages (e.g., higher grades, and higher pathologic stages) showed significantly higher abundance of infiltrating CD56^bright NK cells, but not CD56^dim NK cells. Additionally, a similar enrichment of CD56^bright NK signatures was present in the BLCA tumour-adjacent normal tissues ( Supplementary Figure S3A ). These findings indicate that CD56^bright NK cells may preferentially infiltrate BLCA tumours compared to CD56^dim NK cells.

We also analysed an independent single-cell RNAseq dataset which includes transcriptomes of 7 muscle-invasive bladder cancer patients. We found that signatures for epithelial cells followed by fibroblasts were most abundant within the tumour, consistent with our findings from the deconvolution of TCGA-BLCA dataset ( Supplementary Figure S3B ). We then selected only the CD45⁺ cells and performed clustering to identify the tumour-infiltrating immune cell populations. From this, we obtained 8 clusters representing different immune cell populations ( Supplementary Figure S3C ). Upon annotating the phenotypic characteristics of the clusters, cluster 6 bore the highest resemblance to CD8⁺ T and NK cells ( Supplementary Figure S3D ). Additionally, amongst all the clusters, cluster 6 contained the phenotypic signatures of both NK subsets ( Supplementary Figure S3E ). We also detected that the presence of NK subset signatures in cluster 6 was significantly enriched compared to the average of background signatures for other immune cells used in the deconvolution of the TCGA bulk transcriptome data ( Supplementary Figures S3F, G ).

We next assessed the potential implications of the infiltrating NK subsets for BLCA patient prognosis, which revealed an association between a higher abundance of CD56^bright NK cells, but not CD56^dim NK cells, with improved BLCA patient survival ( Figure 1G ). Also, CD56^bright NK cells were predicting a trend, although not significant, towards the favourable patient survival even after associating other important clinical variables such as high mutational burden and CD8⁺ T cell abundance ( Supplementary Figure S4 ). Additionally, the favourable prognostic implications from CD56^bright NK cell infiltration can greatly differ between patients with varying clinicopathologic tumour status. Interestingly, this positive prognostic association of CD56^bright NK cells was more prominent in patients with moderately developed non-metastatic (e.g., pathologic T2, T3, N0, M0, stage II, and stage III) tumours; as opposed to the patients with more advanced and metastatic, suggesting a potential anti-tumour role of CD56^bright NK cells in the BLCA TME ( Supplementary Figures S5A–D ).

Activation of NK cell functionality is dictated by the signaling through activating and inhibitory surface receptors. However, which specific NK family receptor(s) drive NK cell anti-tumour responses in human cancer patients is understudied. Since the signature of CD56^bright NK cells is correlated with favourable BLCA patient survival, we next determined whether expression of NK family receptors is similarly important. Tumour-specific expression of the NK receptors NCAM1, IL2RG, IL18R1, CXCR3, KIR2DL4, KLRD1, NCR1, KIT, CD7, CD2, and CD226 showed a positive correlation with the infiltrating CD56^bright NK subset, but not with CD56^dim NK cells ( Figure 2A ). Intriguingly, many of these receptors were also positively correlated with the signature of CD8⁺ effector memory T cells (CD8⁺ T_EM) ( Figure 2A ). The expression of KLRK1 (NKG2D), NCR2 (NKp44), CD2, and CD160 receptor genes were associated with favourable prognosis ( Figures 2B–E ). Moreover, patients harbouring higher expression of these receptors alongside higher CD56^bright NK abundance showed trends towards better survival compared to others ( Figures 2B–E ). However, involvement of these receptors with the CD56^dim NK subset was not conclusive ( Figures 2B–E ). Multivariate analysis also highlighted the prognostic associations of these receptors with the signatures of CD56^bright NK cell ( Supplementary Figure S6A ). Similar survival trends for KLRK1 (NKG2D), CD2, and CD160 expression were also predicted in combination with the CD8⁺ T_EM cells ( Supplementary Figure S6B ).

Next, we asked whether the ligands for NK family receptors are also positively associated with the infiltrating immune cell signatures. Correlation analysis indicated positive associations of CD160 ligands such as TNFRSF14 and HLA-C, and NKG2D (KLRK1) ligands MICA and MICB with the tumour-infiltrating CD56^bright NK and CD8⁺ T_EM cells ( Supplementary Figure S6C ). Previously, NKG2D (KLRK1) ligands MICA, and MICB were shown to have associations with beneficial BLCA patient prognosis (24). Expression of the CD160 ligands trended towards favorable prognosis, whereas the ligands of NKp44 (NCR2) were associated with improved survival ( Supplementary Figure S6D ). Compared to the BLCA tumour-adjacent normal tissues, tumour expression of transcripts encoding CD2, KLRK1, and NCR2 were higher ( Supplementary Figure S6E ). Also, in the analysed BLCA scRNA-seq data, most of the NK cell markers were expressed by the cells of cluster 6 to some extent; except for CD44, expression which was also observed in other clusters ( Supplementary Figure S6F ). These results highlight that signaling through NKG2D (KLRK1), NKp44 (NCR2), and CD160 receptors may be critical in the downstream activation of tumour infiltrating CD56^bright NK and/or CD8⁺ T_EM cells, which may contribute to a more favourable BLCA prognosis.

The anti-tumour functions of NK cells may be directed by other cells in the TME. Conversely, NK cells can also influence the recruitment and functional activation of other tumour-infiltrating leukocytes (52). We next sought to determine the immune cell types that may be critical to BLCA patient prognosis, and their potential association with the CD56^bright and CD56^dim NK subsets.

Pair-wise correlation between the signature abundances of the immune and stromal cell subsets revealed an intriguing inverse correlation between the two NK subsets ( Figure 3A ). Additionally, the signature of CD56^bright NK cells was positively associated with immature and mature myeloid dendritic cells (imDCs and mDCs), Memory B, M1 macrophage, Naïve CD8⁺, and CD8⁺ T_EM cells. In contrast, CD56^dim NK was positively correlated with infiltrating central CD8⁺ memory T cells (CD8⁺ T_CM), but negatively with CD8⁺ T_EM and T helper cells ( Figure 3A ). This finding suggests that the infiltrating CD56^bright, and CD56^dim NK cells are associated with distinct immune signaling networks within the tumour.

We next asked whether the associations between the NK subsets and other cell types might influence patient survival. First, we determined the prognostic association of each cell type separately and observed favourable survival associations for CD8⁺ T_EM and mDCs, and poor survival association for fibroblast cells ( Figures 3B–D ). Next, we investigated the survival associations of these cells in combination with the NK subsets and found that BLCA patients with high CD56^bright or CD56^dim NK with high CD8⁺ T_EM and mDCs signatures had improved survival ( Figures 3B, C ). Cancer-associated fibroblasts (CAFs) can potentially dampen the anti-tumour functions of tumour-infiltrating NK cells (53). We observed that patients with high fibroblasts but low CD56^bright NK had the worst survival outcome, with this effect not as pronounced for CD56^dim NK cells ( Figure 3D ). This highlights that the protective effect that correlates with these immune signatures are tumour environment -dependent. To support this finding, we estimated each of the cell type signatures using gene set enrichment approach independently and mapped it onto the BLCA patient clusters that was generated from the estimated abundances from the deconvolution. Interestingly, we observed that a group of BLCA patients were highly enriched for the CD56^bright NK signature as well as for CD8⁺ T_EM and mDCs signatures. Also, a collinearity was present between the abundances of these cell signatures in the BLCA patients ( Supplementary Figures S7A–E ). These further supports the putative NK-mDC-CD8⁺ T signalling axis in BLCA tumour microenvironment. We also sought to detect the putative immune crosstalk in the analysed BLCA scRNA-seq data. Despite many potential ligand-receptor interactions between cluster 6 cells with the immune cells from other clusters, communication between the cells was slightly ambiguous due to the limited number of immune cells in the data ( Supplementary Figures S3H–J ).

In the TME, complex interactions between immune cells can occur through a wide-variety of secreted factors (52). Considering our identified NK-mDC-CD8⁺ T signaling axis, we extended our analysis to encompass the effect of genes encoding secreted factors in the BLCA TME. We detected a strong positive correlation between the expression of cytokine-chemokines, and their receptor encoding genes and myeloid dendritic cells, CD56^bright NK and CD8⁺ T_EM cell phenotypes ( Figure 3E ). Conversely, expression of fibroblast-associated factors, such as VEGFB, VEGFC, TGFB1, TGFB2, TGFB3, and IL6 were positively associated with the fibroblast signature, but were negatively correlated with the infiltrating myeloid dendritic cell, NK, and T-cell subset signatures ( Figure 3F ); indicating that cancer-associated fibroblasts might be negatively modulating the anti-tumour immune cell functions (54). Overall, these results highlight the presence of complex anti- and pro-tumorigenic interactions in the BLCA TME that may shape the anti-tumour function of the CD56^bright NK-mDC-CD8⁺T cell axis.

To obtain insights into the functional mechanisms associated with the BLCA-infiltrating NK subsets, we sought to understand the molecular pathways associated with our CD56^bright and CD56^dim signatures. We performed functional pathway analysis with genes that are positively correlated with the abundance of the NK cell subsets. Various immune signaling pathways were positively associated with the CD56^bright NK signature, including pathways associated with NK cell-mediated responses, interferon signaling, T and B cell maturation and signaling, and cytokine-chemokine signaling ( Figure 4A ). In contrast, the CD56^dim NK signature was weakly associated with interferon and other cytokine-chemokine signaling, and dendritic cell recruitment pathways ( Figure 4B ). Surprisingly, several pro-tumour pathways such as TGF-β signaling, ECM remodelling, and IL-6 signaling were associated with CD56^dim NK cell abundance, suggesting these pathways suppress the functions of this NK subset.

We also sought out the differences in the functional responses between the patients with high and low/no NK cell infiltration. To achieve this, we split the patients into three groups (e.g., high, medium, and low) based on total NK signature abundance. Then we compared the transcriptomes of the high total NK against low total NK group through a differential gene expression analysis. The high total NK group was found to over-express key NK function associated genes such as SELL, IFNG, KLRC2 (NKG2C), KIR2DS4, CXCL10, ZNF683, and CXCL11 ( Figure 4C ). Remarkably, this patient group also expressed high amounts of NK family receptors such as KLRC1 (NKG2A), KIR2DL1, KIR2DL3, KIR2DL4, KIR3DL3 and HAVCR1 (TIM1) ( Figure 4C ). Pathway analysis with the differentially expressed genes revealed that immune signaling pathways (such as NK cell mediated cytotoxicity, antigen processing and presentation, IL-7 interactions in immune responses, Type II interferon signaling, and chemokine mediated signaling pathways) were mostly associated with the high NK infiltration group in contrast to the low NK infiltration group, which had upregulated expression of pro-tumour pathways such as fibroblast-associated pathways (e.g., collagen production), and ECM remodelling pathways ( Figure 4D ). Many of these pathways were also associated with the cells of NK/CD8⁺ T cell cluster (Cluster 6) of the analysed BLCA scRNA-seq dataset ( Supplementary Figure S8 ). These findings indicate that BLCA patients with low overall NK infiltration had more pro-tumorigenic responses as opposed to the patients with high overall NK abundances where more anti-tumorigenic pathways were present.

Transcription factors (TFs) can often play key roles in modulating the pro-/anti-tumour functions of immune cells (55–57). Therefore, we searched for the transcription factors that are positively correlated with the NK signatures and explored their possible impact on patient survival. Positive associations of several TF-encoding genes such as ZNF683 (HOBIT), STAT2, and IRF3 with the abundances of tumour-infiltrating CD56^bright NK and CD8⁺ T_EM cells were observed ( Figure 5A ). Transcript abundance of these TFs was significantly higher in the BLCA tumour tissues compared to the adjacent normal tissue ( Supplementary Figure S9A ). Additionally, these TFs were significantly associated with favourable patient survival, and patients with higher expression of these TFs alongside higher CD56^bright NK and CD8⁺ T_EM cells, but not CD56^dim NK cells, abundances had better survival probabilities ( Figure 5B , Supplementary Figure S9B ). Furthermore, multivariate analysis suggests the potential co-implications of ZNF683 and CD56^bright NK cell signature in favourable BLCA patient prognosis ( Supplementary Figure S9C ). These results identify TF signatures that may serve as markers of modulation of anti-tumour functions of tumour-infiltrating NK and T cell subsets in BLCA.

HOBIT has been shown to regulate a program of tissue-residency in CD8⁺ T lymphocytes (58). We hypothesized that expression of HOBIT might also dictate tissue-residency programming in tumour-infiltrating NK cell subsets. Interestingly, both the CD56^bright NK and CD8⁺ T_EM cell signatures were strongly correlated with the expression of tissue-residency markers such as ITGAE (CD103), and CXCR6 ( Figure 5C ), while CD69 expression was weakly correlated, and no correlation was observed for ITGA1 (CD49a) with CD56^bright NK and CD8⁺ T_EM cells ( Figure 5C ). Overall, these findings suggest the presence of a tissue-resident phenotype of NK and/or T cells in the BLCA TME.

Next, we detected numerous lymphocyte development and function associated pathways that are significantly associated with ZNF683 (HOBIT), STAT2, and IRF3 expression ( Figure 5D , Supplementary Figures S9D, E ), indicating that these TFs may regulate the priming, and their downstream functionalities of tumour-infiltrating NK and CD8⁺ T cells. Also, expression of the NK associated receptors was correlated with the expression of these TFs, indicating the potential influence of these TFs on the functional responses on NK and CD8⁺ T cell subsets ( Supplementary Figure S9F ). Also, expression of ZNF683 (HOBIT) as well as markers of tissue-residency were detected in a fraction of the cluster 6 population in the BLCA scRNA-seq dataset ( Supplementary Figure S9G ) supporting the presence of a tissue-resident NK or T cell population in the bladder TME.

Furthermore, we detected numerous genes known to be regulated by the HOBIT transcription factor that were co-expressed with ZNF683 (HOBIT) in the BLCA patients ( Figure 5E ). EOMES regulates the tissue-residency program of NK and CD8⁺ T cells alongside HOBIT (59, 60). We observed a positive correlation between ZNF683 (HOBIT) and tissue-retention markers such as CD69 and CD103 (ITGAE) that are essential for the establishment of tissue-resident phenotypes in NK and T cells (61–63). On the other hand, tissue-egress factors such as S1PR1 and CCR7 (64) were found to be negatively correlated with ZNF683 (HOBIT) expression. These results suggest that HOBIT may regulate TIL phenotype in the BLCA TME by driving a tissue-residency program in tumour-infiltrating lymphocytes.

Chronic activation/inhibition stimuli in the TME can lead to a dysfunctional/exhausted state in tumour-infiltrating NK and T cells, characterized by the expression of various immune checkpoint receptors (65, 66). Hence, we examined the potential association of immune checkpoint receptors with tumour infiltrating NK and T cell signatures in the BLCA TME.

Firstly, we compared the expression of several immune checkpoint receptor genes between the tumour and adjacent normal tissues. While comparing between matched and unmatched tumour and normal tissue samples, upregulation of checkpoint receptor genes, notably CD47, CTLA4, HAVCR2, PVRIG, SIGLEC7, and SIGLEC9, were specifically detected in malignant compared to normal tissues ( Figures 6A, B ). Furthermore, some of these checkpoint receptors were also positively correlated with tumour infiltrating CD56^bright NK and CD8⁺ T_EM cell signatures; intriguingly, a similar association with T_reg cells and checkpoint receptors was also observed ( Figure 6C ). We next explored the association between immune checkpoint receptor gene expression in patient survival in combination with tumour-infiltrating NK subsets. Combining the NK subset abundances with the expression of checkpoint receptor genes showed little or no significant impact in BLCA patient survival ( Figure 6D , Supplementary Figure S10A ). Conversely, there were significant associations with the expression of several checkpoint receptors in combination with the CD8⁺ T_EM cell signature and survival, with a marked inverse relationship between higher expression of KLRC1 (NKG2A), TIGIT, SIGLEC7, SIGLEC9, and HAVCR2 checkpoint receptors and the CD8⁺ T_EM signature on BLCA patient survival ( Figure 6D , Supplementary Figure S10A ). Also, most immune checkpoint receptor genes were expressed by cells in the NK/CD8⁺ T cluster of the BLCA scRNA-seq dataset ( Supplementary Figure S10B ). Collectively, these results suggest that expression of immune checkpoint receptors may negatively impact the prognostic implications of CD8⁺ T cells but not NK subsets in BLCA patients.