The analysis was done with R version 4.2.0. Basic data transformation tasks were accomplished with the tidyverse package bundle (18) and the development package trafo (https://github.com/PiotrTymoszuk/trafo). Exploratory data analysis and hypothesis testing was done with the rstatix (19), ExDA (https://github.com/PiotrTymoszuk/ExDA) and microViz (https://github.com/PiotrTymoszuk/microViz) packages. Network analysis was accomplished with igraph (20). Results were visualized with the packages ggplot2 (bubble plots) (21), ExDA (box plots), microViz (Forest plots, Volcano plots, bar plots of regulation estimates and p values), survminer (Kaplan-Meier plots) (22) and ggnetwork (network plots) (23).

The following publicly available RCC data sets were re-analyzed: TCGA Clear Cell Carcinoma (KIRC) project (24), whole-genome subsets of the CheckMate 010 (CM 010) and CheckMate 025 (CM 025) studies (25), GSE73731 (26), GSE167093 (27), RECA-EU (28) and E-MTAB 1980 (29). Author-provided expression data and clinical information for the TCGA cohort were extracted from the GDC Data Portal with the TCGA-Assembler script (https://github.com/compgenome365/TCGA-Assembler-2/blob/master/TCGA-Assembler/). Author-provided expression and clinical data for GSE73731 and GSE167093 were fetched with the GEOquery package (30). E-MTAB 1980 and RECA-EU data sets were imported from the ArrayExpress and ICGC Data Portal repositories, respectively, with in-house developed R scripts. To investigate possible effects of therapy with checkpoint inhibitors, the CheckMate 025 everolimus (CM 025 EVER) and nivolumab (CM 025 NIVO) treatment arms were analyzed separately.

For the microarray expression studies (GSE73731, GSE167093, E-MTAB 1980), integration of multiple probes targeting the same gene was accomplished by geometric mean. Expression values were transformed with log ₂(transcriptcount + 1) (RNA sequencing: TCGA, CM 010, CM 025, RECA-EU) of log ₂. Immune infiltration estimates were calculated using the QuanTIseq and xCell algorithms (immunedeconv package) (31–33).

Gene signatures corresponding to Reactome pathways were extracted from the MSig database, version 7.5.1, and the signatures’ single sample gene set enrichment analysis scores (ssGSEA) were calculated with the GSVA algorithm (34). Genes associated with T cell exhaustion were retrieved from four recent papers (35–38).

The chemokine genes of interest were identified within the set of 28 chemokine genes available for all investigated datasets (CCL2, CCL7, CCL8, CCL11, CCL13, CCL17, CCL20, CCL21, CCL22, CCL24, CCL25, CCL26, CCL28, CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL9, CXCL10, CXCL11, CXCL12, CXCL13, CXCL14, CXCL16, CXCL17, CX3CL1, XCL1). The chemokine genes were screened in the TCGA KIRC cohort for differences in expression levels between the RCC and normal kidney tissue (paired normal/tumor samples, two-tailed paired T test, function test_two_groups(), microViz package), and for correlation with CD8⁺ T cell content in RCC predicted with the QuanTIseq algorithm (31) (Spearman’s correlation, correlate_variables(), package ExDA).

Among 28 chemokine genes measured in all investigated collectives, expression of CXCL9, CXCL10, CXCL11, CXCL13 and XCL1 was found to be highly enriched in the malignant tissue (paired two-tailed T test) and strongly significantly associated with CD8⁺ T cell levels (Spearman’s correlation) in the TCGA KIRC cohort ( Supplementary Figure S1 ). These chemokine genes and their cognate receptors CXCR3 (CXCL9/10/11), CXCR5 (CXCL13) and XCR1 (XCL1) were investigated further in the current report.

log₂-transformed expression was compared between the donor-matched normal kidney and cancer samples by two-tailed paired T test and Cohen’s d effect size statistic. Correlation of log₂-transformed gene expression with QuantIseq- and xCell-predicted non-malignant cell content (31, 33) was investigated with Spearman’s test. The expression comparison and correlation results were corrected for multiple testing with the false discovery rate (FDR) method (39). The scaled Spearman’s correlation matrices for the genes of interest and cell types associated with them with moderate-to-large strength (correlation coefficient ρ > 0.4) were visualized as force-directed network plots (packages igraph and ggnetwork) (20, 23).

Correlation of log₂-transformed gene expression with overall (OS) or relapse-free survival (RFS) by uni-variable Cox proportional hazard modeling (40) including the linear and spline term for the gene expression variable. Significance for model terms was corrected with FDR for multiple testing. Linear term model estimates are presented as hazard ratios (HR) with 95% confidence intervals.

Multi-parameter modeling of OS was done by Ridge regularized Cox regression (package glmnet) (41). The explanatory factors included age, squared age, sex, tumor grade, tumor stage, metastasis stage) variables and log₂-transformed expression values of the genes of interest. Two models were trained in the TCGA KIRC cohort: a model with demographic/grade/stage variables only and a model with the complete set of explanatory factors and both models were subsequently validated in the E-MTAB 1980 and RECA-EU collectives. To assess the add-on effect of gene expression on the survival prediction accuracy as compared with the demographic/grade/stage information, concordance indexes of the full models and the demographic/grade/stage-only models were compared.

By PAM (partition around medoids) - cosine distance clustering (42, 43) in respect to to normalized log₂-transformed expression of the genes of interest, two clusters of RCC samples, the chemokine^high and chemokine^low cluster were defined in the training TCGA KIRC cohort. The clustering algorithm was chosen based on its excellent reproducibility in 10-fold cross validation (44) and good explanatory performance as compared with hierarchical and KMEANS clustering algorithms ( Supplementary Figure S2A ). Choice of the cluster number (k = 2) was motivated by the bend of the within-cluster sum of squares curve, the peak mean silhouette statistic and visual inspection of the heat map of distances between cancer samples ( Supplementary Figures S2B, C ).

Prediction of the chemokine cluster assignment in the test collectives (E-MTAB 1980, GSE73731, GSE167093, RECA-EU and CheckMate cohorts) was done with a 15-nearest neighbor (15-NN) classifier with inverse distance weighting. The prediction yielded clustering structures with comparable fractions of explained variance and comparable cluster sizes in the training and the test cohorts indicative of high reproducibility ( Supplementary Figure S3A ).

xCell estimates of non-malignant cell infiltration (43) were compared between the chemokine clusters by FDR-corrected Mann-Whitney U test. Differences in distribution of tumor stages, MSKCC risk strata and frequencies of best overall therapy response (complete/partial response vs. stable/progressive disease) between the clusters were compared by FDR-corrected χ ² test. Differences in OS and RFS between the two clusters were assessed by Peto-Peto test (22).

Differences in ssGSEA scores (34) of Reactome pathway gene signatures between the chemokine clusters were investigated by FDR-adjusted two-tailed T test ( Supplementary Table S1 ). Genes differentially expressed between the chemokine clusters were identified by FDR-corrected two-tailed T test and > 1.25 fold-regulation in the chemokine high versus chemokine low cluster ( Supplementary Tables S2, S3 ). Modulation of KEGG-listed signaling pathways in chemokine^high versus chemokine^low cluster cancers based on the differential gene expression was investigated with SPIA (45) ( Supplementary Table S4 ).

Rules of assignment of genes to biochemical reactions were retrieved from the Recon2 human metabolism model available via the BiGG database (46). Regulation of biochemical reactions between the chemokine^high and chemokine^low clusters based on whole-genome estimates of differential gene expression and their standard errors was done by evaluation of the gene assignment rules with the ‘gene - protein - reaction’ (GPR) principle (47). Standard deviation, 95% confidence intervals and p values for the predicted reaction regulation estimates were obtained by a Monte Carlo simulation (n = 1000 draws) ( Supplementary Table S5 ). The analysis was done with the package BiGGR (47) and the development package biggrExtra (https://github.com/PiotrTymoszuk/biggrExtra).

Peripheral blood mononuclear cells, adjacent normal kidney tissues, tissues from the center and periphery of ccRCC tumors ( Supplementary Table S6 ; n=4) were obtained freshly at the day of surgery. Patients were all treated at the University Hospital Regensburg (Urology) between 2019 and 2020. Tissues were processed as described previously (48). Peripheral blood was enriched using leukocyte reduction system cones and processed using Ficoll density gradient centrifugation. After processing, single-cell suspensions were cryopreserved until the day of analysis. For flow cytometric analyses, cells were permeabilized using the Cytofix/Cytoperm kit (BD) and stained with antibodies against: CD3, CD14, CD56, CD19, CD25, CD4 (all BD); CD8 (BioLegend), CXCL9, (Biotechne); and with Fixable Viability Dye eFluor 708 (eBioscience). Flow cytometry was performed using LSRFortessa (BD).

The respective dataset consisting of samples obtained from 12 patients with RCC was downloaded as AnnData object (h5ad) from previously published studies of Li R. et al. (49) [Dataset (50)] and imported in Scanpy version 1.9.1. The dataset was controlled for the quality with scanpy by thresholding the number of detected genes (200), counts (2000) and the fraction of mitochondrial reads (<30%).

Cell transcriptomes were embedded into a batch-corrected low-dimensional latent space using scVI (51, 52) treating each sample as a separate batch. The scVI model was trained on the 2000 most ‘highly variable genes’ as determined with scanpy’s “pp.highly_variable_genes’ with parameters ‘flavor=“seurat’, and batch_key=‘orig.ident’. A neighborhood graph and UMAP embedding was computed based on the scVI latent space. All cell-type annotations were used from the original study (49). Annotated cell types were confirmed by a set of cell type–specific markers such as, CD3E, CD68, CD8A, CD4, CD79A, KIT, CDH5, ACTA2, EPCAM. For more detailed analysis ‘tumor’ and ‘normal-kidney’ samples were extracted from the dataset and merged into a sparate AnnData object.

Data analysis and graphical visualization was performed with scanpy v.1.9.1, anndata v.0.8.0, umap v.0.5.3, numpyv.1.21.5, scipy v.1.7.3, pandas v.1.4.2, scikit-learn v.1.02.2, statsmodels v.0.13.2, pynndescent v.0.5.7, and python-igraph v.0.10.2.  A method and samples overview of scRNA dataset is shown in Supplementary Figure S4 .

Altogether, we analyzed clinical and transcriptome data from 8 publicly available RCC cohorts including a total of 1819 RCC samples. Males constituted between 57 and 76% of investigated patients and the gender distribution was comparable between the cohorts. The median age ranged between 60 and 64 years and was similar in the study collectives. The RECA-EU, GSE167093 and E-MTAB 1980 cohort individuals tended towards lower tumor grades (G2) as compared with the remaining collectives with the grade data available. The majority of tumors was stage T1. In the TCGA, RECA-EU and E-MTAB 1980, information on initial metastasis status was available; only a small fraction of participants had lymph node or distant metastases at RCC diagnosis. The median observation time was clearly shorter in the CM 010 and CM 025 collectives. In these cohorts, the rate of relapses and mortality within the observation time window was substantially higher. In the samples from the CM studies, between 5.1 (CM 025 everolimus) and 25% patients (CM 010 and CM 025, nivolumab) displayed an overall therapy response defined as CR or PR ( Table 1 ).

First, we sought to identify RCC-specific chemokines which may mediate CD8⁺ T cells recruitment. Expression of 28 chemokine genes was compared between RCC and normal kidney, and correlated with CD8⁺ T cell infiltration predicted by the QuanTIseq algorithm (46) in the TCGA KIRC cohort. CXCL9/10/11, CXCL13 and XCL1 were found to be at least four-fold upregulated in the cancer tissue and correlated strongly with CD8⁺ T cell infiltration with ρ > 0.6 in Spearman test ( Supplementary Figure S1 ). These chemokine genes along with genes coding for the cognate receptors CXCR3 (CXCL9/10/11), CXCR5 (CXCL13) and XCR1 (XCL1) were analyzed further.

CXCL9/10/13 and XCL1 chemokine transcripts and the CXCR3, CXCR5 receptor transcripts were significantly enriched in RCC as compared with the non-malignant tissue in an analysis of donor-matched samples in the TCGA KIRC and GSE167093. This upregulation was particularly strong for CXCL9, CXCL10 and CXCL13. In turn, significantly elevated expression of CXCL11 and XCR1 in RCC as compared with normal kidney tissue could be observed only in the TCGA KIRC cohort but not in the GSE167093 collective ( Figure 1 ).

As revealed by network analysis of Spearman’s correlations between gene expression and non-malignant cell content predicted by the xCell algorithm, expression of CXCL9/10/11, CXCR3, CXCL13, XCL1 and XCR1 was tightly associated with increased infiltration by panCD8⁺, central memory and effector memory CD8⁺ T cells, activated myeloid and plasmacytoid DC (mDC and pDC), TAMs, CD4⁺ T cells and B cells ( Figure 2 ; Supplementary Figure S5 ). The association of CD8⁺ T cells with CXCL9/10/13 expression was the strongest. Substantially weaker correlations were detected for CD8⁺ T cells and CXCR5 or XCR1 ( Supplementary Figures S6, S7 ). Likewise, highly reproducible correlation could be observed between RCC expression of the genes of interest and CD8⁺ T cell infiltration levels as predicted by the QuanTIseq algorithm (31) ( Supplementary Figures S8, S9 ). Interestingly, levels of chemokines and the cognate receptors were not associated with Treg infiltration.

To further delineate expression of the chemokines of interest and the cognate receptors at the single-cell level, we reanalyzed a publicly available PDAC dataset (49), comprising single-cell RNAseq data from 12 RCC samples. Following cell type annotation, re-analysis revealed CXCL9 as well as CXCL10 expression mainly in myeloid cells, whereas CXCL11 was exclusively expressed in tumor cells. T cells were the major source of CXCL13, while XCL1 was predominantly expressed in NK cells. Focusing on the cognate receptors, CXCR3 was mainly expressed in T cells, CXCR5 in B cells and finally, XCR1 in myeloid cells ( Figure 3A ). In-depth analysis of CXCL9/10 and XCR1 expression in different myeloid subclusters showed that myeloid/conventional DC (cDC) were the majour source of CXCR1. On the contrary, CXCL9 and CXCL10 were predominantly expressed in M1 TAMs. In line with this finding, our FACS of tissue samples from periphery and center of RCC tumors ( Appendix Table S6 ) also identified CD14⁺ tumor-infiltrating myeloid cells as the population with the highest levels of CXCL9 protein ( Figure 3B ). In-depth analysis of CXCR3 expression in T cell subsets revealed CXCR3 expression in CD4⁺ and CD8⁺ T cells. Immunofluorescence experiments of in-house RCC samples again corroborated that CXCR3 is mainly co-expressed on CD4⁺ and CD8⁺ T cells ( Figure 3C ).

PAM clustering of the TCGA KIRC samples in respect to normalized expression levels of CXCL9/10/11, CXCR3, CXCL13, XCL1 and XCR1 yielded two distinct subsets of RCC termed the chemokine^high and chemokine^low cluster ( Supplementary Figure S2 ). This clustering structures were highly conserved in the remaining investigated collectives as evident from comparable fractions of explained clustering variance and similar distribution of chemokine^high and chemokine^low sample frequencies within the cohorts ( Supplementary Figure S3A ). The strongest differences in expression of the clustering genes between the chemokine^high and chemokine^low cluster were detected for CXCL9/10/11, CXCR3 and CXCL13 ( Figure 4 ; Supplementary Figure S3B ). In line with the correlation results, xCell-predicted infiltration of central memory CD8⁺ T cells was significantly higher in the chemokine^high than in the chemokine^low cluster. Levels of xCell-predicted TAM, activated mDC, RCC and B cell infiltration were significantly upregulated in chemokine^high cancers in the majority of analyzed cohorts ( Figure 5 ; Supplementary Figure S10 ).

Except for the TCGA KIRC cohort displaying a significant enrichment of late stage tumors in the chemokine^high cluster, no significant differences in distribution of tumor stages and MSKCC risk groups could be observed neither between the two clusters ( Supplementary Figure S11 ) nor at the level of single chemokine genes (data not shown). No differences in best therapy response between the chemokine clusters could be detected neither for everolimus nor for nivolumab in the CM collectives ( Supplementary Figure S12 ). Similarly, OS and RFS times in the two clusters did not differ significantly ( Figure 6 ; Supplementary Figure S13 ).

Finally, in multi-parameter Ridge Cox modeling of OS as a function of age, sex, grade, tumor and metastasis stage together with gene expression, the influence of chemokine levels was found to be marginal as compared with age, grade, tumor and metastasis stage ( Supplementary Figure S14A ). As evident from a comparison of concordance indexes of the multi-parameter model consisting of solely clinical and demographic variables and the model including gene expression parameters, there was no add-on value of the investigated gene expression levels to the prediction accuracy of the OS in the TCGA KIRC, E-MTAB 1980 or RECA-EU collectives ( Supplementary Figure S14B ).

To elucidate such mechanisms of activation and suppression of CD8⁺ T cell response in RCC, we delved into differences in transcriptomes between chemokine^high and chemokine^low clusters. In gene set variance analysis (34), ssGSEA scores of gene signatures of 26 Reactome pathways were found significantly upregulated in the chemokine^high cluster as compared with the chemokine^low cluster. These included PD-1 signaling involved in T cell exhaustion and signatures of pathways employing JAK/STAT signaling: IFN (IFN)α/β, IFNγ, interleukin 9 (IL9) and IL21 signaling ( Figure 7 ; Supplementary Table S1 ). Of note, the upregulation of JAK/STAT along with an activated chemokine and cytokine signaling, and increased cytotoxicity in the chemokine^high cluster were identified by the SPIA algorithm modelling activity of signaling pathways based on differential gene expression (45) ( Supplementary Figure S15 ).

To investigate the JAK/STAT signaling and T cell exhaustion process in more details, we resorted to identification of genes differentially expressed in the chemokine^high cluster as compared with chemokine^low RCCs ( Supplementary Tables S2, S3 ). Out of 140 transcripts associated with the IFNα/β and IFNγ Reactome pathways, between 21 and 71 genes were upregulated in chemokine^high RCCs of the TCGA KIRC, E-MTAB 1980, GSE73731, GSE167093, RECA-EU and CM 025 collectives. Only two IFN signaling/related genes were differentially regulated in the smallest CM 010 cohort. The top regulated IFN signaling-related genes encompassed IFNγ (IFNG), genes of signaling mediators (JAK2, STAT1, ISG15), canonical transcriptional targets of IFN signaling (IRF and GPB family), anti-viral defense genes (OAS1, OASL, SAMHB), antibody receptors (FCGR1B, TRIM21) and genes of immunoproteasome and antigen-presenting complex (PSMB8, B2M, HLA family). In particular, high levels of IFNG and direct JAK2/STAT1 transcriptional targets suggest high activity and functionality of IFNγ-mediated signaling in chemokine^high RCCs ( Supplementary Figures S16, S17 ). Of note, CXCL9/10/11 belong to classical IFNγ-stimulated genes and as such may perpetuate CD8⁺ T cell recruitment to the IFN-rich milieu of chemokine^high RCCs.

A sizable fraction of genes related to T cell exhaustion (18–21) was upregulated in the chemokine^high cluster as compared with chemokine^low cluster of the TCGA KIRC, E-MTAB 1980, GSE73731, GSE167093, RECA-EU and CM 025 collectives (total genes: 87, upregulated genes: 17 - 64). These upregulated genes included transcripts of surface co-stimulation/inhibition molecules (TIGIT, CD27, ICOS, LAG3, TNFSRF9), cytotoxic proteins (GZMB) and transcription factors driving T cell differentiation, persistence and exhaustion (TOX, EOMES, BATF) ( Supplementary Figures S18, S19 ). Genes coding for immune checkpoint proteins, PDCD1, CD274 (PD-L1), PDCD1LG2 (PD-L2), CTLA4 and HAVCR2 (TIM3), were found to be expressed at significantly higher levels in the chemokine^high than chemokine^low cluster ( Figure 8 , Supplementary Figure S20 ). In sum, the abundant CD8⁺ T cell infiltration and high activity of potentially anti-tumorigenic IFN signaling in the chemokine^high cluster is strongly counteracted by multiple redundant immune checkpoint and T cell exhaustion processes.

Finally, we sought to investigate potential alterations of cellular metabolism in the chemokine^high cluster as compared with chemokine^low cluster based on whole-genome differences in gene expression (46, 47) ( Supplementary Table S5 ). Substantial fractions of reactions involved in oxidative energy metabolism such as fatty acid oxidation, citric acid cycle (also known as TCA or Krebs cycle) and oxidative phosphorylation (OXPHOS) were found to be inhibited in chemokine^high RCCs. In turn, one-third of enzymatic reactions of tryptophan metabolism (TRP) implicated in T cell immunosuppression were significantly upregulated in chemokine^high RCCs of the TCGA KIRC collective ( Figure 9 ). In more detail, widespread significant inhibition of multiple fatty acid oxidation reactions was observed in chemokine^high cancers of the TCGA KIRC, E-MTAB 1980, GSE73731, GSE167093 and CM 025 nivolumab cohorts ( Supplementary Figure S21 ). An overall reduced activity of the citric acid cycle was particularly evident in the chemokine^high cluster of the four largest collectives (TCGA KIRC, E-MTAB 1980, GSE73731 and GSE167093) and the 2-oxoglutarate - malate section of the cycle ( Supplementary Figure S22 ). In OXPHOS, sustained inhibition of the complexes I and IV could be detected ( Supplementary Figure S23 ).

In the chemokine^high cluster of the TCGA KIRC, E-MTAB 1980, GSE73731, GSE167093 and RECA-EU cancers, significant increase (>50%) in activity of IDO was detected, a reaction responsible for TRP depletion and the first step of synthesis of immunosuppressive kynurenine and quinolinic acid. In the CM collectives, an increased predicted IDO activity was evident, yet missed statistical significance following multiple testing correction. Furthermore, we could observe an elevated activity of kynureninase (KYNU), a key enzyme in quinolinate biosynthesis, in chemokine^high RCCs as compared with the chemokine^low cluster ( Supplementary Figure S24 ). These findings were corroborated at the gene expression level in the TCGA KIRC, E-MTAB 1980, GSE73731, GSE167093, where both isoforms of IDO, IDO1 and IDO2, were significantly upregulated in the chemokine^high cluster. Significantly higher levels of IDO1 could be detected in the RECA-EU and CM 025 everolimus chemokine^high cancers as well. KYNU gene was also expressed at significantly increased levels in chemokine^high RCCs in four out of eight investigated cohorts ( Figure 10 ; Supplementary Figure S25 ).