A chemokine network of T cell exhaustion and metabolic reprogramming in renal cell carcinoma

Abstract

Renal cell carcinoma (RCC) is frequently infiltrated by immune cells, a process which is governed by chemokines. CD8⁺ T cells in the RCC tumor microenvironment (TME) may be exhausted which most likely influence therapy response and survival. The aim of this study was to evaluate chemokine-driven T cell recruitment, T cell exhaustion in the RCC TME, as well as metabolic processes leading to their functional anergy in RCC. Eight publicly available bulk RCC transcriptome collectives (n=1819) and a single cell RNAseq dataset (n=12) were analyzed. Immunodeconvolution, semi-supervised clustering, gene set variation analysis and Monte Carlo-based modeling of metabolic reaction activity were employed. Among 28 chemokine genes available, CXCL9/10/11/CXCR3, CXCL13/CXCR5 and XCL1/XCR1 mRNA expression were significantly increased in RCC compared to normal kidney tissue and also strongly associated with tumor-infiltrating effector memory and central memory CD8⁺ T cells in all investigated collectives. M1 TAMs, T cells, NK cells as well as tumor cells were identified as the major sources of these chemokines, whereas T cells, B cells and dendritic cells were found to predominantly express the cognate receptors. The cluster of RCCs characterized by high chemokine expression and high CD8⁺ T cell infiltration displayed a strong activation of IFN/JAK/STAT signaling with elevated expression of multiple T cell exhaustion-associated transcripts. Chemokine^high RCCs were characterized by metabolic reprogramming, in particular by downregulated OXPHOS and increased IDO1-mediated tryptophan degradation. None of the investigated chemokine genes was significantly associated with survival or response to immunotherapy. We propose a chemokine network that mediates CD8⁺ T cell recruitment and identify T cell exhaustion, altered energy metabolism and high IDO1 activity as key mechanisms of their suppression. Concomitant targeting of exhaustion pathways and metabolism may pose an effective approach to RCC therapy.

1 Introduction

Renal cell carcinoma (RCC) differs from other cancer entities by its high immunogenicity, which includes the efficient recruitment of tumor-infiltrating immune cells. As a result, the RCC tumor microenvironment (TME) harbors heterogenous mixtures of leukocyte subsets (1). Tumor-infiltrating lymphocytes (TILs) have previously attracted clinical research interest. In a form of adoptive T cell therapy, TILs were expanded ex vivo with high-dose IL-2 and subsequently re-infused into the patient (2). However, clinical efficacy of TIL therapy remained low. More recent studies demonstrated that CD8⁺ T cells in the RCC TME may be exhausted, express heterogenous phenotypes and express immune-evasive molecules (PD1, PD-L1, PD-L2, and CTLA4) (3, 4). In addition, the elevated CD8⁺ T cell presence corresponded to a higher frequency of BAP1 mutations (3), which is a key regulator of cancer-associated pathways (5). Exhausted T cells in the TME, which express multiple immune checkpoints, have been proposed to mediate resistance to immunotherapy with low therapeutic response rates (6–8). In RCC patients treated within the CheckMate-010 trial, a high percentage of CD8⁺ T cells expressing PD-1 but not TIM-3 and LAG-3 were positively associated with longer progression-free survival on the anti-PD-1 antibody nivolumab, suggesting that the higher predictive value of these cells might be related to their less exhausted phenotype and their ability to be more efficiently reactivated during PD-1 blockade (9). Additional immune cells such as CD4⁺CD25⁺FoxP3⁺ regulatory T cells (Tregs), myeloid-derived suppressor cells and tumor-associated macrophages (TAMs) (1) as well as dendritic cells (DCs) (10) have been detected in the RCC TME and may also affect anti-tumor immunity in different ways. Thus, the composition of TME-infiltrating immune cells and the functional phenotype of the individual infiltrating immune cell types are critical determinant factors in cancer prognosis and outcome.

The recruitment of immune cells to the TME is governed by chemokines (11), a large family of chemotactic cytokines, which are grouped into four subfamilies: CXC, CC, (X)C, and CX3C (12). Two members of the CXC family, CXCL9 (also known as MIG) and CXCL10 (also known as IP-10) are associated with Th1-type immune response by recruiting natural killer (NK) cells, CD4⁺ Th1 and CD8⁺ cytotoxic lymphocytes, which all contribute to anti-tumoral responses (13).

The migratory capacity of immune and non-immune cells is based on chemokine receptors that allow cells to migrate along chemokine gradients. In addition to being agonists of their cognate receptors, chemokines can also act as antagonists at other chemokine receptors. CXCL9, CXCL10 and CXCL11, for example are natural antagonists for CCR3 (14).

In addition to their direct effects on anti-tumor immune responses, chemokines affect angiogenesis, cancer cell proliferation, stemness and invasiveness and can thus be key determinants of disease progression (15–17). All these observations make the chemokine/chemokine receptor network an attractive target for cancer immunotherapy. Moreover, chemokines and their receptors could be useful biomarkers of response and/or survival stratification.

Herein, published transcriptome data from eight publicly available bulk RCC collectives including 1819 cancer samples and a single cell RNAseq dataset were examined. We investigated the association between chemokine expression and T cell infiltration, metabolic changes and T cell in the RCC TME.

2 Methods

Detailed description of analysis procedures is provided in Supplementary Material.

2.1 Software

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).

2.2 Data import and transformation

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).

2.3 Genes of interest

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.

2.4 Comparison of gene expression between the normal kidney and tumor samples, correlation with non-malignant cell infiltration

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).

2.5 Survival modeling

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.

2.6 Semi-supervised clustering

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).

2.7 Differences in Reactome pathways, gene expression and signaling modulation between the clusters

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).

2.8 Biochemical reaction modulation in the CXCL9 expression strata

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).

2.9 In-house flow cytometry

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).

2.10 Single-cell RNAseq analysis

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.

2.10.1 Multiplex immunofluorescence

In-House RCC samples (n = 4) after surgery were fixed in 4% paraformaldehyde, embedded in paraffin and five-micrometer sections were used for the immunofluorescence staining. Multiplex IHC was performed using Opal 6-plex Detection Kit (cat: NEL821001KT, Akoya Biosciences, Menlo Park, USA). A multiplex panel of immune markers was developed with antibodies against CXCR3 (clone EPR25373-32, cat: ab288437, dilution 1:200, Abcam, Cambridge, MA, USA), CD4 (clone EP204, cat: 104R-26, dilution 1:50, Cell Marque), CD8 (clone C8/144B, cat: M710301-2, dilution 1:200, Dako/Agilent, Santa Clara, CA, USA), CD68 (clone PG-M1, cat: M087601-2, dilution 1:250, Dako/Agilent), cytokeratin (clone AE1/AE3, cat: MA5-13156, dilution 1:500, Thermo-Fisher). The staining procedure was performed using an automated staining system (BOND-RX; Leica Biosystems, Vienna, Austria). To visualize cell nuclei, the tissue was stained with 4’,6-diamidino-2-phenylindole (spectral DAPI, Akoya Biosciences). Slides were scanned at 20x magnification using Mantra 2 Quantitative Pathology Workstation (Akoya Biosciences) and representative images from each tissue were acquired with the Mantra Snap software version 1.0.4. Image spectral deconvolution, multispectral image analysis and cell phenotyping was carried out using the InForm Tissue Analysis Software version 2.4.10 (Akoya Biosciences).

3 Results

3.1 Characteristic of the study cohorts

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).

3.2 CXCL9/10/11, CXCL13 and XCL1 govern CD8⁺ T cell recruitment in RCC

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).

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.

Figure 2

3.3 Chemokine expression is enriched in immune and tumor cells, whereas the cognate receptors are exclusively produced by immune cells

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).

Figure 3

3.4 Clustering according to chemokine expression

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).

Figure 4

Figure 5

3.5 Chemokine^high expression shows no influence on survival and therapy response

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).

Figure 6

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).

3.6 Chemokine^high cluster is associated with IFN signaling and T cell exhaustion

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).

Figure 7

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.

Figure 8

3.7 Altered energy metabolism and immunosuppressive IDO1-mediated tryptophan degradation in Chemokine^high RCC

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).

Figure 9

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).

Figure 10

4 Discussion

RCC is considered as an immunogenic tumor with frequent infiltration of immune cells. However, infiltrating anti-tumor immune cells may become dysfunctional in the TME. In addition, regulatory immune cells such as Tregs and MDS cells can infiltrate into the TME resulting in the impairment of tumor immunity (53).

Chemokines play an essential role within the TME by differentially regulating the infiltration of various immune cell subsets (16), guiding the trafficking behavior of T cells (54) and thus, influencing therapeutic outcomes in cancer patients (16). Consequently, targeting chemokine signaling pathways might be an innovative concept to improve efficacy of current cancer therapies including immunotherapy (16). We previously investigated a possible role of different chemokine receptors in RCC showing that CCR3 was predominantly expressed in RCC cells and correlated with a higher tumor grading (55). Moreover, its ligand eotaxin-1 (CCL11), possibly up-regulated as a result of tumor-associated inflammation, might be involved in the development and progression of RCC (55). In addition to CCR3, increased expression of CXCR3 and its ligands in RCC tissue compared to normal kidney tissue have already been reported (56–58). To more systematically examine the role of chemokines and their cognate receptors in RCC we performed a comprehensive analysis of published bulk RNAseq data from a total of 1819 RCC samples and identified, within the set of 28 chemokine genes available for all eight datasets, an overexpression of CXCL9/10/11/CXCR3, CXCL13/CXCR5 and XCL1/XCR1 in RCC compared with normal renal tissue. Collectively, the results of correlation, scRNAseq and semi-supervised clustering put forward CXCL9/10/11, CXCL13 and XCL1 chemokine genes along with their receptors CXCR3, CXCR5 and XCR1 as the key molecules involved in recruitment of effector memory and central memory CD8⁺ T cells to the RCC TME. Focusing on the cellular localization of these chemokines of interest and cognate receptors, M1 TAMs, T cells, NK cells and tumor cells were the major sources of chemokines. In contrast, the cognate receptors CXCR3, CXCR5 and XCR1 were predominantly expressed on T cells (CXCR3), B cells (XCR5) and denritic cells (XCR1).

RCC shows a high frequency of metabolic reprogramming (59–61) and importantly, such local and systemic metabolic alterations can affect anti-tumor immune responses (62). We could previously demonstrate that tumor-infiltrating T cells in RCC are functionally impaired due to mitochondrial and glycolytic dysfunction (48). In accordance, in our current comprehensive analysis we found a significant correlation between the chemokine^high cluster and reduced OXPHOS. The entire process from fatty acid oxidation and TCA cycle to the electron transport chain (ETC) was significantly suppressed. Importantly, Complex I and Complex IV of the ETC were most strongly inhibited in the chemokine^high RCC cluster. This phenomenon of “nutrient competition exhaustion” by reduced OXPHOS might be explained by a recent study which has shown that CD4⁺ T cells preferentially differentiate towards Tregs when OXPHOS is inhibited (63).

Among transcripts differentially regulated in the chemokine^high cluster, multiple STAT-activated genes were found, suggestive of activation of IFN/JAK/STAT signaling. This is of interest as this specific signaling pathway drives for example the expression of CXCL9/10/11 and CXCR3 (64). Moreover, the same pathway promotes IDO-1-mediated degradation of tryptophan resulting in the accumulation of the immunosuppressive metabolite kynurenine (65). We recently described that expression of IDO-1, an enzyme that catalyzes TRP and induces the accumulation of kynurenine metabolites, was predominantly expressed in tumor endothelial cells and was mostly absent from RCC tumor cells (66). Of note, in this current analysis, the TRP degrading pathway was also upregulated in the chemokine^high RCC cluster. Moreover, a series of genes related to T cell exhaustion was upregulated in the chemokine^high cluster including transcripts of surface co-stimulation/inhibition molecules, cytotoxic proteins and transcription factors driving T cell differentiation, persistence and exhaustion as well as immune checkpoint proteins such as PDCD1, CD274 (PD-L1), PDCD1LG2 (PD-L2), CTLA4 and HAVCR2 (TIM3). Thus, results of metabolic reaction activity modeling suggest profound inhibition of OXPHOS in the chemokine^high RCC cluster along with an increased activity of the immunosuppressive IDO1-mediated degradation pathway. These processes may pose another hurdle to effective anti-tumor T cell response in RCCs hallmarked by high expression of the T cell-attracting chemokines CXCL9/10/11, CXCL13 and XCL1. In sum, the T cell-inflamed TME in RCC is characterized by CD8⁺ T cell infiltration, which is mediated by these specific chemokines, and an IFNγ signature, indicating a strong interplay between tumor cells and immune cells. However, CD8⁺ T cell-inflamed RCCs may activate different immunosuppressive redundant immune checkpoint and T cell exhaustion pathways such as IDO1 and/or PD-L1/2, CTLA4, TIGIT, LAG3 and TIM3 (Figure 11). Thus, upregulation of immunosuppressive pathways within the TME is more intrinsically driven by immune cells itself rather than by tumor cells (67).

Figure 11

Focusing on the predictive and prognostic role of our identified chemokine signature, we could neither observe significant differences in survival nor in response to immunotherapy between the chemokine^high and chemokine^low RCC cluster. None of the chemokine genes used for cluster definition correlated consistently and significantly with OS or RFS. This fact may suggest, that potential beneficial effects of CD8⁺ T cells recruited via CXCL9/10/11, CXCL13 and XCL1 are likely outweighed by potent activated immunosuppressive properties of the RCC TME as described here. An additional explanation relates to the fact that chemokines may differentially regulate the biological function of T cells, resulting in distinct anti- and pro-tumoral effects in the TME (64, 68). While CXCL9/10/11 may exhibit indirect anti-tumor effects through the recruitment of CXCR3-expressing T cells, these chemokines and their cognate receptor may also contribute to tumor development and metastasis via T cell-independent mechanisms (56). Activation of CXCR3 on different tumor cell types has been shown to prevent apoptosis and promote proliferation. Moreover, immunohistochemical detection of CXCR3 on localized RCC correlated with poor disease prognosis (69). In addition to promoting and sustaining tumor development, CXCR3 may facilitate tumor cell dissemination, for instance, to the lymph nodes (64, 68). Thus, when RCC acquires CXCR3 expression and possibly also CXCL9 expression, this fact can generate an autocrine CXCL9/CXCR3 axis that supports tumor progression and metastasis. CXCL13, originally identified as a B cell chemoattractant, and its receptor CXCR5 have also emerged as key players of carcinogenesis and cancer progression (70). They can also act via autocrine and paracrine signals between the TME and the tumor cells itself (70). In RCC, CXCL13 promotes tumor cell proliferation and migration by activation of PI3K/Akt/mTOR signaling. Thus, CXCL13 up-regulation has been shown to correlate with advanced disease stage and poor prognosis (71). In addition, XCL1 positively correlated with the expression of several inhibitory receptors and cytokines that are associated with T cell exhaustion and several immune checkpoints in ccRCC, thus resulting in a significantly poorer survival (72). However, 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 chemokine gene expression levels to the prediction accuracy of OS or response to immunotherapy. Despite the fact that chemokines included in our chemokine^high cluster seem to be good candidates to predict response to immunotherapy (73, 74), RCCs with chemokine^high expression did not show better survival or response to immunotherapy in the present study. One reason we show for this is an increased activation of immunosuppressive pathways such as IDO, which is known to lead to resistance to immunotherapy (75). In RCC, we also showed in a previous work that IDO1 expression correlates positively with increased CD8⁺ T cells reflecting a T cell-inflamed TME in RCC (66). IDO1 is responsible for immune escape mechanisms for tumor cells contributing to T cell exhaustion (76). As a consequence, accelerated breakdown of TRP with increased IDO-1 activity is associated with disease progression, decreased OS and poor prognosis in different cancer entities (77). Thus, blocking more than one immunosuppressive pathway combining PD-1/PD-L1 inhibitors with IDO1 inhibitors may improve the therapeutic response to immunotherapy especially in chemokine^high RCCs. Early data already show a promising response of the combination pembrolizumab and the IDO inhibitor epacadostat in metastatic RCC (78). The corresponding phase 3 study (KEYNOTE-679/ECHO-302) is currently ongoing.

There are some limitations of our research. The analysis of published bulk RNAseq data does not allow us to link CXCL9/10/11, CXCL13 and XCL1 expression with T cell priming, differentiation and exhaustion. Analogically, it does not explain how CXCL9/10/11, CXCL13 and XCL1 overexpression induces reprogramming of energy metabolism and increased IDO-mediated tryptophan degradation. To address that and validate our bioinformatic analysis results, ex vivo assays and protein level analyses including both chemokine^high and chemokine^low RCC are considered as a next research step.

5 Conclusion

XCL9/10/11/CXCR3, CXCL13/CXCR5 and XCL1/XCR1 expression defines a chemokine^high subset of RCC characterized by abundant CD8+ T cell infiltration and active IFN/JAK/STAT signaling. These favorable prognostic features are overweighed by redundant immunosuppression processes involving T cell exhaustion, suppression of oxidative energy metabolism and OXPHOS, and IDO1-mediated tryptophan degradation. As a result, no survival or therapy response benefit could be observed for chemokine^high RCCs. Our results stress the importance of therapeutic combinations targeting multiple immunosuppressive pathways and metabolic reprogramming in RCC.

References

• 1

  MierJW. The tumor microenvironment in renal cell cancer. Curr Opin Oncol (2019) 31(3):194–9. doi: 10.1097/CCO.0000000000000512

  • CrossRef
  • Google Scholar

• 2

  RosenbergSA. IL-2: The first effective immunotherapy for human cancer. J Immunol (2014) 192(12):5451–8. doi: 10.4049/jimmunol.1490019

  • CrossRef
  • Google Scholar

• 3

  ClarkDJDhanasekaranSMPetraliaFPanJSongXHuYet al. Clinical proteomic tumor analysis consortium. Integrated proteogenomic characterization of clear cell renal cell carcinoma. Cell (2019) 179(4):964–983.e31.

  • Google Scholar

• 4

  WuXJiangDLiuHLuXLvDLiangL. CD8⁺T cell-based molecular classification with heterogeneous immunogenomic landscapes and clinical significance of clear cell renal cell carcinoma. Front Immunol (2021) 12:745945. doi: 10.3389/fimmu.2021.745945

  • CrossRef
  • Google Scholar

• 5

  TestaJRCheungMPeiJBelowJETanYSementinoEet al. Germline BAP1 mutations predispose to malignant mesothelioma. Nat Genet (2011) 43(10):1022–5. doi: 10.1038/ng.912

  • CrossRef
  • Google Scholar

• 6

  WherryEJKurachiM. Molecular and cellular insights into T cell exhaustion. Nat Rev Immunol (2015) 15(8):486–99. doi: 10.1038/nri3862

  • CrossRef
  • Google Scholar

• 7

  BlackburnSDShinHFreemanGJWherryEJ. Selective expansion of a subset of exhausted CD8 T cells by alphaPD-L1 blockade. Proc Natl Acad Sci U.S.A. (2008) 105(39):15016–21. doi: 10.1073/pnas.0801497105

  • CrossRef
  • Google Scholar

• 8

  XuWAtkinsMBMcDermottDF. Checkpoint inhibitor immunotherapy in kidney cancer. Nat Rev Urol (2020) 17(3):137–50. doi: 10.1038/s41585-020-0282-3

  • CrossRef
  • Google Scholar

• 9

  PignonJCJegedeOShuklaSABraunDAHorakCEWind-RotoloMet al. irRECIST for the evaluation of candidate biomarkers of response to nivolumab in metastatic clear cell renal cell carcinoma: Analysis of a phase II prospective clinical trial. Clin Cancer Res (2019) 25(7):2174–84. doi: 10.1158/1078-0432.CCR-18-3206

  • CrossRef
  • Google Scholar

• 10

  RadmayrCBöckGHobischAKlockerHBartschGThurnherM. Dendritic antigen-presenting cells from the peripheral blood of renal-cell-carcinoma patients. Int J Cancer (1995) 63(5):627–32. doi: 10.1002/ijc.2910630505

  • CrossRef
  • Google Scholar

• 11

  Mollica PoetaVMassaraMCapucettiABonecchiR. Chemokines and chemokine receptors: New targets for cancer immunotherapy. Front Immunol (2019) 10:379. doi: 10.3389/fimmu.2019.00379

  • CrossRef
  • Google Scholar

• 12

  ZlotnikAYoshieO. The chemokine superfamily revisited. Immunity (2012) 36(5):705–16. doi: 10.1016/j.immuni.2012.05.008

  • CrossRef
  • Google Scholar

• 13

  WangRHuangK. CCL11 increases the proportion of CD4+CD25+Foxp3+ treg cells and the production of IL-2 and TGF-β by CD4+ T cells via the STAT5 signaling pathway. Mol Med Rep (2020) 21(6):2522–32.

  • Google Scholar

• 14

  XanthouGDuchesnesCEWilliamsTJPeaseJE. CCR3 functional responses are regulated by both CXCR3 and its ligands CXCL9, CXCL10 and CXCL11. Eur J Immunol (2003) 33(8):2241–50. doi: 10.1002/eji.200323787

  • CrossRef
  • Google Scholar

• 15

  BulePAguiarSIAires-Da-SilvaFDiasJNR. Chemokine-directed tumor microenvironment modulation in cancer immunotherapy. Int J Mol Sci (2021) 22(18):9804. doi: 10.3390/ijms22189804

  • CrossRef
  • Google Scholar

• 16

  NagarshethNWichaMSZouW. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol (2017) 17(9):559–72. doi: 10.1038/nri.2017.49

  • CrossRef
  • Google Scholar

• 17

  KohliKPillarisettyVGKimTS. Key chemokines direct migration of immune cells in solid tumors. Cancer Gene Ther (2022) 29(1):10–21. doi: 10.1038/s41417-021-00303-x

  • CrossRef
  • Google Scholar

• 18

  WickhamHAverickMBryanJChangWMcGowanLFrançoisRet al. Welcome to the tidyverse. J Open Source Softw (2019) 4:1686. doi: 10.21105/joss.01686

  • CrossRef
  • Google Scholar

• 19

  KassambaraA. Rstatix: Pipe-friendly framework for basic statistical tests (2021). Available at: https://cran.r-project.org/package=rstatix.

  • Google Scholar

• 20

  CsardiGNepuszT. The igraph software package for complex network research. InterJournal (2006) 1695.

  • Google Scholar

• 21

  WickhamH. ggplot2: Elegant graphics for data analysis. 1st ed. New York: Springer-Verlag (2016). Available at: https://ggplot2.tidyverse.org.

  • Google Scholar

• 22

  KassambaraAKosinskiMBiecekP. Survminer: Drawing survival curves using ‘ggplot2’ (2016). Available at: https://cran.r-project.org/package=survminer.

  • Google Scholar

• 23

  BriatteFBojanowskiMCanouilMCharlop-PowersZFisherJCJohnsonKet al. Ggnetwork: Geometries to plot networks with ‘ggplot2’ (2021). Available at: https://cran.r-project.org/package=ggnetwork.

  • Google Scholar

• 24

  CreightonCJMorganMGunaratnePHWheelerDAGibbsRARobertsonGet al. Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature (2013) 499:43–9. doi: 10.1038/nature12222

  • CrossRef
  • Google Scholar

• 25

  BraunDAHouYBakounyZFicialMSant’ AngeloMFormanJet al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med (2020) 26:909. doi: 10.1038/s41591-020-0839-y

  • CrossRef
  • Google Scholar

• 26

  WeiXChoudhuryYLimWKAnemaJKahnoskiRJLaneBet al. Recognizing the continuous nature of expression heterogeneity and clinical outcomes in clear cell renal cell carcinoma. Sci Rep (2017) 7. doi: 10.1038/S41598-017-07191-Y

  • CrossRef
  • Google Scholar

• 27

  LaskarRSLiPEcsediSAbedi-ArdekaniBDurandGRobinotNet al. Sexual dimorphism in cancer: insights from transcriptional signatures in kidney tissue and renal cell carcinoma. Hum Mol Genet (2021) 30:343–55. doi: 10.1093/hmg/ddab031

  • CrossRef
  • Google Scholar

• 28

  RECA. RECA-EU, ICGC data portal. Available at: https://dcc.icgc.org/projects/RECA-EU (Accessed August 16, 2022).

  • Google Scholar

• 29

  SatoYYoshizatoTShiraishiYMaekawaSOkunoYKamuraTet al. Integrated molecular analysis of clear-cell renal cell carcinoma. Nat Genet (2013) 45:860–7. doi: 10.1038/ng.2699

  • CrossRef
  • Google Scholar

• 30

  SeanDMeltzerPS. GEOquery: A bridge between the gene expression omnibus (GEO) and BioConductor. Bioinf (Oxford England) (2007) 23:1846–7. doi: 10.1093/bioinformatics/btm254

  • CrossRef
  • Google Scholar

• 31

  FinotelloFMayerCPlattnerCLaschoberGRiederDHacklHet al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med (2019) 11:34. doi: 10.1186/s13073-019-0638-6

  • CrossRef
  • Google Scholar

• 32

  SturmGFinotelloFListM. Immunedeconv: An r package for unified access to computational methods for estimating immune cell fractions from bulk RNA-sequencing data. Methods Mol Biol (Clifton NJ) (2020) 2120:223–32. doi: 10.1007/978-1-0716-0327-7_16

  • CrossRef
  • Google Scholar

• 33

  AranDHuZButteAJ. xCell: Digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1

  • CrossRef
  • Google Scholar

• 34

  HänzelmannSCasteloRGuinneyJ. GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 35

  CaiMCZhaoXCaoMMaPChenMWuJet al. T-Cell exhaustion interrelates with immune cytolytic activity to shape the inflamed tumor microenvironment. J Pathol (2020) 251:147–59. doi: 10.1002/path.5435

  • CrossRef
  • Google Scholar

• 36

  WoronieckaKChongsathidkietPRhodinKKemenyHDechantCHarrison FarberSet al. T-Cell exhaustion signatures vary with tumor type and are severe in glioblastoma. Clin Cancer Res (2018) 24:4175–86. doi: 10.1158/1078-0432.CCR-17-1846

  • CrossRef
  • Google Scholar

• 37

  KocherFPucciniAUntergasserGMartowiczAZimmerKPircherAet al. Multi-omic characterization of pancreatic ductal adenocarcinoma relates CXCR4 mRNA expression levels to potential clinical targets. Clin Cancer Res (2022) 28:4957–67. doi: 10.1158/1078-0432.CCR-22-0275

  • CrossRef
  • Google Scholar

• 38

  ZhangZChenLChenHZhaoJLiKSunJet al. Pan-cancer landscape of T-cell exhaustion heterogeneity within the tumor microenvironment revealed a progressive roadmap of hierarchical dysfunction associated with prognosis and therapeutic efficacy. eBioMedicine (2022) 83:104207. doi: 10.1016/j.ebiom.2022.104207

  • CrossRef
  • Google Scholar

• 39

  BenjaminiYHochbergY. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc: Ser B (Methodological) (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar

• 40

  TherneauTMGrambschPM. Modeling survival data: Extending the cox model. 1st ed. New York: Springer Verlag (2000).

  • Google Scholar

• 41

  FriedmanJHastieTTibshiraniR. Regularization paths for generalized linear models via coordinate descent. J Stat Softw (2010) 33:1–22. doi: 10.18637/jss.v033.i01

  • CrossRef
  • Google Scholar

• 42

  SchubertERousseeuwPJ. Faster k-medoids clustering: Improving the PAM, CLARA, and CLARANS algorithms. In: Lecture notes in computer science (including subseries lecture notes in artificial intelligence and lecture notes in bioinformatics). Springer (2019). p. 171–87.

  • Google Scholar

• 43

  DrostH-G. Philentropy: Information theory and distance quantification with r. J Open Source Softw (2018) 3:765. doi: 10.21105/joss.00765

  • CrossRef
  • Google Scholar

• 44

  LangeTRothVBraunMLBuhmannJM. Stability-based validation of clustering solutions. Neural Comput (2004) 16:1299–323. doi: 10.1162/089976604773717621

  • CrossRef
  • Google Scholar

• 45

  TarcaALDraghiciSKhatriPHassanSSMittalPKimJSet al. A novel signaling pathway impact analysis. Bioinformatics (2009) 25:75–82. doi: 10.1093/bioinformatics/btn577

  • CrossRef
  • Google Scholar

• 46

  KingZALuJDrägerAMillerPFederowiczSLermanJAet al. BiGG models: A platform for integrating, standardizing and sharing genome-scale models. Nucleic Acids Res (2016) 44:D515–22. doi: 10.1093/nar/gkv1049

  • CrossRef
  • Google Scholar

• 47

  GavaiAKSupandiFHettlingHMurrellPLeunissenJAMVan BeekJHGM. Using bioconductor package BiGGR for metabolic flux estimation based on gene expression changes in brain. PloS One (2015) 10:e0119016. doi: 10.1371/journal.pone.0119016

  • CrossRef
  • Google Scholar

• 48

  SiskaPJBeckermannKEMasonFMAndrejevaGGreenplateARSendorABet al. Mitochondrial dysregulation and glycolytic insufficiency functionally impair CD8 T cells infiltrating human renal cell carcinoma. JCI Insight (2017) 2(12):e93411. doi: 10.1172/jci.insight.93411

  • CrossRef
  • Google Scholar

• 49

  LiRFerdinandJRLoudonKWBowyerGSLaidlawSMuyasFet al. Mapping single-cell transcriptomes in the intra-tumoral and associated territories of kidney cancer. Cancer Cell (2022) 40(12):1583–1599.e10. doi: 10.1016/j.ccell.2022.11.001

  • CrossRef
  • Google Scholar

• 50

  LiR. Mapping single cell transcriptomes in kidney cancer. Mendeley Data (2022) V1. doi: 10.17632/g67bkbnhhg

  • CrossRef
  • Google Scholar

• 51

  XuCLopezRMehlmanERegierJJordanMIYosefN. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models. Mol Syst Biol (2021) 17(1):e9620. doi: 10.15252/msb.20209620

  • CrossRef
  • Google Scholar

• 52

  GayosoALopezRXingGBoyeauPValiollah Pour AmiriVHongJet al. A Python library for probabilistic analysis of single-cell omics data. Nat Biotechnol (2022) 40(2):163–6. doi: 10.1038/s41587-021-01206-w

  • CrossRef
  • Google Scholar

• 53

  Díaz-MonteroCMRiniBIFinkeJH. The immunology of renal cell carcinoma. Nat Rev Nephrol (2020) 16(12):721–35. doi: 10.1038/s41581-020-0316-3

  • CrossRef
  • Google Scholar

• 54

  SteinJVNombela-ArrietaC. Chemokine control of lymphocyte trafficking: a general overview. Immunology (2005) 116(1):1–12. doi: 10.1111/j.1365-2567.2005.02183.x

  • CrossRef
  • Google Scholar

• 55

  JöhrerKZelle-RieserCPerathonerAMoserPHagerMRamonerRet al. Up-regulation of functional chemokine receptor CCR3 in human renal cell carcinoma. Clin Cancer Res (2005) 11(7):2459–65. doi: 10.1158/1078-0432.CCR-04-0405

  • CrossRef
  • Google Scholar

• 56

  OldhamKAParsonageGBhattRIWallaceDMDeshmukhNChaudhriSet al. T Lymphocyte recruitment into renal cell carcinoma tissue: a role for chemokine receptors CXCR3, CXCR6, CCR5, and CCR6. Eur Urol (2012) 61:385–94. doi: 10.1016/j.eururo.2011.10.035

  • CrossRef
  • Google Scholar

• 57

  GriffithsRWElkordEGilhamDERamaniVClarkeNSternPLet al. Frequency of regulatory T cells in renal cell carcinoma patients and investigation of correlation with survival. Cancer Immunol Immunother (2007) 56:1743–53. doi: 10.1007/s00262-007-0318-z

  • CrossRef
  • Google Scholar

• 58

  UtsumiTSuyamaTImamuraYFuseMSakamotoSNiheiNet al. The association of CXCR3 and renal cell carcinoma metastasis. J Urol (2014) 192(2):567–74. doi: 10.1016/j.juro.2014.01.100

  • CrossRef
  • Google Scholar

• 59

  ZhangYUdayakumarDCaiLHuZKapurPKhoEYet al. Addressing metabolic heterogeneity in clear cell renal cell carcinoma with quantitative Dixon MRI. JCI Insight (2017) 2(15):e94278. doi: 10.1172/jci.insight.94278

  • CrossRef
  • Google Scholar

• 60

  ReinfeldBIMaddenMZWolfMMChytilABaderJEPattersonARet al. Cell-programmed nutrient partitioning in the tumour microenvironment. Nature (2021) 593(7858):282–8. doi: 10.1038/s41586-021-03442-1

  • CrossRef
  • Google Scholar

• 61

  LindnerAKTulchinerGSeeberASiskaPJThurnherMPichlerR. Targeting strategies in the treatment of fumarate hydratase deficient renal cell carcinoma. Front Oncol (2022) 12:906014. doi: 10.3389/fonc.2022.906014

  • CrossRef
  • Google Scholar

• 62

  LercherABaazimHBergthalerA. Systemic immunometabolism: Challenges and opportunities. Immunity (2020) 53(3):496–509. doi: 10.1016/j.immuni.2020.08.012

  • CrossRef
  • Google Scholar

• 63

  ShinBBenavidesGAGengJKoralovSBHuHDarley-UsmarVMet al. Mitochondrial oxidative phosphorylation regulates the fate decision between pathogenic Th17 and regulatory T cells. Cell Rep (2020) 30(6):1898–1909.e4. doi: 10.1016/j.celrep.2020.01.022

  • CrossRef
  • Google Scholar

• 64

  ReyndersNAbboudDBaragliANomanMZRogisterBNiclouSPet al. The distinct roles of CXCR3 variants and their ligands in the tumor microenvironment. Cells (2019) 8(6):613. doi: 10.3390/cells8060613

  • CrossRef
  • Google Scholar

• 65

  CampiaIBuondonnoICastellaBRolandoBKopeckaJGazzanoEet al. An autocrine Cytokine/JAK/STAT-signaling induces kynurenine synthesis in multidrug resistant human cancer cells. PloS One (2015) 10(5):e0126159. doi: 10.1371/journal.pone.0126159

  • CrossRef
  • Google Scholar

• 66

  SeeberAKlinglmairGFritzJSteinkohlFZimmerKCAignerFet al. High IDO-1 expression in tumor endothelial cells is associated with response to immunotherapy in metastatic renal cell carcinoma. Cancer Sci (2018) 109(5):1583–91. doi: 10.1111/cas.13560

  • CrossRef
  • Google Scholar

• 67

  SprangerSSpaapenRMZhaYWilliamsJMengYHaTTet al. Up-regulation of PD-L1, IDO, and t(regs) in the melanoma tumor microenvironment is driven by CD8(+) T cells. Sci Transl Med (2013) 5(200):200ra116. doi: 10.1126/scitranslmed.3006504

  • CrossRef
  • Google Scholar

• 68

  KarinN. CXCR3 ligands in cancer and autoimmunity, chemoattraction of effector T cells, and beyond. Front Immunol (2020) 11:976. doi: 10.3389/fimmu.2020.00976

  • CrossRef
  • Google Scholar

• 69

  KlatteTSeligsonDBLeppertJTRiggsSBYuHZomorodianNet al. The chemokine receptor CXCR3 is an independent prognostic factor in patients with localized clear cell renal cell carcinoma. J Urol (2008) 179(1):61–6. doi: 10.1016/j.juro.2007.08.148

  • CrossRef
  • Google Scholar

• 70

  KazanietzMGDurandoMCookeM. CXCL13 and its receptor CXCR5 in cancer: Inflammation, immune response, and beyond. Front Endocrinol (Lausanne) (2019) 10:471. doi: 10.3389/fendo.2019.00471

  • CrossRef
  • Google Scholar

• 71

  ZhengZCaiYChenHChenZZhuDZhongQet al. CXCL13/CXCR5 axis predicts poor prognosis and promotes progression through PI3K/AKT/mTOR pathway in clear cell renal cell carcinoma. Front Oncol (2019) 8:682. doi: 10.3389/fonc.2018.00682

  • CrossRef
  • Google Scholar

• 72

  ChenZWuRMaJZhengJ. C chemokines are prognostic biomarkers correlated with diverse immune cell infiltrations in clear cell renal cell carcinoma. Transl Cancer Res (2022) 11(8):2501–22. doi: 10.21037/tcr-22-424

  • CrossRef
  • Google Scholar

• 73

  AuLHatipogluERobert de MassyMLitchfieldKBeattieGRowanAet al. Determinants of anti-PD-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell (2021) 39(11):1497–1518.e11. doi: 10.1016/j.ccell.2021.10.001

  • CrossRef
  • Google Scholar

• 74

  OzgaAJChowMTLusterAD. Chemokines and the immune response to cancer. Immunity (2021) 54(5):859–74. doi: 10.1016/j.immuni.2021.01.012

  • CrossRef
  • Google Scholar

• 75

  KocherFAmannAZimmerKGeislerSFuchsDPichlerRet al. High indoleamine-2,3-dioxygenase 1 (IDO) activity is linked to primary resistance to immunotherapy in non-small cell lung cancer (NSCLC). Transl Lung Cancer Res (2021) 10(1):304–13. doi: 10.21037/tlcr-20-380

  • CrossRef
  • Google Scholar

• 76

  MellorALBabanBChandlerPMarshallBJhaverKHansenAet al. Cutting edge: induced indoleamine 2,3 dioxygenase expression in dendritic cell subsets suppresses T cell clonal expansion. J Immunol (2003) 171(4):1652–5. doi: 10.4049/jimmunol.171.4.1652

  • CrossRef
  • Google Scholar

• 77

  BrochezLChevoletIKruseV. The rationale of indoleamine 2,3-dioxygenase inhibition for cancer therapy. Eur J Cancer (2017) 76:167–82. doi: 10.1016/j.ejca.2017.01.011

  • CrossRef
  • Google Scholar

• 78

  MitchellTCHamidOSmithDCBauerTMWasserJSOlszanskiAJet al. Epacadostat plus pembrolizumab in patients with advanced solid tumors: Phase I results from a multicenter, open-label phase I/II trial (ECHO-202/KEYNOTE-037). J Clin Oncol (2018) 36(32):3223–30. doi: 10.1200/JCO.2018.78.9602

  • CrossRef
  • Google Scholar