Skip to main content
Download PDF

Spatiotemporal evolution of the clear cell renal cell carcinoma microenvironment links intra-tumoral heterogeneity to immune escape

Genome Medicine volume 14, Article number: 143 (2022) Cite this article

Abstract

Background

Intratumoral heterogeneity (ITH) is a hallmark of clear cell renal cell carcinoma (ccRCC) that reflects the trajectory of evolution and influences clinical prognosis. Here, we seek to elucidate how ITH and tumor evolution during immune checkpoint inhibitor (ICI) treatment can lead to therapy resistance.

Methods

Here, we completed a single-arm pilot study to examine the safety and feasibility of neoadjuvant nivolumab in patients with localized RCC. Primary endpoints were safety and feasibility of neoadjuvant nivolumab. Then, we spatiotemporally profiled the genomic and immunophenotypic characteristics of 29 ccRCC patients, including pre- and post-therapy samples from 17 ICI-treated patients. Deep multi-regional whole-exome and transcriptome sequencing were performed on 29 patients at different time points before and after ICI therapy. T cell repertoire was also monitored from tissue and peripheral blood collected from a subset of patients to study T cell clonal expansion during ICI therapy.

Results

Angiogenesis, lymphocytic infiltration, and myeloid infiltration varied significantly across regions of the same patient, potentially confounding their utility as biomarkers of ICI response. Elevated ITH associated with a constellation of both genomic features (HLA LOH, CDKN2A/B loss) and microenvironmental features, including elevated myeloid expression, reduced peripheral T cell receptor (TCR) diversity, and putative neoantigen depletion. Hypothesizing that ITH may itself play a role in shaping ICI response, we derived a transcriptomic signature associated with neoantigen depletion that strongly associated with response to ICI and targeted therapy treatment in several independent clinical trial cohorts.

Conclusions

These results argue that genetic and immune heterogeneity jointly co-evolve and influence response to ICI in ccRCC. Our findings have implications for future biomarker development for ICI response across ccRCC and other solid tumors and highlight important features of tumor evolution under ICI treatment.

Trial registration

The study was registered on ClinicalTrial.gov (NCT02595918) on November 4, 2015.

Background

Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype of kidney cancer and demonstrates a high response rate to immune checkpoint inhibitors such as nivolumab, pembrolizumab, and ipilimumab [1,2,3,4]. However, only a subset of ccRCC patients respond to ICI, and biomarkers for ICI response in other disease settings such as tumor mutation burden, neoantigen load, and mismatch repair deficiency do not associate with ICI response in ccRCC [5,6,7,8]. Recently, several studies have identified transcriptomic microenvironmental features including angiogenic gene expression, T cell infiltration, and myeloid activation that correlate with response or resistance to ICI and combination therapies in ccRCC [8,9,10,11,12,13,14]. This suggests that the ccRCC microenvironment, in addition to genomic factors, influences ICI response.

In parallel, recent work has demonstrated the prevalence of ITH in untreated ccRCC [15]. This study has largely focused on heterogeneity in the presence of key driver mutations and copy number alterations and has demonstrated that ccRCC tumors follow one of a small number of evolutionary trajectories, each of which is associated with distinct patterns of genomic ITH and clinical prognosis. However, the potential for non-genomic heterogeneity in the tumor microenvironment, including but not limited to variability in the amount and identity of immune cells in spatially distinct regions of the same tumor is overlooked. Recently, we and others described substantial heterogeneity in the tumor microenvironment (TME) in several small cohorts of ccRCC tumors both in the treatment-naïve and treatment-exposed settings, raising the possibility that heterogeneity in the TME may itself shape the evolution of the tumor and its likelihood to respond to therapy [16, 17].

In this study, we hypothesized ccRCC tumors with elevated ITH constitute a genomically and immunologically distinct class of tumors, with distinguishing clonal/subclonal genomic alterations, immunologic profiles, and therapeutic response trajectories. To test this hypothesis, we utilize whole-exome sequencing (WES), whole-transcriptome sequencing (WTS), TCRseq, and histopathologic multi-regional data across a cohort of untreated and ICI-exposed patients from a phase 2 clinical trial to reveal the molecular determinants of therapy response in ccRCC (Fig. 1 and Additional file 1: Table S1). Our integrated analysis demonstrated that ITH is highly correlated among genomic, transcriptomic, and TME characteristics. ITH-high tumors are enriched for features including SETD2 and PBRM1 mutations, HLA loss of heterozygosity (HLA LOH), and CDKN2A/B loss. Immunologically, ITH-high tumors display a depletion of putative neoantigens, elevated myeloid activation, and reduced T cell diversity that are in aggregate associated with escape from the anti-tumor immune response. Premised on these observations, we developed a transcriptional signature for immune escape which correlates with distinct histopathologic patterns and is associated with ICI resistance across several diverse clinical trial cohorts.

Fig. 1
figure 1

Patient characteristics and study design. A Multi-regional multi-omics was performed on 29 patients. Serial and multi-regional sampling strategies shown schematically. B In total, 6 out of 29 patients were untreated and the rest were treated with ICI or in combination with TKI. TCRseq of PBMC was performed at 4 time points on therapy for a subset of patients. In addition, pathological review was performed to assign N-TIL (tumors sparsely infiltrated by TILs), S-TIL (tumors dominated by stromal TILs), and ES-TIL (tumors with substantial levels of both epithelial and stromal TILs) classes to a subset of patients

Full size image

Methods

Sample acquisition

After acquiring informed consent and institutional review board approval from Memorial Sloan Kettering Cancer Center (MSK), partial or radical nephrectomies were performed at MSK (New York) and stored at the MSK Translational Kidney Research Program (TKRCP). Samples were flash frozen and stored at − 80 °C prior to molecular characterization. Clinical metadata was recorded for all tumor samples. All patients represent clear cell histology and were treated via ICI alone or in combination with tyrosine kinase inhibitor (TKI). All treatments were administered prior to surgery in a neoadjuvant setting and biopsies were collected. Detailed clinical data and treatment regimen for each patient is included in Additional file 1: Table S2. After obtaining informed consent for tissue collection, samples were directly obtained from the operating room during nephrectomy. At the time of specimen extraction, samples of around 1–1.5 cm were obtained by the treating surgeon (A. A. H.) from spatially distinct tumor regions (at least 1 cm apart) and each one labeled according to its spatial location (relative to the adjacent kidney). For untreated patients (MR…) Additional file 1: Table S2, 3–5 tumor regions and an adjacent normal kidney sample (at least 2 cm away from the tumor) were obtained. For all ICB-treated patients (Nivo…SC…), three to five tumor regions, an adjacent normal kidney sample (at least 2 cm away from the tumor), and PBMC sample (in cellular preparation tubes, CPTs) were obtained.

Untreated cohort

Using and institutional database, we identified six patients with advanced or metastatic ccRCC that underwent nephrectomy with multi-regional data available, MR01,02,03,05,06, SC03. Clinical and pathologic data is available in Additional file 1: Table S2.

Neoadjuvant multi-regional cohort

This open-label, single-arm, pilot prospective study was done at Memorial Sloan Kettering Cancer Center and funded through the National Institute of Health’s Cancer Therapy Evaluation Program (CTEP) [1]. Primary endpoints were safety and feasibility of neoadjuvant Nivolumab. Patients received nivolumab (dose initially 3 mg/kg, then protocol amended to 240 mg flat dose) every 2 weeks for 4 treatments. Surgery was planned 7–14 days after the last dose. Prior to starting therapy, all patients had a kidney biopsy to confirm ccRCC, and tumor staging with renal protocol MRI and CT of the chest. After 4 doses and prior to surgery, patients also had a renal protocol MRI. Changes in primary tumor size were assessed according to Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1. Resection of the primary tumor and lymph nodes was done according to standard institutional procedures. From May 27, 2016, to September 9, 2019, 21 patients were screened and 18 were enrolled into the study of which 17 had available genomic data. Baseline patient characteristics are in Additional file 1: Table S2. All patients had localized disease at time of enrollment and biopsy-proven clear cell RCC. Perioperative and pathological details are included in Additional file 1: Table S2. Median time to nephrectomy after the last dose of nivolumab was 10.5 days (range, 9–13 days). In this small group of patients who had cancer confined to the kidney, this approach appeared safe and feasible. The details and clinical outcome of this trial are currently published in [1].

Metastatic multi-regional cohort

Using an institutional database, we identified 6 additional patients who had received ICI prior to nephrectomy (Additional file 1: Table S2). All patients had metastatic disease at time of ICI; two received anti-VEGF therapies before ICI.

Therapy response collection

All ICI-treated patients (NIVO) were only neoadjuvant treated and no responses were seen on imaging for any of the patients.

SC06: ipilimumab/nivolumab resistant (i.e., stable disease).

SC08: ipilimumab/nivolumab mixed response (partial response).

SC09: sutent then ipilimumab/nivolumab—progression of disease (progressed disease).

SC10: ipilimumab/nivolumab (complete response).

SC12: ipilimumab/nivolumab resistant (i.e., stable disease).

SC14: ipilimumab/nivolumab resistant then lenvatinib/pembrolizumab-resistant (i.e., stable disease).

Multi-regional sampling

For the prospective neoadjuvant trial and the “MR” samples, single-region biopsies were obtained preoperatively. Following nephrectomy, tumors were bivalved and 5 regions were chosen: One region from the tumor center and 4 from each quadrant (upper medial, upper later, lower medial, lower lateral). Grossly necrotic or hemorrhagic regions were avoided. For the remaining samples (those treated with definitive immunotherapy “SC”) regions were taken from distinct regions of tumors separated by 1–2 cm avoiding grossly necrotic or hemorrhagic regions).

Whole-exome sequencing

Libraries for whole-exome sequencing were generated with TruSight Oncology DNA Library Prep Kit (V1) with 40 ng input DNA per sample. TruSight Oncology index PCR products were directly used for enrichment, and target exome enrichment was performed using the IDT xGen Universal Blockers and IDT xGen Exome Research panel V1. A single-plex hybridization was done overnight at 65 °C. Accuclear dsDNA Ultra High Sensitivity assay (Biotium) was used for library quantification of the post-enriched libraries. Post enrichment libraries were normalized using bead-based normalization and pooled. Samples were sequenced with 101 bp paired-end reads on Illumina NovaSeq™ 6000 S4 flow cell using the XP workflow for individual lane loading (12-plex per lane). On average, each sample yielded 500 million reads and MEDIAN_TARGET_COVERAGE depth of 360X.

Whole-transcriptome sequencing

Libraries for whole-transcriptome RNAseq were generated with Illumina TruSeq Stranded Total RNA. One hundred nanograms RNA was used as input for Ribo-Zero rRNA Removal Kit, with Illumina TruSeq RNA UD Indexes (96 indexes) for sample indexing. Qubit dsDNA High Sensitivity assay (Thermo Fisher Scientific) was used for library quantification. Sequencing was done on Illumina NovaSeq™ 6000 S2 (36-plex) or S4 (72-plex) flow cell with 76 bp paired-end sequencing to produce ~ 200 million paired reads per library.

T cell repertoire sequencing

Libraries for T cell repertoire sequencing were generated with AmpliSeq for Illumina Library PLUS paired with AmpliSeq cDNA Synthesis for Illumina with 100 ng RNA input per cDNA synthesis reaction. The TCR beta-SR Panel was used for generating amplicons, and AmpliSeq CD Indexes Set A for Illumina were used for sample barcodes. Qubit dsDNA High Sensitivity assay (Thermo Fisher Scientific) was used for library quantification. Sequencing was done on the NextSeq 550 (41-plex) with 151 bp paired-end sequencing to produce ~ 5 million paired reads per library.

WTS pipeline

WTS raw read sequences were aligned against human genome assembly hg19 by STAR 2-pass alignment [18]. QC metrics, for example general sequencing statistics, gene feature, and body coverage, were then calculated based on the alignment result through RSeQC. WTS gene-level count values were computed by using the R package GenomicAlignments [19] over aligned reads with UCSC KnownGene [20] in hg19 as the base gene model. The union counting mode was used and only mapped paired reads after alignment quality filtering were considered. Finally, gene-level FPKM (fragments per kilobase million) and raw read count values were computed by the R package DESeq2 [21].

ESTIMATE

The ESTIMATEScore, which is the estimate of the presence of stromal and immune cells in tumor tissue, is calculated through the ESTIMATE R package [22] based on a given gene expression profile in FPKM.

Immune deconvolution analysis

Two distinct popular computational methods, ssGSEA [23] and CIBERSORT [24], were chosen for immune deconvolution analysis. Signature gene lists of immune cell types for ssGSEA were obtained from Bindea et al. [25] and Senbabaoglu et al. [4]. ssGSEA takes the sample FPKM WTS expression values as the input and computes an enrichment score for the given gene list of immune cell type relative to all other genes in the transcriptome. On the other hand, CIBERSORT also takes FPKM WTS expression values as the input but uses a signature gene expression matrix of interest immune cell types instead to compute the infiltration level of each immune cell type. The LM22 immune cell signature which was validated and published along with CIBERSORT is used. We also used FRICTION [26] to deconvolute WTS into absolute CD8 and CD4 T cells as well as CD19 B cells.

HERV quantification

We used WTS to quantify HERVs as described before [26]. Briefly, all WTS reads were aligned (using STAR aligner with optimized multi-mapping options) to a custom genome build where human reference (hg19) and HERV-specific reference are combined. Then reads aligned to non-HERV genes are removed and the rest are annotated. Three samples contained super high median HERVs (Grubbs test P < 0.05) and removed for better visualization.

WES analysis pipeline

Raw sequencing data were aligned to the hg19 genome build using the Burrows-Wheeler Aligner (BWA) version 0.7.17 [27]. Further, indel realignment, base-quality score recalibration, and duplicate-read removal were performed using the Genome Analysis Toolkit (GATK) version 3.8 [28] following raw read alignment guidelines [29]. VarScan 2 [30], Strelka v2.9.10 [31], Platypus 0.8.1 [32], Mutect2—part of GATK 4.1.4.1 [29], Somatic Sniper version 1.0.5.0 (SNVs only), and [33] were used for small variant calling, and combination of 2 out 5 callers is reported as per Cancer Genome Atlas Research Network recommendations [34]. Variants were filtered using the following criteria:

  1. (1)

    Tcov (tumor coverage) > 10 and Taf (tumor allele frequency) ≥ 0.04 and Ncov (normal coverage) > 7 and Naf (normal allele frequency) ≤ 0.01 and Tac > 4 are set to PASS

  2. (2)

    Common SNPs are eliminated by comparison to snp142.vcf

  3. (3)

    Rare variants found in dbSNP are kept if Naf = 0

  4. (4)

    Variants with Tcov < 20 or Tac < 4 are marked as low_confidence

  5. (5)

    Only variants called by more than 1 caller are reported.

  6. (6)

    Common variables gnomAD v 2.1.1 are excluded.

Variants were annotated using Ensembl Variant Effect Predictor (VEP) [35]. Additional optimization and filtering are applied for INDELS. INDELS in blacklisted regions (https://www.encodeproject.org/annotations/ENCSR636HFF/) and low mappability regions (such as repeat maskers) are excluded as per [36]. Combination of filtered SNV and INDELS are used by maftools R package to generate oncoplots and summary plots, as per author’s recommendations (https://www.bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html).

All nonsynonymous point mutations identified as above were translated into strings of 17 amino acids with the mutant amino acid situated centrally using a bioinformatics tool called NAseek. A sliding window method is used to identify the 8–11 amino acid substrings within the mutant 17-mer that had a predicted MHC Class I binding affinity of ≤ 2%Rank to one (or more) of the patient-specific HLA alleles. Binding affinity for the mutant and corresponding wild type nonamer is analyzed using NetMHCpan4.0 software. Only neoantigens with TPM > 1 are considered to be expressed.

Allele-specific copy number analysis is done by the FACETS v.6.1 [37]. Allele-specific HLA loss is determined using LOHHLA as described before [38].

RNA and TCR ITH scores

Gene- and patient-wise intra-patient heterogeneity scores were calculated using multi-region data. Data was first median-centered to remove any gene-level bias. For each gene, the difference between each pair of samples from the same tumor was calculated. The median difference between the paired-differences was then taken, yielding a gene-specific, patient-specific measure of heterogeneity. This was repeated for all genes, across all tumors, generating a matrix of gene by patient values. Gene intratumor heterogeneity values are summarized as the median value per gene across all tumors in the cohort. Patient intratumor heterogeneity values are summarized as the median value per tumor across all genes. Patient intratumor heterogeneity values represent the expected value of the absolute log2-fold change for a randomly chosen gene within a given tumor.

TCR ITH score is defined as 1 − percentage of shared clonotypes across multiple regions of tumor based on WTS. T cell clones are estimated using MiXCR application on Illumina BaseSpace (http://basespace.illumina.com/apps/). Furthermore, all ITH scores are classified as high versus low using the median as threshold.

Distinction between dedicated TCRseq and TCR clones inferred from RNAseq using MiXCR

All TCR-associated data analysis in this study (including tissue or PBMC) are based on ultra-deep T cell repertoire sequencing (targeted TCRseq) to mitigate undersampling of T cell clones except TCR ITH analysis in Fig. 3B where ITH associated with multi-regional sequencing is derived from MiXCR T cell estimates from RNAseq data due to the lack of multi-regional TCRseq for all patients.

ccRCC evolutionary subtypes and intratumor DNA heterogeneity score

DNA ITH score is calculated as the ratio of subclonal to clonal driver genomic alterations including SNVs, INDELs, and SCNA [15]. A genomic alteration is defined to be subclonal if it is present in less than half of the regions collected in each patient. Patients with enough DNA biopsies collected are classified into 1 of the 7 ccRCC evolutionary subtypes as described before [15]. We used neighbor joining tree construction in ape package in R [39] for reconstruction of tumor clones. TCGA ITH score was obtained from a previous study as measured by the number of clones estimated per sample using PhyloWGS [40]. Briefly, PhyloWGS is a method to infer tumor evolution evolutionary using the relationships between tumor subpopulations based on variant allele frequencies while considering copy number alterations.

HLA and TCR diversity

Shannon entropy is calculated to define TCR diversity [41]. We used MiXCR application on Illumina BaseSpace (http://basespace.illumina.com/apps/) for alignment and T cell clonotype identification. Immunarch (https://immunarch.com/) [42] was used for downstream analysis including visualization and data analysis. Morisita index [43] was used to measure clonotype overlap. HLA diversity index is measured as adopted from [44] as described in [26].

Neoantigen depletion

The fraction of neoantigens depleted is defined for each sample where pre-treatment data was available. We first calculated the neoantigen depletion as the number of neoantigens that were undetectable after therapy but were detected pre-treatment. The fraction of neoantigens depleted was then defined as the ratio of the total number of depleted neoantigens over total pre-treatment neoantigens. To distinguish neoantigen depletion due to contraction (immune elimination) from evasion, we exclude any neoantigens that were depleted without the presence of HLA LOH (defects in antigen presentation machinery), or reduced expression, i.e., log2(FC) <  − 1 where FC is the fold change defined as the ratio of posttreatment TPM over pre-treatment TPM after correction for tumor purity. Conversely, a neoantigen is annotated and was deleted due to immune elimination if log2(FC) >  = 0, and no HLA LOH was detected. Likewise, HERV editing is defined as the median change in the expression of immunogenic HERVs compared to pre-treatment expression. Immunogenic HERVs refer to HERV loci whose expression strongly correlates with TIL abundance, FDR < 0.05.

Weighted gene co-expression network analysis (WGCNA) and gene signature extraction

We performed WGCNA [44] on all samples where the fraction of neoantigens depleted was available similar to previously described [11]. Briefly, genes with low expression values and invariant genes, that is, genes that were expressed in < 5% of samples or had s.d. ≤ 1 for expression (log2 TPM), were filtered together with non-coding genes. The soft power of 6 was chosen based on goodness of fit to a scale-free network. We first annotate modules as JAVELIN or angiogenesis according to the Spearman correlation between the module eigengene and JAVELIN or angiogenesis ssGSEA scores (highest correlation is classified as JAVELIN or angiogenesis module). Likewise, among all modules, the module with the highest Spearman correlation with the fraction of neoantigens depleted was annotated as immune escape module (85 genes). This 85 genes’ gene signature was strongly associated with PFS of Avelumab plus Axitinib in JAVELIN Renal 101 (HR = 1.45, P = 0.02, Additional file 2: Fig. S1 A). To further refine this gene signature, we first sorted genes based on their pairwise Spearman correlation (Additional file 2: Fig. S1 B) and then selected genes with the highest Spearman correlation such that no genes have a Spearman correlation < 0.6 (Additional file 2: Fig. S1 C). This reduced the number of genes to total of 12 highly correlated genes known as immune escape signature (TIMP1, PXDN, COL15A1, OLFML2B, COL5A2, DLX5, SOX11, KLHDC8A, UNC5A, ADAMTS14, MMP11, FN1). Several genes (ADAMTS14, MMP11, FN1, COL5A1, COL5A2, and TIMP1) in this signature have previously been described as TGF-β-associated extracellular matrix genes that are linked to immune evasion and immunotherapy failure [45].

TCGA validation of CDKN2A/B loss with ITH and myeloid enrichment

We downloaded WES data from TCGA KIRC cohort from https://www.cancer.gov/tcga and processed the raw data using the same pipeline used to process data generated for the multi-regional cohort as described earlier. CDKN2A/B loss calls extracted and the association with ITH assessed. TCGA ITH score was obtained from a previous study as measured by the number of clones estimated per sample using PhyloWGS [40]. Likewise, to evaluate the association between CDKN2A/B loss and myeloid score, gene expression count data corresponding to TCGA KIRC cohort was downloaded from https://www.cancer.gov/tcga and ssGSEA of myeloid signature calculated as described earlier.

scRNAseq data analysis pipeline

scRNAseq data was downloaded from https://www.ncbi.nlm.nih.gov/sra/PRJNA705464. This repository contains a Seurat object containing the raw counts, normalized counts, dimensionality reduction post-batch correction, and cell type identities based on Krishna et al. study [16], which was used for further analysis in this study.

Statistical analysis

All statistical tests were performed in R. To calculate correlations, cor.test with Spearman’s method was used. Tests comparing distributions were performed using wilcox.test. All statistical analyses were two-sided, and p-values were Benjamini–Hochberg corrected.

Results

The landscape of microenvironmental ITH in ccRCC

To study ITH in ccRCC, we completed ultra-deep (median coverage of 360X) multi-regional whole-exome sequencing and whole-transcriptome sequencing across 142 tumor regions from 29 patients, including 6 untreated and 23 post ICI (see “Methods” and Additional file 1: Table S2. The details and clinical outcome of this trial are currently published in [1]). Tumor biopsies were extracted from different regions of the same primary tumor unless specified (Fig. 1A, B, Additional file 1: Table S2). While intratumoral genetic heterogeneity in ccRCC is well-described [46], comparatively little is known about the extent of microenvironmental heterogeneity and its relationship to other molecular features of the tumor. To measure the extent of intratumoral microenvironmental heterogeneity, we leveraged multi-regional WTS of up to 5 regions from 29 patients. Using single-sample gene set enrichment analysis (ssGSEA) of established gene signatures, we quantified the expression of several TME gene expression signatures recently proposed as biomarkers of response to ICIs and antiangiogenic agents [47] (myeloid signature [9], JAVELIN signature [11], and angiogenesis signature, see “Methods” and Additional file 1: Table S3). We confirmed that these RNA signatures accurately quantified the abundance of key immune populations using matched immunofluorescence data, including statistically significant associations between CD31/angiogenesis (p = 0.0003), CD8/JAVELIN T cell signature (p = 0.02), and CD68/Myeloid infiltration (p = 0.0013) (Additional file 2: Fig. S2).

For each RNA signature, we normalized scores to capture the magnitude of expression relative to all other profiled regions in our cohort. We then investigated the variability of microenvironmental RNA signatures across regions, finding that they demonstrated extensive heterogeneity across tumor regions from the same patient (Fig. 2A). While a small number of patients showed relatively uniform immune infiltration (e.g., NIVO02, Fig. 2A), the significantly more common phenomenon was for patients to exhibit regions both above and below the median score for a microenvironmental feature of interest (e.g., angiogenesis in MR03, JAVELIN/T-effector signatures in NIVO22).

Fig. 2
figure 2

TME ITH in ccRCC. A Intratumoral heterogeneity of several gene expression signatures across multiple tumor regions. For each patient, each tumor region is assigned a unique color. For each vertex on the radar plot, we plot the Z-score of the relevant RNA signature relative to all other samples in the cohort. Min and max radius for each feature in each panel represent min and max of that feature across the cohort. A wide (narrow) spread of dots for a given feature in a single patient corresponds to qualitatively large (small) heterogeneity of that microenvironmental RNA signature in that patient. B For each gene signature, the number of patients who were classified as high or low or a mixture of high and low across tumor regions are shown. Two patients (MR05 and NIVO10) were excluded since WTS data of only one region was available. Also, pre-treatment regions of ICI-treated patients were excluded to avoid treatment-related effects in these signatures. C Intratumoral heterogeneity of myeloid score observed across multiple regions of tumors of patients in this study

Full size image

ITH in expression signatures had substantial consequences on the accuracy of stratifying patients into high/low expressing groups based on a single-region biopsy. Using the myeloid signature (which has previously been associated with poor response to ICI [9]), we calculated the median myeloid score across all samples in the cohort and then assigned each tumor region to either a high myeloid or low myeloid group. Consequently, we observed that the majority of patients had tumor regions classified into both myeloid-high and myeloid-low regions (Fig. 2C). Given that several of these signatures are under active investigation as biomarkers of response to ICI, we investigated more generally how classification of regions into high/low was affected by ITH. Remarkably, in more than half of the patients, clinically relevant signatures (Angiogenesis, T-effector, Myeloid, and JAVELIN) could not be consistently classified as high or low (Fig. 2B, 2 patients, i.e., MR05 and NIVO10 were excluded since WTS data of only one region was available).

We hypothesized elevated microenvironmental heterogeneity may reflect the presence of underlying genomic driver alterations. To test this, we leveraged multi-regional WES data collected for these patients. Frequencies of established ccRCC driver alterations were in agreement with a previous multi-regional study by TRACERx Renal [15] (Fig. 3A). We performed unsupervised hierarchical clustering of major ccRCC driver mutations (with mutation frequency > 10% across TCGA KIRC cohort and TRACERx Renal [15]) that were previously shown to be associated with ITH [15] (i.e., VHL, PBRM1, SETD2, BAP1) and genomic alterations enriched with metastatic disease and ICI response (HLA LOH and CDKN2A/B copy number loss) [38, 44, 48], ultimately identifying two clusters (Fig. 3B, Additional file 1: Table S4). We compared the results of these clusters to aggregate, univariate measures of intratumoral DNA, RNA, and T cell receptor (TCR) heterogeneity. Interestingly, one cluster was characterized both by an enrichment of specific genomic alterations (SETD2 mutations, Fisher exact test P = 0.002; CDKN2A/B copy number loss, Fisher exact test P = 0.0001; HLA LOH, Fisher exact test P = 0.0007). This same cluster of patients, which we refer to herein as “ITH-high”, had comparable levels of tumor purity to the other “ITH-low” cluster, but demonstrated elevated ITH at the level of somatic DNA alterations, RNA, and TCR (combined Fisher exact test P = 0.0495). Moreover, by classifying patients into previously described ccRCC evolutionary subtypes (Additional file 2: Fig. S3), we observed that PBRM1-driven tumors were enriched in the ITH-high cluster (on sample level, Fisher exact test P = 0.0018), in agreement with TRACERx Renal [15]. We did not find any association between ITH and other gene mutations (Fisher exact test P > 0.5). However, this finding must be treated with caution due to our relatively small cohort size as well as low number of regions collected in some patients. These findings were robust to the number of regions collected per tumor, and we found no significant association between ITH and exposure to ICI (Fisher exact test P = 0.65, Fig. 3B); however, due to the small size of our untreated cohort, this analysis might be underpowered. Together, our results demonstrate that (1) ITH is not restricted to genomic events, but rather is pervasive in the transcriptome, microenvironment, and immune compartment of ccRCC tumors, and (2) correlates with specific somatic events at the level of individual patients (i.e., PBRM1 and SETD2 mutations, HLA LOH and CDKN2A/B loss).

Fig. 3
figure 3

Landscape of ITH in ccRCC. A Oncoprint of key ccRCC driver mutations and copy number alterations for all regions of all 29 patients in this cohort. Margin shows comparison between mutation frequency observed in this cohort and TRACERx Renal. B We have performed unsupervised hierarchical clustering of genomic features including patient level presence or absence of a small variant in VHL, PBRM1, SETD2, BAP1 (most commonly mutated genes) as well as loss of heterozygosity in HLA genes as well as 9p (which includes CDKN2A/B) SCNA which are known to affect ICI response. Heatmap shows ITH high vs low classification across data type. Annotation illustrates evolutionary subtypes and treatment status of patients. A patient is annotated as wildtype if all regions are wild type for that alteration. Cases where ITH score could not be calculated due to lack of sufficient number of biopsies are shown in gray pixels. CIN: chromosome instability. C Association between antigen presentation machinery (APM), effector T cell (Teff) and myeloid gene signatures, and ITH. Wilcox P, false discovery rate (FDR), and linear mixed effect (LME) P shown. D Intratumoral heterogeneity and myeloid score are associated with CDKN2A/B loss in TCGA KIRC cohort. E ITH low patients show a significantly higher TCR diversity, richness, and clone count

Full size image

ITH-high ccRCC tumors are immunologically distinct

Comparing the TME characteristics of ITH-high and ITH-low patients, we observed that ITH-high tumors (defined as all regions belonging to a patient who is classified as ITH-high) were characterized by high myeloid and low T cell effector (Teff) signatures (Fig. 3C). Similarly, a signature associated with antigen presentation (APM) [4] was downregulated in ITH-high patients, consistent with elevated levels of HLA LOH in the ITH-high subtype. To validate if genomic features uniquely characterizing ITH-high tumors (HLA LOH and CDKN2A/B loss) might be more generally associated with myeloid infiltration in a large, independent cohort, we obtained DNA and RNA sequencing data from the TCGA KIRC study and scored samples by the presence of CDKN2A/B loss, ITH (as measured by the number of clones estimated per sample using PhyloWGS, see “Methods”), and myeloid infiltration. This analysis confirmed that in ccRCC, CDKN2A/B loss was associated with higher levels of ITH (P = 3 × 10−5) and higher myeloid infiltration (P = 7 × 10−5) (Fig. 3D). However, the association between genomic ITH and myeloid infiltration did not reach statistical significance in TCGA KIRC cohort suggesting the association between myeloid infiltration and ITH is likely indirect through certain genomic events such as CDKN2A/B loss.

The findings above suggested that ITH-high tumors may be distinct in their immunophenotype, including in the diversity of their T cell repertoire. We therefore investigated the association between ITH and T cell diversity both peripherally and within the tumor. To do so, we compared the overlap between tissue-resident and peripheral T cells. Repertoire overlap analysis (Additional file 2: Fig. S4) illustrated a high degree of shared clonotypes across different tumor regions from the same patient, but a lack of shared clonotypes across patients. ITH-high patients demonstrated a significantly lower peripheral TCR diversity, richness, and clone count compared to ITH-low patients (Fig. 3E), suggesting that elevated heterogeneity in the primary tumor is associated with reduced peripheral immunologic diversity in a manner that is consistent with reports in other diseases [49]. However, the association between ITH and TCR diversity remains correlative and future mechanistic studies are required to establish a causal relation between these two features of tumor and immune-phenotype. Together, the above data argue that elevated molecular heterogeneity in ccRCC tumors is associated with a distinct microenvironmental and immunologic phenotype.

ICI therapy is associated with loss of putative neoantigens and HLA LOH

The clinical management of ccRCC (for which pre-surgical biopsies are often not indicated or used) makes serial profiling of primary tumors on therapy challenging, rendering our understanding of how ICI may remodel tumor physiology incomplete. To overcome this challenge, we took advantage of 16 patients from our neoadjuvant nivolumab clinical trial who had WES performed on their pre-treatment biopsies. This offered a unique opportunity to interrogate both genomic adaptations (including both somatic mutations and the expression of potentially immunogenic endogenous retroviral elements, HERVs) to ICI therapy and immunologic changes in the T cell repertoire.

Focusing first on genetic alterations, we anticipated that ICI administration would lead to elimination of some tumor clones and therefore a contraction in total mutation count. However, we observed no consistent trend in the change of either SNV or indel mutational count following ICI therapy (Additional file 2: Fig. S5). Nevertheless, the number of nonsynonymous SNVs that were predicted to bind to MHC complex in silico was consistently reduced across all patients and all biopsies except for NIVO03 (Additional file 2: Fig. S5 and Fig. 4A). An opposite trend was observed in the number of putative non-binders, suggesting a selection in favor of non-neoantigenic mutations by tumor during clonal evolution (Fig. 4A).

Fig. 4
figure 4

The landscape of heterogeneity of neoantigen depletion. A Change in the number of nonsynonymous binder SNVs (predicted in silico) and nonsynonymous non-binder SNVs compared to pre-treatment. Reduction in only putative neoantigens illustrates selective pressure and immunoediting. One sample Wilcox test P (compared to zero) is shown. B Clonality of neoantigen depletion. Only strong binders are shown. C Association between neoantigen depletion and ITH subtypes. D, E Immunoediting in an HLA-intact patient NIVO20 through reduced neoantigen expression. NKA/NKB/NKC (shown in RED) are normal adjacent tissues 1, 2, and 4 cm away from the center of the tumor; BX (shown in blue) represents pre-treatment biopsy; RA/RB/RC/RD/RE (shown in green) illustrate 5 tumor regions from the treated tumor sample. F Immunoediting with amino acid resolution. Higher phenylalanine (F) depletion compared to glutamic acid (E) and arginine (R) suggests immune selection. G Association between putative neoantigen depletion and myeloid activation across all regions of patients where pre-treatment WES data was available (n = 16 patients). H Association between the fraction of expressed putative neoantigens depleted and immune signatures. In H, correlations are calculated across different regions of the same patient, for all patients with > 3 treated, as well as pre-treatment RNA samples were available (n = 7 patients). I, J HERVs are enriched in tumors compared to normal samples and are associated with treatment. K HERV depletion association with myeloid signature. L Clonotype tracking of dominant untreated T cell clones in treated regions of patient NIVO20. The color of each ribbon shows different T cell clones, and the width is scaled corresponding to the frequency of that clone. Tissue data consists of 5 tumor regions after treatment (RA/RB/RC/RD/RE), one single normal adjacent (NKC), and one tumor region pre-treatment (BX). Likewise, PBMC data points on treatment are NIVO20-68, -54, -40, -12. M, N TCR diversity is negatively associated with neoantigen depletion and HLA LOH

Full size image

In order to characterize the clonality of putative neoantigen depletion across distinct tumor regions, we counted all 8–11 amino-acid-long putative neoantigens seen prior to treatment but deleted in at least one biopsy after treatment. Among 7 patients with at least 4 tumor regions sequenced, we observed an enrichment for putative neoantigen depletion across 4 or more sites (Fig. 4B). To assess the link between ITH and neoantigen depletion, we examined the magnitude of putative neoantigen depletion in each patient by measuring the average number of putative neoantigens deleted per biopsy (i.e., the ratio of the deleted neoantigens in a treated region compared to pre-treatment over the total number of pre-treatment neoantigens). Using this metric, we observed a strong association between neoantigen depletion and ITH subtypes (Fig. 4C). Focusing on patient NIVO20, all 6 identified depleted putative neoantigens were deleted in at least 4 regions, suggesting putative neoantigen depletion is a clonal event (Fig. 4D). Genes expressing these depleted neoantigens demonstrated a 2–threefold reduction in expression related to pre-treatment biopsy (NIVO20-RA/RB/RC/RD/RE vs NIVO20-BX) (Fig. 4E). Together with the data above, these observations suggest that ICI therapy in ccRCC is associated with the clonal loss of mutations with elevated immunogenicity.

Premised on prior reports [50] of the increased immunogenicity of hydrophobic residues, we sought to determine whether a selective pressure exists on certain neoantigens. We compared the number of amino acids preserved versus depleted upon immunotherapy and noticed a strong selection against phenylalanine (F, extremely hydrophobic) in favor of arginine (R, extremely hydrophilic) and glutamic acid (E, extremely hydrophilic) in our cohort (Fig. 4F).

Next to elucidate the link between the TME and neoantigen depletion (Additional file 1: Table S5), we compared different TME gene expression signatures and the fraction of neoantigen depleted. We observed that the fraction of neoantigens depleted was strongly associated with myeloid-high regions (n = 16 patients whose pre-ICI treatment WES data was available, Fig. 4G). The association between myeloid activation and neoantigen depletion remained strong when total number of neoantigens depleted was used (instead of fraction) (Additional file 2: Fig. S6B) or when putative neoantigen (transcriptional) expression was taken into account (n = 7 patients whose pre-treatment WTS data was available, Fig. 4H) and was not affected by variation in tumor purity (Additional file 2: Fig. S6). Furthermore, the correlation between the degree of neoantigen depletion and myeloid infiltration was also evident when examining different regions of individual patients, where highly depleted regions were associated with the highest myeloid and lowest ImmuneScore (Fig. 4H).

A recent study [51] identified tumor infiltrating lymphocyte-specific HERV epitopes that are translated, can bind to MHC I complex, and induce high-avidity cytotoxic T cells. In [51] as well as other previous reports [52], overexpression of HERVs on tumor cells has been reported and a link to ICI response has been documented [53]. To interrogate other tumor-intrinsic features associated with immune response in our cohort, we utilized our deep RNA sequencing (~ 200 million reads/library) to quantify HERV expression. HERVs were overexpressed in tumors compared to normal tissues in our cohort (Fig. 4I), and median HERV (median of all HERV loci investigated) was correlated to angiogenic expression (Additional file 2: Fig. S8A). Notably, PBRM1 mutations, which lead to further HIF upregulation [54] and angiogenic expression [55, 56], were also positively associated with HERV (Additional file 2: Fig. S8B), consistent with a recent report [57]. In agreement with [53], we then confirmed the association between the median expression of different HERV loci and TIL abundance (Additional file 2: Fig. S8A). Median HERV was anti-correlated with tumor purity; however, the association between HERV expression and TIL abundance remained valid even when HERV expression was corrected for tumor purity (Additional file 2: Fig. S8A). Conversely, we observed a significant reduction in HERV expression an observation akin to reduction in neoantigens (Fig. 4J). Likewise, we observed a strong correlation between HERV editing (i.e., change in the expression of immunogenic HERV loci after treatment, see “Methods”) and myeloid signature further highlighting the association between neoantigen depletion and myeloid enrichment (Fig. 4K). Due to the limitations of HERV quantification using WTS, we could not rule out that a strong correlation between HERV and TIL abundance might be due to expression of HERV on immune cells. However, the expression of HERV on ccRCC tumor cells has been previously shown [58] and their immunogenicity is well-established [51]. Nevertheless, rigorous determination in future studies of cell-specific expression of HERVs will be critical to understanding their putative association with ICI response.

Finally, using TCRseq of tissue-resident and peripheral T cells, we investigated the impact of ICI and neoantigen depletion on T cell diversity. Focusing again on patient NIVO20 where TCR data of multiple regions of pre-treatment and ICI-treated tumor were available, we evaluated the degree of overlap between T cell clonotypes at different regions and time points, i.e., pre-treatment, on-therapy, and post ICI treatment (Fig. 4L). Tracking dominant tissue-resident T cell clonotypes, we noticed a substantial depletion of dominant T cell clones upon ICI therapy (Fig. 4M). This observation was mirrored across our entire cohort, where we observed a strong negative association between peripheral TCR diversity and neoantigen depletion and allele-specific HLA loss across the entire cohort where PBMC TCRseq data was collected (Fig. 4M, N). Together, if validated using future mechanistic experiments, our findings suggest that neoantigen depletion in primary ccRCC tumors is associated with peripheral loss of neoantigen reactive T cells. However, at this point, no causal relationship between neoantigen loss and TCR diversity can be drawn.

Subclonal evolution underlies immune escape

In order to understand the immunologic mechanisms driving subclonal evolution after ICI, we investigated in detail patients whose tumors underwent subclonal immunoediting in distinct regions. Strikingly, subclonal reconstruction revealed recurrent subclonal evolution of HLA LOH and CDKN2A/B loss following ICI therapy (Fig. 5A). Notably, we observed HLA LOH and CDKN2A/B loss co-occur in 9 patients (Fisher exact test P = 0.003) and most tumor regions (Fisher exact P = 5 × 10−7) (Fig. 5B). Strikingly, comparing the untreated and treated regions, we only observed a significant immunological response (as measured by Th1 response) in regions without CDKN2A/B loss or HLA LOH (Fig. 5C), suggesting that HLALOH or CDKN2A/B loss are subclonal determinants of response to ICI [38, 44, 48]. This is consistent with recently published data [49] indicating the loss of 9p21—encompassing CDKN2A/B—confers a cold tumor immune microenvironment and resistance to ICI. In that study, Han et al. [49] linked 9p21 loss to a decreased abundance of B, T, CD8 T, NK cells and cytotoxic lymphocytes, and an increased fractions of macrophages, as well as reduced TCR CDR3 repertoire abundance and diversity. We interpret our observations to mean that immunoediting occurs under selective pressure by which certain tumor subclones transform to a less immunogenic phenotype through HLA LOH and CDKN2A/B loss, and that this subclonal selection can produce a highly heterogenous TME.

Fig. 5
figure 5

Branch evolution demonstrates immune evasion. A Evolutionary tree illustrates tumors can exploit concurrent HLA LOH and CDKN2A/B loss to escape immune surveillance. B Co-occurrence of HLA LOH and CDKN2A/B can be seen both across regions and patients. C Differential immune response to ICI therapy in patients with CDKN2A/B loss or HLA LOH or belonging to ITH high subtype

Full size image

To further shed light on the how tumor evolution can transform TME, we sought to analyze the spatial distribution of TILs within the TME and their interaction with the stromal compartment using immunohistochemical data. Following A.W. Zhang and colleagues [59], a dedicated genitourinary pathologist classified tumor regions into 3 subtypes according to the co-localization of tumor infiltrating lymphocytes and tumor cells based on lymphocyte morphology: N-TIL (tumors sparsely infiltrated by TILs), S-TIL (tumors dominated by stromal TILs), and ES-TIL (tumors with substantial levels of both epithelial and stromal TILs) (Additional file 2: Fig. S9, Additional file 1: Table S6). We observed that an ES-TIL enriched TME is strongly associated with regions with HLA LOH (ES = 4, N = 7, S = 5 compared to ES = 2, N = 32, S = 21 in HLA-intact regions, Fisher’s exact test P = 0.036) or loss of CDKN2A/B (ES = 4, N = 7, S = 9 compared to ES = 2, N = 32, S = 17 in regions without loss of CDKN2A/B, Fisher’s exact test P = 0.03) whereas N-TIL pathology is linked with regions with no HLA LOH and no CDKN2A/B loss across the cohort. These findings suggest that despite abundant TILs, post ICI ES-TIL are associated with tumor clones that have evolved genetic mechanisms for evasion of the immune response (HLA LOH and/or CDKN2A/B loss). However, future mechanistic studies are needed to pinpoint the primary genomic event that transforms the ccRCC TME into a cold niche.

An adverse ccRCC TME is enriched stroma and myeloid signatures

We hypothesized that neoantigen depletion could be associated with a specific transcriptional signature, akin to those identified in clinical trial settings as biomarkers for response to ICI in ccRCC. To identify such a signature, we performed unsupervised Weighted Gene Co-expression Network Analysis (WGCNA) [60] to reconstruct modules from our transcriptomic samples similar to [11] (Fig. 6A). Reassuringly, we identified two gene expression modules #7 and #4 reflecting established microenvironmental features associated with therapeutic response in ccRCC: immune inflammatory response (“JAVELIN-like” signature) and “angiogenesis-like” (Fig. 6A, B). We next assessed the correlation between the expression of each WGCNA gene module and neoantigen depletion. While the JAVELIN-like and angiogenesis-like modules showed no association with neoantigen depletion, module 16 demonstrated the strongest association (Fig. 6A). Correlation analysis with previously known gene expression signatures illustrated that module 16 (which we refer to as an “Immune Escape” signature) was strongly associated with myeloid and stroma features of TME. The Immune Escape signature also resembled a recently described pan-cancer TGF-β signature derived in a previous study [45] which was linked to cancer-associated fibroblasts enriched in immune evasion and immunotherapy failure. However, no association between the Immune Escape signature and treatment status was observed (Wilcox P = 0.79) (Additional file 2: Fig. S10).

Fig. 6
figure 6

Immunoediting correlates with stroma and myeloid signatures. A) WGCNA identifies gene expression modules associated with inflammation (“JAVELIN-like”), angiogenesis, and Immune Escape. Gene dendrogram was first generated and then modules were extracted using dynamic tree cutting (top). Modules were annotated by comparing the correlation between the module eigengenes and previously known gene signatures describing different phenotypes (bottom). B Modules 7 (black), 4 (salmon), and 16 (magenta) are associated with previously described signatures, JAVELIN, angiogenesis, and myeloid/stroma. C scRNAseq demonstrates the cell type enrichment of Immune Escape signature in ccRCC patients. Different colors represent different cell types inferred from scRNAseq data. UMAP plot illustrates single cells collected from all 6 patients including treated and untreated patients. Computational extracted clusters were annotated as previously described [15]. D, E Association between ITH groups, Immune Escape signature, and N/S/ES pathologies. F These regions demonstrate an elevated Immune Escape gene signature in NIVO21. RA, RB, RC, RD, and RE denote different regions of a tumor sample

Full size image

To reveal the primary cellular populations driving the Immune Escape signature in the ccRCC TME, we leveraged scRNAseq from multiple tumor regions, lymph node, normal kidney, and peripheral blood of two ICI-naïve and four ICI-treated patients [16] (n = 167,283 single cells). We identified 28 clusters (Fig. 6C) using Louvain clustering [61, 62], and each cluster was annotated based on our previous study [16]. As expected, scRNAseq revealed enrichment of this signature in renal epithelium, tumor stroma, and tumor-associated macrophages (TAMs) and monocytes (Fig. 6C). Hence, scRNAseq, histopathological evaluation, and immunofluorescence (Additional file 2: Fig. S10B) further confirmed the association between Immune Escape and neoantigen depletion (Fig. 6A, Spearman correlation = 0.6), ITH (Fig. 6D, Wilcox P = 0.003), myeloid activation (Fig. 6B, Spearman correlation = 0.8, Additional file 2: Fig. S10B) and with stroma, and renal epithelium histopathology (Additional file 2: Fig. S9 and Fig. 6E, F).

Immune Escape correlates with clinical outcome to ICI therapy

Several previous studies have associated signatures of Immune Escape with poor clinical outcome in ICI-treated patients [63]. Thus, we evaluated whether our Immune Escape signature can correlate with clinical outcome to ICI treatment. We obtained publicly available RNAseq data for several clinical trials including phase 3 JAVELIN Renal 101 trial [11]—a phase III randomized anti-PD-L1 (avelumab) plus tyrosine kinase inhibitor (TKI, axitinib) versus multi-target TKI (sunitinib), IMmotion151 [64]—a phase III trial comparing anti-PDL1 (atezolizumab) plus anti-angeniogenesis agent (bevacizumab) versus TKI (sunitinib) in first-line metastatic renal cell carcinoma, CheckMate 009/010—a phase I/II, aPD-1 (nivolumab) treated, and CheckMate 025—a phase III randomized mTOR inhibitor (everolimus) versus aPD-1 [10]. We stratified patients by the median score (see “Methods”) of the 3 gene signatures obtained in our study (i.e., module 4/JAVELIN_like, 7/angiogenesis-like, and 16/immune escape).

The Immune Escape signature was strongly associated with the response to all three ICI regimens (avelumab plus axitinib HR = 1.53 P = 0.008, atezolizumab plus bevacizumab HR = 1.35 P = 0.019, and nivolumab HR = 1.45 P = 0.02, Fig. 7 and Additional file 2: Fig. S11). In contrast, the JAVELIN-like inflammatory signature was strongly associated with clinical outcome to avelumab plus axitinib (HR = 0.64 P = 0.006), but no association with clinical benefit was found between atezolizumab plus bevacizumab (HR = 0.82 P = 0.126) or nivolumab treatment (HR = 0.97 P = 0.823) (Fig. 7). Similarly, the angiogenesis-like signature was strongly correlated with the response to sunitinib in both IMmotion151 (HR = 0.48 P < 0.001) and JAVELIN Renal 101 (HR = 0.68 P = 0.008) as expected, but not with ICI-associated regimens. Associations between the Immune Escape signature and therapeutic response remained valid even when thresholds other than median were used to define immune escape high and low (Additional file 2: Fig. S12). Even though the Immune Escape signature was also associated with response to sunitinib in JAVELIN Renal 101, no association between sunitinib response or mTOR inhibition was observed in IMmotion151 and CheckMate 025. Overall, this analysis suggests that a transcriptional signature associated the tendency to lose putative neoantigens after ICI is associated with response to combination ICI therapy and nominates a new potential biomarker for this therapeutic regimen.

Fig. 7
figure 7

Association between immune escape and clinical outcome to checkpoint blockade. Survival analysis shows the association between gene signatures obtained in this study and clinical outcome of different independent retrospective trials. A Immune Escape and JAVELIN-like signatures are associated with PFS in patients treated Avelumab plus Axitinib in JAVELIN Renal 101 cohort. B Immune Escape signature, but not the JAVELIN-like signature, correlates with the response to atezolizumab plus bevacizumab in IMmotion151 but not JAVELIN signature. C Immune Escape signature, but not the JAVELIN-like signature, correlates with the efficacy of anti-PD1-treament in CheckMate 009, 010, 025

Full size image

Discussion

Here we used spatiotemporal, multimodal profiling to investigate the link tumor genomics, microenvironmental heterogeneity, peripheral immune response, and eventual immune escape in advanced and metastatic ccRCC. The fundamental discovery of our analysis is that ITH manifests well beyond the tumor genome and produces highly heterogeneous immune microenvironments in the tumor.

Emerging data on biomarkers of response to ICI in ccRCC has identified two potentially paradoxical observations: first, that TIL abundance alone is an insufficient predictor of ICI response [9], and second, that the presence of myeloid cells correlate with resistance to both ICI and anti-VEGF treatments. Strikingly, we observed that high myeloid score tumors were associated with neoantigen depletion which could, in principle, render ICI treatment ineffective. In agreement with this, we derived a transcriptomic signature associated with neoantigen depletion and Immune Escape, which was expressed in renal epithelium, tumor stroma as well as tumor-associated macrophages (TAMs) and monocytes. This Immune Escape signature was associated with response to several ICI regiments in published clinical trials. In total, these findings suggest that myeloid cells are associated with tumor clones that have evolved mechanisms to escape anti-tumor immune responses. Critically, such a hypothetic model requires detailed work and mechanistic validation in immunocompetent systems, which we are actively developing.

Why do regions with neoantigen depletion demonstrate elevation of myeloid cells but not cytotoxic T cells that would presumably eliminate tumor clones? Cancer immunoediting proceeds through three phases: elimination, equilibrium, and escape [65]. Throughout these phases, tumor immunogenicity evolves, and thereby, despite possible initial response to therapy, acquires immunosuppressive mechanisms that may enable disease progression. Our data suggests that myeloid-high, stroma-enriched (Fig. 6D), neoantigen-depleted tumor regions historically experienced a cytotoxic T cell response, which prompted the selection of tumor clones losing neoantigens and/or HLA/CDKN2A/B. Such a loss of target antigens through HLA LOH or neoantigen depletion would result in loss of antigen-TCR stimulation, leading to death of the corresponding neoantigen reactive T cells (Fig. 8). The selection of hydrophobic residues during the course of neoantigen depletion is another intriguing finding of our study. As demonstrated by others [49], hydrophobicity can result in an increased immunogenicity of the neoantigens. In agreement with this notion, we showed that the hydrophobic epitopes are the most frequently depleted neoantigens likely due to the higher immunogenicity. Likewise, HERVs were pruned by immune selection which might imply that HERVs are associated with immune recognition of HERV expressing tumor cells and thus, ICI-treated regions of the tumor may reduce HERV expression. Importantly, as with other findings in this analysis, the association between neoantigen loss and myeloid activation observed in our data remains purely correlative, and future studies will be necessary to mechanistically establish how immune evasion spatiotemporally evolves in ccRCC following ICI therapy.

Fig. 8
figure 8

A hypothetical model for spatiotemporal evolution of ccRCC links ITH to immune escape and adverse TME. Cancer cell death, potentially by cytotoxic killing, early in tumor evolution selects for tumor clones with HLA LOH and/or CDKN2A/B loss. This promotes the evolution of a TME depleted of antigen-specific T cells and enriched for myeloid cells

Full size image

Our multi-regional data also has significant implications for biomarker development. We demonstrated that TME markers of response such as JAVELIN and myeloid scores can be heterogenous within tumor regions (Fig. 2, Additional file 1: Additional file 2: Fig. S13). This underscores the importance of accounting for ITH when these signatures are used for patient selection for a specific therapy and longitudinal monitoring of therapies. Given recent data that ICI may have a role in adjuvant therapy following nephrectomy for high-risk disease, our data would suggest that several regions of the primary tumor should be sampled specially in the presence of ITH-associated genomic alterations (e.g., HLA LOH and CDKN2A/B loss). An intriguing finding was a trend towards lower ITH in ICI-treated tumors, even though this observation did not reach statistical significance. If validated in other studies, this in part can be attributed to outgrowth of few nonimmunogenic tumor subclones that managed to escape immune surveillance upon ICI treatment.

An important limitation of this study is that TME heterogeneity of metastatic disease was not assessed and may be less of an issue in biomarker development. Our study has several other potential limitations including its small sample size. To overcome this shortcoming, we validated several of our major findings in several independent cohorts. Another potential limitation of our study is the unavailability pre-treatment multi-regional sequencing data. However, inclusion of multi-regional data from 6 untreated patients allowed us to account for ITH in untreated tumors. Moreover, our neoadjuvant cohort was treated with single agent nivolumab over a short course which may not reflect the TME, and genomic changes induced by more potent combination strategies. Finally, we portrayed the characteristics of an adverse TME which may contribute to ICI resistance. Our study clearly demonstrates the interplay between genomic events and TME transformation from a cytotoxic to a cold immunophenotype. However, these findings remain purely an association of several contributing factors to ICI resistant and the exact causative hierarchy of events requires further investigation.

Conclusions

Our findings clearly suggest that the ccRCC genome and microenvironment co-evolve, and that loss of putative neoantigens (including SNVs, indels, and HERVs) is associated with a qualitatively myeloid-high environment and the loss of HLA and CDKN2A/B. These distinct genomic alterations are also associated with more peripheral changes, i.e., reduced T cell clonal diversity in the peripheral circulation. In conclusion, we find distinct genomic event enriched in immune escape tumor microenvironment in ccRCC both across and within tumors. Our findings have implications for future biomarker development for ICI response across ccRCC and other solid tumors.

Availability of data and materials

No new code was generated. All data generated in this study are provided in Additional file 1 and Additional file 2. Raw sequencing reads are available from the database of Genotypes and Phenotypes (dbGaP) phs003079.v1.p1 (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs003079.v1.p1) [66] for WES, WTS, and TCRseq.

Abbreviations

TCGA:

The Cancer Genome Atlas

WES:

Whole-exome sequencing

WTS:

Whole-transcriptome sequencing

SCNA:

Somatic copy number alterations

LOH:

Loss of heterozygosity

ITH:

Intratumoral heterogeneity

TCR:

T cell receptor

ccRCC:

Clear cell renal cell carcinoma

TME:

Tumor microenvironment

ICI:

Immune checkpoint inhibitor

TKI:

Tyrosine kinase inhibitor

HERV:

Human endogenous retrovirus

PBMC:

Peripheral blood mononuclear cells

TIL:

Tumor infiltrating lymphocytes

WGCNA:

Weighted gene co-expression network analysis

ssGSEA:

Single-sample gene set enrichment analysis

NKA:

Normal kidney 1 cm away from tumor

NKB:

Normal kidney 2 cm away from tumor

NKC:

Normal kidney 4 cm away from tumor

RA:

Region A of tumor (tumor center)

RB:

Region B of tumor (tumor superolateral)

RC:

Region C of tumor (tumor inferolateral)

RD:

Region D of tumor (tumor superomedial)

BX:

Before treatment tumor region

TC:

Tumor center

TN:

Tumor region nearest to kidney

TF:

Tumor region furthest from kidney

TEM:

Tumor medial

TLAT:

Tumor lateral

LN:

Lymph node

ADR:

Adrenal gland with tumor

TUP:

Tumor region upper

TLOW:

Tumor region lower

TMED:

Tumor region medial

References

  1. Carlo MI, Attalla K, Mazaheri Y, Gupta S, Yildirim O, Murray SJ, Coskey DT, Kotecha R, Lee CH, Feldman DR, Russo P. Phase II study of neoadjuvant nivolumab in patients with locally advanced clear cell renal cell carcinoma undergoing nephrectomy. Eur Urol. 2022;81:570–3.

    Article  CAS  Google Scholar 

  2. Liu X-D, Hoang A, Zhou L, Kalra S, Yetil A, Sun M, Ding Z, Zhang X, Bai S, German P. Resistance to antiangiogenic therapy is associated with an immunosuppressive tumor microenvironment in metastatic renal cell carcinoma. Cancer Immunol Res. 2015;3:1017–29.

    Article  CAS  Google Scholar 

  3. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160:48–61.

    Article  CAS  Google Scholar 

  4. Şenbabaoğlu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, Miao D, Ostrovnaya I, Drill E, Luna A. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17:1–25.

    Article  Google Scholar 

  5. Snyder A, Makarov V, Merghoub T, Yuan J, Zaretsky JM, Desrichard A, Walsh LA, Postow MA, Wong P, Ho TS. Genetic basis for clinical response to CTLA-4 blockade in melanoma. N Engl J Med. 2014;371:2189–99.

    Article  Google Scholar 

  6. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, Lu S, Kemberling H, Wilt C, Luber BS. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Sci. 2017;357:409–13.

    Article  CAS  Google Scholar 

  7. Samstein RM, Lee C-H, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, Barron DA, Zehir A, Jordan EJ, Omuro A. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51:202–6.

    Article  CAS  Google Scholar 

  8. Havel JJ, Chowell D, Chan TA. The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat Rev Cancer. 2019;19:133–50.

    Article  CAS  Google Scholar 

  9. McDermott DF, Huseni MA, Atkins MB, Motzer RJ, Rini BI, Escudier B, Fong L, Joseph RW, Pal SK, Reeves JA. Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Nat Med. 2018;24:749–57.

    Article  CAS  Google Scholar 

  10. Braun DA, Hou Y, Bakouny Z, Ficial M. Sant’Angelo M, Forman J, Ross-Macdonald P, Berger AC, Jegede OA, Elagina L: 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–18.

    Article  CAS  Google Scholar 

  11. Motzer RJ, Robbins PB, Powles T, Albiges L, Haanen JB, Larkin J, Mu XJ, Ching KA, Uemura M, Pal SK. Avelumab plus axitinib versus sunitinib in advanced renal cell carcinoma: Biomarker analysis of the phase 3 JAVELIN Renal 101 trial. Nat Med. 2020;26:1733–41.

    Article  CAS  Google Scholar 

  12. Motzer RJ, Escudier B, McDermott DF, George S, Hammers HJ, Srinivas S, Tykodi SS, Sosman JA, Procopio G, Plimack ER. Nivolumab versus everolimus in advanced renal-cell carcinoma. N Engl J Med. 2015;373:1803–13.

    Article  CAS  Google Scholar 

  13. Yarchoan M, Hopkins A, Jaffee EM. Tumor mutational burden and response rate to PD-1 inhibition. N Engl J Med. 2017;377:2500.

    Article  Google Scholar 

  14. Rini BI, Plimack ER, Stus V, Gafanov R, Hawkins R, Nosov D, Pouliot F, Alekseev B, Soulières D, Melichar B. Pembrolizumab plus axitinib versus sunitinib for advanced renal-cell carcinoma. N Engl J Med. 2019;380:1116–27.

    Article  CAS  Google Scholar 

  15. Turajlic S, Xu H, Litchfield K, Rowan A, Horswell S, Chambers T, O’Brien T, Lopez JI, Watkins TB, Nicol D. Deterministic evolutionary trajectories influence primary tumor growth: TRACERx renal. Cell. 2018;173(595–610):e511.

    Google Scholar 

  16. Krishna C, DiNatale RG, Kuo F, Srivastava RM, Vuong L, Chowell D, Gupta S, Vanderbilt C, Purohit TA, Liu M. Single-cell sequencing links multiregional immune landscapes and tissue-resident T cells in ccRCC to tumor topology and therapy efficacy. Cancer Cell. 2021;39:662.

    Article  CAS  Google Scholar 

  17. Au L, Hatipoglu E, de Massy MR, Litchfield K, Beattie G, Rowan A, Schnidrig D, Thompson R, Byrne F, Horswell S. Determinants of anti-PD-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell. 2021;39(1497–1518):e1411.

    Google Scholar 

  18. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  Google Scholar 

  19. Lawrence M, Huber W, Pages H, Aboyoun P, Carlson M, Gentleman R, Morgan MT, Carey VJ. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9:e1003118.

    Article  CAS  Google Scholar 

  20. Karolchik D, Baertsch R, Diekhans M, Furey TS, Hinrichs A, Lu Y, Roskin KM, Schwartz M, Sugnet CW, Thomas DJ. The UCSC genome browser database. Nucleic Acids Res. 2003;31:51–4.

    Article  CAS  Google Scholar 

  21. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Article  Google Scholar 

  22. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:1–11.

    Article  CAS  Google Scholar 

  23. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462:108–12.

    Article  CAS  Google Scholar 

  24. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7.

    Article  CAS  Google Scholar 

  25. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39:782–95.

    Article  CAS  Google Scholar 

  26. Golkaram M, Salmans ML, Kaplan S, Vijayaraghavan R, Martins M, Khan N, Garbutt C, Wise A, Yao J, Casimiro S. HERVs establish a distinct molecular subtype in stage II/III colorectal cancer with poor outcome. NPJ Genom Med. 2021;6:1–11.

    Article  Google Scholar 

  27. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.

    Article  CAS  Google Scholar 

  28. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.

    Article  CAS  Google Scholar 

  29. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, Del Angel G, Rivas MA, Hanna M. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43:491.

    Article  CAS  Google Scholar 

  30. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22:568–76.

    Article  CAS  Google Scholar 

  31. Kim S, Scheffler K, Halpern AL, Bekritsky MA, Noh E, Källberg M, Chen X, Kim Y, Beyter D, Krusche P. Strelka2: fast and accurate calling of germline and somatic variants. Nat Methods. 2018;15:591–4.

    Article  CAS  Google Scholar 

  32. Rimmer A, Phan H, Mathieson I, Iqbal Z, Twigg SR, Wilkie AO, McVean G, Lunter G. Integrating mapping-, assembly-and haplotype-based approaches for calling variants in clinical sequencing applications. Nat Genet. 2014;46:912–8.

    Article  CAS  Google Scholar 

  33. Larson DE, Harris CC, Chen K, Koboldt DC, Abbott TE, Dooling DJ, Ley TJ, Mardis ER, Wilson RK, Ding L. SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics. 2012;28:311–7.

    Article  CAS  Google Scholar 

  34. Ellrott K, Bailey MH, Saksena G, Covington KR, Kandoth C, Stewart C, Hess J, Ma S, Chiotti KE, McLellan M. Scalable open science approach for mutation calling of tumor exomes using multiple genomic pipelines. Cell Syst. 2018;6(271–281):e277.

    Google Scholar 

  35. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GR, Thormann A, Flicek P, Cunningham F. The ensembl variant effect predictor. Genome Biol. 2016;17:1–14.

    Article  Google Scholar 

  36. Amemiya HM, Kundaje A, Boyle AP. The ENCODE blacklist: identification of problematic regions of the genome. Sci Rep. 2019;9:1–5.

    Article  Google Scholar 

  37. Shen R, Seshan VE. FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing. Nucleic Acids Res. 2016;44:e131–e131.

    Article  Google Scholar 

  38. McGranahan N, Rosenthal R, Hiley CT, Rowan AJ, Watkins TB, Wilson GA, Birkbak NJ, Veeriah S, Van Loo P, Herrero J. Allele-specific HLA loss and immune escape in lung cancer evolution. Cell. 2017;171(1259–1271):e1211.

    Google Scholar 

  39. Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90.

    Article  CAS  Google Scholar 

  40. Raynaud F, Mina M, Tavernari D, Ciriello G. Pan-cancer inference of intra-tumor heterogeneity reveals associations with different forms of genomic instability. PLoS Genet. 2018;14:e1007669.

    Article  Google Scholar 

  41. Wu TD, Madireddi S, de Almeida PE, Banchereau R. Chen Y-JJ, Chitre AS, Chiang EY, Iftikhar H, O’Gorman WE, Au-Yeung A: Peripheral T cell expansion predicts tumour infiltration and clinical response. Nature. 2020;579:274–8.

    Article  CAS  Google Scholar 

  42. Vadim Nazarov, immunarch.bot, Eugene Rumynskiy. immunomind/immunarch: 0.6.5: Basic single-cell support (0.6.5). Zenodo. 2020. https://doi.org/10.5281/zenodo.3893991.

  43. Horn HS. Measurement of “overlap” in comparative ecological studies. Am Nat. 1966;100:419–24.

    Article  Google Scholar 

  44. Chowell D, Krishna C, Pierini F, Makarov V, Rizvi NA, Kuo F, Morris LG, Riaz N, Lenz TL, Chan TA. Evolutionary divergence of HLA class I genotype impacts efficacy of cancer immunotherapy. Nat Med. 2019;25:1715–20.

    Article  CAS  Google Scholar 

  45. Chakravarthy A, Khan L, Bensler NP, Bose P, De Carvalho DD. TGF-β-associated extracellular matrix genes link cancer-associated fibroblasts to immune evasion and immunotherapy failure. Nat Commun. 2018;9:1–10.

    Article  Google Scholar 

  46. Turajlic S, Xu H, Litchfield K, Rowan A, Chambers T, Lopez JI, Nicol D, O’Brien T, Larkin J, Horswell S. Tracking cancer evolution reveals constrained routes to metastases: TRACERx renal. Cell. 2018;173(581–594):e512.

    Google Scholar 

  47. Hakimi AA, Voss MH, Kuo F, Sanchez A, Liu M, Nixon BG, Vuong L, Ostrovnaya I, Chen Y-B, Reuter V. Transcriptomic profiling of the tumor microenvironment reveals distinct subgroups of clear cell renal cell cancer: data from a randomized phase III trial. Cancer Discov. 2019;9:510–25.

    Article  CAS  Google Scholar 

  48. Horn S, Leonardelli S, Sucker A, Schadendorf D, Griewank KG, Paschen A. Tumor CDKN2A-associated JAK2 loss and susceptibility to immunotherapy resistance. JNCI: J Nat Cancer Institute. 2018;110:677–81.

    Article  CAS  Google Scholar 

  49. Han G, Yang G, Hao D, Lu Y, Thein K, Simpson BS, Chen J, Sun R, Alhalabi O, Wang R. 9p21 loss confers a cold tumor immune microenvironment and primary resistance to immune checkpoint therapy. Nat Commun. 2021;12:1–19.

    Article  CAS  Google Scholar 

  50. Riley TP, Keller GL, Smith AR, Davancaze LM, Arbuiso AG, Devlin JR, Baker BM. Structure based prediction of neoantigen immunogenicity. Front Immunol. 2019;10:2047.

    Article  CAS  Google Scholar 

  51. Bonaventura P, Alcazer V, Mutez V, Tonon L, Martin J, Chuvin N, Michel E, Boulos RE, Estornes Y, Valladeau-Guilemond J. Identification of shared tumor epitopes from endogenous retroviruses inducing high-avidity cytotoxic T cells for cancer immunotherapy. Sci Adv. 2022;8:eabj3671.

    Article  CAS  Google Scholar 

  52. Panda A, de Cubas AA, Stein M, Riedlinger G, Kra J, Mayer T, Smith CC, Vincent BG, Serody JS, Beckermann KE. Endogenous retrovirus expression is associated with response to immune checkpoint blockade in clear cell renal cell carcinoma. JCI Insight. 2018;3:e121522.

    Article  Google Scholar 

  53. Smith CC, Beckermann KE, Bortone DS, De Cubas AA, Bixby LM, Lee SJ, Panda A, Ganesan S, Bhanot G, Wallen EM. Endogenous retroviral signatures predict immunotherapy response in clear cell renal cell carcinoma. J Clin Investig. 2019;128:4804–20.

    Article  Google Scholar 

  54. Nargund AM, Pham CG, Dong Y, Wang PI, Osmangeyoglu HU, Xie Y, Aras O, Han S, Oyama T, Takeda S. The SWI/SNF protein PBRM1 restrains VHL-loss-driven clear cell renal cell carcinoma. Cell Rep. 2017;18:2893–906.

    Article  CAS  Google Scholar 

  55. Hakimi AA, Attalla K, DiNatale RG, Ostrovnaya I, Flynn J, Blum KA, Ged Y, Hoen D, Kendall SM, Reznik E. A pan-cancer analysis of PBAF complex mutations and their association with immunotherapy response. Nat Commun. 2020;11:1–11.

    Article  Google Scholar 

  56. Liu X-D, Kong W, Peterson CB, McGrail DJ, Hoang A, Zhang X, Lam T, Pilie PG, Zhu H, Beckermann KE. PBRM1 loss defines a nonimmunogenic tumor phenotype associated with checkpoint inhibitor resistance in renal carcinoma. Nat Commun. 2020;11:1–14.

    CAS  Google Scholar 

  57. Zhou M, Leung JY, Gessner KH, Hepperla AJ, Simon JM, Davis IJ, Kim WY. PBRM1 inactivation promotes upregulation of human endogenous retroviruses in a HIF-dependent manner. Cancer Immunol Res. 2022;10:285.

    Article  CAS  Google Scholar 

  58. Zhang Y, Narayanan SP, Mannan R, Raskind G, Wang X, Vats P, Su F, Hosseini N, Cao X, Kumar-Sinha C: Single-cell analyses of renal cell cancers reveal insights into tumor microenvironment, cell of origin, and therapy response. Proc Nat AcadSci 2021;118.

  59. Zhang AW, McPherson A, Milne K, Kroeger DR, Hamilton PT, Miranda A, Funnell T, Little N, de Souza CP, Laan S. Interfaces of malignant and immunologic clonal dynamics in ovarian cancer. Cell. 2018;173(1755–1769):e1722.

    Google Scholar 

  60. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:1–13.

    Article  Google Scholar 

  61. Levine JH, Simonds EF, Bendall SC, Davis KL, El-ad DA, Tadmor MD, Litvin O, Fienberg HG, Jager A, Zunder ER. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell. 2015;162:184–97.

    Article  CAS  Google Scholar 

  62. Xu C, Su Z. Identification of cell types from single-cell transcriptomes using a novel clustering method. Bioinformatics. 2015;31:1974–80.

    Article  CAS  Google Scholar 

  63. Hegde PS, Chen DS. Top 10 challenges in cancer immunotherapy. Immunity. 2020;52:17–35.

    Article  CAS  Google Scholar 

  64. Motzer RJ, Banchereau R, Hamidi H, Powles T, McDermott D, Atkins MB, Escudier B, Liu L-F, Leng N, Abbas AR. Molecular subsets in renal cancer determine outcome to checkpoint and angiogenesis blockade. Cancer Cell. 2020;38(803–817):e804.

    Google Scholar 

  65. O’Donnell JS, Teng MW, Smyth MJ. Cancer immunoediting and resistance to T cell-based immunotherapy. Nat Rev Clin Oncol. 2019;16:151–67.

    Article  Google Scholar 

  66. Hakimi AA: Spatiotemporal evolution of the ccRCC microenvironment links intra-tumoral heterogeneity to immune escape CINOMA. phs003079.v1.p1, database of Genotypes and Phenotypes. 2022. https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs003079.v1.p1

Download references

Acknowledgements

We thank members of the Chan lab for their suggestions and critical reading of the manuscript.

Funding

We acknowledge funding sources including NIH R01 CA205426 (T.A.C.), NIH R35 CA232097 (T.A.C.), DOD grant KC180165, NIH R01 DE027738 (to L.G.T.M.), the NIH/NCI Cancer Center Support Grant P30 CA008748 (to MSKCC), P30 core grants (to MSKCC), Ludwig institute (A.A.H.), Weiss family fund (A.A.H.), Department of Defense (A.A.H.), and Illumina Inc.

Author information

Authors and Affiliations

  1. Illumina, Inc., 5200 Illumina Way, San Diego, CA, 92122, USA

    Mahdi Golkaram, Michael L. Salmans, Raakhee Vijayaraghavan, Chen Zhao, Rami Mehio, Shile Zhang, Jim Godsey, Traci Pawlowski & Li Liu

  2. Immunogenomics and Precision Oncology Platform, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Fengshen Kuo, Vlad Makarov, Kyle A. Blum, Renzo G. DiNatale, Timothy A. Chan & A. Ari Hakimi

  3. Department of Pathology, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Sounak Gupta, Ying-Bei Chen & Victor E. Reuter

  4. Department of Medicine, Genitourinary Oncology, Memorial Sloan Kettering Cancer Center, New York, New York, NY, 10065, USA

    Maria I. Carlo, Ritesh R. Kotecha, Robert J. Motzer & Martin H. Voss

  5. Human Oncology and Pathogenesis Program, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Cerise Tang, Vlad Makarov, Phillip Rappold, Kyle A. Blum & Renzo G. DiNatale

  6. Urology Service, Department of Surgery, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Kyle A. Blum, Renzo G. DiNatale, Paul Russo, Jonathan Coleman, Ed Reznik & A. Ari Hakimi

  7. Department of Surgery, Head & Neck Service, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Luc G. T. Morris

  8. Interventional Radiology, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Jeremy Durack

  9. Computational Oncology Service, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Ed Reznik

  10. Marie-Josée and Henry R. Kravis Center for Molecular Oncology, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Ed Reznik

  11. Department of Radiation Oncology, Memorial Sloan Kettering Cancer Center, New York, NY, 10065, USA

    Timothy A. Chan

  12. Center for Immunotherapy and Precision Immuno-Oncology, Cleveland Clinic, Cleveland, OH, 44195, USA

    Timothy A. Chan

  13. Lerner Research Institute, Cleveland Clinic, Cleveland, OH, 44195, USA

    Timothy A. Chan

  14. National Center for Regenerative Medicine, Cleveland Clinic, Cleveland, OH, 44195, USA

    Timothy A. Chan

Authors
  1. Mahdi Golkaram
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Fengshen Kuo
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Sounak Gupta
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Maria I. Carlo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Michael L. Salmans
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Raakhee Vijayaraghavan
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Cerise Tang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Vlad Makarov
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Phillip Rappold
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Kyle A. Blum
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Chen Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Rami Mehio
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Shile Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. Jim Godsey
    View author publications

    You can also search for this author in PubMed Google Scholar

  15. Traci Pawlowski
    View author publications

    You can also search for this author in PubMed Google Scholar

  16. Renzo G. DiNatale
    View author publications

    You can also search for this author in PubMed Google Scholar

  17. Luc G. T. Morris
    View author publications

    You can also search for this author in PubMed Google Scholar

  18. Jeremy Durack
    View author publications

    You can also search for this author in PubMed Google Scholar

  19. Paul Russo
    View author publications

    You can also search for this author in PubMed Google Scholar

  20. Ritesh R. Kotecha
    View author publications

    You can also search for this author in PubMed Google Scholar

  21. Jonathan Coleman
    View author publications

    You can also search for this author in PubMed Google Scholar

  22. Ying-Bei Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  23. Victor E. Reuter
    View author publications

    You can also search for this author in PubMed Google Scholar

  24. Robert J. Motzer
    View author publications

    You can also search for this author in PubMed Google Scholar

  25. Martin H. Voss
    View author publications

    You can also search for this author in PubMed Google Scholar

  26. Li Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  27. Ed Reznik
    View author publications

    You can also search for this author in PubMed Google Scholar

  28. Timothy A. Chan
    View author publications

    You can also search for this author in PubMed Google Scholar

  29. A. Ari Hakimi
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

M.G. performed data analysis and wrote the manuscript with input from all authors. M.G., A.A.H., T.A.C., and E.R. conceived the study and contributed to data interpretation. F.K., E.R., C.T., V.M., S.Z., C.Z., and R.M. assisted with analytical methodology development and bioinformatics support. M.L.S., R.V., L.L., T.P., and J.G. contributed to DNA, RNA, and TCR sequencing. R.J.M., M.I.C., P.R., J.C., L.G.T.M., and M.H.V. cared for patients analyzed in the study. K.A.B. and R.G.N. handled patient samples. R.G.N. collected and analyzed clinical metadata. S.G. and Y.C. analyzed H&E images. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Li Liu, Ed Reznik, Timothy A. Chan or A. Ari Hakimi.

Ethics declarations

Ethics approval and consent to participate

All patients involved in this study provided informed consent to participate. Informed consent and institutional review board approval were acquired at Memorial Sloan Kettering Cancer Center (MSK) (protocol number 16–195). This research conformed to the Declaration of Helsinki.

Consent for publication

Consent was obtained from all participants to publish supporting data including Additional file 1: Tables S1-S6.

Competing interests

M.G., R.V., M.L.S., J.G., T.P., R.M., C.Z., S.Z., and L.L. are current employees and shareholders of Illumina Inc. T.A.C. and L.G.T.M. are inventors on a patent held by Memorial Sloan Kettering related to the use of TMB in cancer immunotherapy. MSK has licensed the use of TMB for the identification of patients that benefit from immune checkpoint therapy to PGDx. L.G.T.M. reports laboratory research funding from AstraZeneca. T.A.C. is a co-founder of Gritstone Oncology and holds equity. T.A.C. holds equity in An2H. T.A.C. acknowledges grant funding from Bristol-Myers Squibb, AstraZeneca, Illumina, Pfizer, An2H, and Eisai. T.A.C. has served as an advisor for Bristol Myers Squibb, Illumina, Eisai, and An2H. R.J.M. reports consulting fees from Aveo, Calithera, Eisai, Eli Lilly, EMD Serono, Genentech, Merck, Novartis AG, Pfizer, and Roche, and contracted research to employer MSKCC for Bristol Myers Squibb, Eisai, Exelixis, Genentech, Merck, Pfizer, and Roche. A.A.H. is on the advisory board for Merck. M.H. receives commercial research grants from Bristol Myers Squibb, Pfizer, and Genentech/Roche, honoraria from Novartis, Bristol Myers Squibb, travel/accommodation from Astra Zeneca, Eisai, Novartis, and Takeda, and is a consultant/advisory board member for Alexion Pharmaceuticals, Aveo, Calithera Biosciences, Corvus Pharmaceuticals, Exelixis, Eisai, GlaxoSmithKline, Merck, Natera; Onquality Pharmaceuticals; Novartis and Pfizer. The remaining authors declare that they have no competing interests.

Additional information

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Li Liu, Ed Reznik, Timothy A. Chan, and A. Ari Hakimi are the senior authors.

A. Ari Hakimi is the lead contact.

Supplementary Information

Additional file1 Supplementary Tables

: Cohort characteristics, gene signature scores, ITH classification and histopathology review results

Additional file2 Supplementary Figures

: Provides additional data analysis supporting claims and conclusions drawn throughout the paper

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Golkaram, M., Kuo, F., Gupta, S. et al. Spatiotemporal evolution of the clear cell renal cell carcinoma microenvironment links intra-tumoral heterogeneity to immune escape. Genome Med 14, 143 (2022). https://doi.org/10.1186/s13073-022-01146-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13073-022-01146-3

Genome Medicine

ISSN: 1756-994X

Contact us