Metastasis remains the main reason for renal cell carcinoma (RCC)–associated mortality. Tyrosine kinase inhibitors (TKI) impart clinical benefit for most patients with RCC, but the determinants of response are poorly understood. We report an integrated genomic and transcriptomic analysis of patients with metastatic clear cell RCC (ccRCC) treated with TKI therapy and identify predictors of response. Patients in the COMPARZ phase III trial received first-line sunitinib or pazopanib with comparable efficacy. RNA-based analyses revealed four distinct molecular subgroups associated with response and survival. Characterization of these subgroups identified mutation profiles, angiogenesis, and macrophage infiltration programs to be powerful predictors of outcome with TKI therapy. Notably, predictors differed by the type of TKI received. Our study emphasizes the clinical significance of angiogenesis and immune tumor microenvironment and suggests that the critical effects its various aspects have on TKI efficacy vary by agent. This has broad implications for optimizing precision treatment of RCC.

Significance:

The determinants of response to TKI therapy in metastatic ccRCC remain unknown. Our study demonstrates that key angiogenic and immune profiles of the tumor microenvironment may affect TKI response. These findings have the potential to inform treatment personalization in patients with RCC.

This article is highlighted in the In This Issue feature, p. 453

Tyrosine kinase inhibitors (TKI) have revolutionized the treatment of metastatic renal cell carcinoma (RCC) and constitute a mainstay of care for patients with advanced disease (1). However, currently there are no reliable biomarkers to predict treatment response (2). The function of the tumor microenvironment (TME) in modifying the response to TKI therapy in metastatic RCC remains poorly understood. With the introduction of immune-checkpoint blockade, there is growing interest in unraveling the role of the TME to develop biomarkers that may help personalize treatment selection, understand resistance mechanisms, and inform the development of novel combination therapies.

Clear cell RCC (ccRCC), the most common and lethal form of RCC, is driven by distinct driver gene mutations. Inactivating mutations or methylation of von Hippel–Lindau (VHL) occur in all ccRCC tumors and lead to accumulation of hypoxia-inducible factor (HIF; ref. 3). HIF stabilization leads to activation of genes that promote angiogenesis (VEGF, PDGF). Proangiogenic tyrosine kinase receptors (VEGFR, PDGFR) are the primary targets for TKIs approved in this disease. Recently, other common molecular drivers of ccRCC have been identified. These mutations commonly involve chromatin-remodeling genes (PBRM1, KDM5C, SETD2, and BAP1; ref. 4). The predictive and prognostic significance of recurrent somatic mutations and gene-expression patterns in tumor tissue has previously been explored in TKI-treated patients, but an understanding of the determinants of treatment response is incomplete (5, 6), Here, we report the results of the largest molecular analysis of patients with TKI-treated metastatic ccRCC. We identify genomic and transcriptional determinants of treatment response, involving relevant aspects of angiogenesis and the immune TME. Furthermore, we propose molecular subgroups that reflect distinctly different survival outcomes after TKI therapy. Our findings set the stage for molecularly defined strategies to personalize therapy in metastatic ccRCC.

Gene-Expression Clustering Identifies Four Biologically Distinct Molecular Subgroups Associated with Differences in TKI Efficacy

We analyzed clinical outcomes and molecular features of archival tumor specimens prospectively collected for patients enrolled on the phase III COMPARZ trial (7). This international randomized study confirmed progression-free survival (PFS) noninferiority of first-line pazopanib versus sunitinib in patients with advanced ccRCC and reported comparable overall survival (OS) for both agents in a secondary analysis (7, 8).

We performed unsupervised consensus nonnegative matrix factorization (cNMF) clustering on the expression microarray data from 409 patients (212 sunitinib and 197 pazopanib treated) in the COMPARZ trial and identified four biologically distinct clusters based on the 1,500 most variable genes (Fig. 1A). The clusters were associated with significant differences in median OS (P = 2.00E–4) and PFS (P = 0.03), with cluster 4 demonstrating the worst OS (HR, 2.09; 95% CI, 1.47–2.97; P = 4.58E–5) and PFS (HR 1.54; 95% CI, 1.13–2.09; P = 5.72E–3; Fig. 1B) when compared with clusters 1–3. Cluster 4 demonstrated worse OS compared with each cluster individually. For PFS, cluster 4 was worse only when compared with cluster 3 (data not shown). Due to the similar median PFS and OS between clusters 1 and 3, we decided to compare cluster 4 with clusters 1 to 3 for the remainder of our analyses.

Figure 1.

Consensus NMF clustering identifies four biologically distinct clusters associated with different survival outcomes after TKI therapy. A, Unsupervised cNMF clustering from 412 patients identified four robust clusters based on the 1,500 most variable annotated genes. B, Kaplan–Meier curves depict OS (top) and PFS (bottom) by cluster. Censored data are indicated by vertical tick marks in the curves. All P values are calculated by the log-rank test. HR and CI values for OS and PFS were extracted from Cox proportional hazard regression models comparing cluster 4 (worst survival) to cluster 3 (best survival). B–D, Sample number per group is indicated below each graph. C, The percent distribution of clusters within IMDC risk groups. D, The percent distribution of ClearCode34 ccRCC molecular subtype within each cluster. E, Comparison of angiogenesis scores across clusters. TKI, tyrosine kinase inhibitor; cNMF, consensus nonnegative matrix factorization; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval; IMDC, International Renal Cell Carcinoma Database Consortium.

Figure 1.

Consensus NMF clustering identifies four biologically distinct clusters associated with different survival outcomes after TKI therapy. A, Unsupervised cNMF clustering from 412 patients identified four robust clusters based on the 1,500 most variable annotated genes. B, Kaplan–Meier curves depict OS (top) and PFS (bottom) by cluster. Censored data are indicated by vertical tick marks in the curves. All P values are calculated by the log-rank test. HR and CI values for OS and PFS were extracted from Cox proportional hazard regression models comparing cluster 4 (worst survival) to cluster 3 (best survival). B–D, Sample number per group is indicated below each graph. C, The percent distribution of clusters within IMDC risk groups. D, The percent distribution of ClearCode34 ccRCC molecular subtype within each cluster. E, Comparison of angiogenesis scores across clusters. TKI, tyrosine kinase inhibitor; cNMF, consensus nonnegative matrix factorization; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval; IMDC, International Renal Cell Carcinoma Database Consortium.

Close modal

The International Metastatic Database Consortium (IMDC) risk model is a validated tool used to prognosticate patients with metastatic RCC utilizing criteria including time from surgery to systemic therapy, performance status, and laboratory values (hemoglobin, calcium, neutrophil, and platelet counts). The number of adverse features is added as a score, which then serves to categorize patients as favorable, intermediate, or poor risk (9). Patients in the IMDC poor-risk group were enriched with cluster 4 (45.7%) compared with clusters 1 to 3 (Fisher exact test, P = 0.009; Fig. 1C). Furthermore, patients with progressive disease (PD), as defined by RECIST 1.0 criterion (10), were enriched with cluster 4 (χ2 test, P = 0.017) compared with those with stable disease (SD), complete response (CR), or partial response (PR; Supplementary Fig. S1A). We additionally assessed differentially expressed genes (DEG) and Ingenuity pathways analysis (IPA) to compare PD with CR/PR and found that patients with PD to TKI therapy demonstrated significant enrichment of upregulated genes involved in inflammatory pathways (IFNγ, IFNα, inflammatory response, IL6, and TNFα signaling; Supplementary Fig. S1B and S1C). These findings are further supported by the IPA results demonstrating enrichment in canonical pathways involved in IL8, macrophage phagocytosis, and enrichment in immune cell programs (Supplementary Fig. S1D–S1E).

The distribution of ClearCode34 (ccA/ccB) molecular subtypes, RNA signatures closely associated with ccRCC prognosis and previously validated in the metastatic setting (11, 12), was significantly different among the four clusters (Fisher exact test, P = 1.44E−30; Fig. 1D). Notably, cluster 3, which performed most favorably per comparison of PFS, OS, and IMDC distribution, was strongly enriched for ccA (89%), previously associated with improved survival outcomes and distinct from ccB by an upregulation of angiogenesis genes (FLT4, FLT1, VEGFB, ENG, KDR, and BAI1; ref. 11). To assess overall differences in angiogenesis gene expression among clusters, we utilized an established angiogenesis (Angio) gene-expression signature (13). Angiogenesis gene expression was significantly different among the four clusters (Kruskal–Wallis test, P = 2.20E−16), where cluster 3 had the highest angiogenesis gene-expression levels and cluster 1 had the lowest (Fig. 1E).

TP53 and BAP1 Mutations Are Enriched in Cluster 4

Next, we explored differences in the frequency of mutations in common ccRCC cancer driver genes among the clusters. Cluster 4 had similar overall mutational load compared with clusters 1 to 3 (Wilcoxon test, P = 0.12; Supplementary Fig. S1F). Somatic PBRM1 mutations, which previously demonstrated favorable prognostic effects among TKI-treated patients (5), were less frequent in cluster 4 (Fisher exact test, P = 0.003), whereas cluster 4 was enriched for the presence of TP53 (Fisher exact test, P = 0.03) and BAP1 (Fisher exact test, P = 0.05) mutations (Supplementary Fig. S1G), which conferred adverse outcomes in prior analyses (5, 14). There were no significant differences in SETD2, TERT promoter, KDM5C, or VHL mutation frequencies among clusters (data not shown).

Angiogenesis Gene-Expression Program Is Associated with TKI Response and Survival

Because significant differences in angiogenesis gene expression occurred among clusters, we explored its impact on TKI response and survival outcomes. Among the entire COMPARZ cohort, higher angiogenesis gene-expression levels were associated with improved objective response rate (by RECIST 1.0, ref. 10; Kruskal–Wallis test, P = 0.03; Fig. 2A). The association of angiogenesis gene expression with TKI response (by RECIST 1.0) was validated in the Beuselinck cohort (6) comprised of 53 patients with metastatic ccRCC treated with sunitinib (Wilcoxon test, P = 0.017; Fig. 2B). Furthermore, the COMPARZ Angiohi group (relative to the median angiogenesis score) demonstrated improved OS (HR, 0.68; 95% CI, 0.52–0.90; P = 6.11E−3) and PFS (HR, 0.68; 95% CI, 0.53–0.88; P = 2.49E−3) compared with the Angiolo group (Fig. 2C). Among IMDC risk groups, there was no significant difference in angiogenesis gene expression (Kruskal–Wallis test, P = 0.64; Fig. 2D). Within cluster 4 alone, the Angiohi group (relative to the median) did not demonstrate improvement in either OS (P = 0.05) or PFS (P = 0.08) when compared with the Angiolo group, although there was a similar trend (Supplementary Fig. S2A and S2B). Therefore, we hypothesized that angiogenesis programs alone do not explain the poor outcomes demonstrated in cluster 4.

Figure 2.

High angiogenesis gene expression is associated with improved TKI response and survival. A, Demonstrates objective response by RECIST 1.0 to TKI therapy based on angiogenesis score. B, Beuselinck cohort validation of objective response by RECIST 1.0 to TKI therapy based on angiogenesis score. C, Kaplan–Meier analyses demonstrating the impact of angiogenesis gene expression (Angiohi vs. Angiolo, based on median score) on OS (left) and PFS (right) among all patients in COMPARZ. D, Demonstrates angiogenesis score by IMDC risk group. E, Demonstrates angiogenesis score by mutation status of PBRM1 and BAP1. TKI, tyrosine kinase inhibitor; OS, overall survival; PFS, progression-free survival; IMDC, International Renal Cell Carcinoma Database Consortium; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models. Sample number per group indicated below each graph. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Figure 2.

High angiogenesis gene expression is associated with improved TKI response and survival. A, Demonstrates objective response by RECIST 1.0 to TKI therapy based on angiogenesis score. B, Beuselinck cohort validation of objective response by RECIST 1.0 to TKI therapy based on angiogenesis score. C, Kaplan–Meier analyses demonstrating the impact of angiogenesis gene expression (Angiohi vs. Angiolo, based on median score) on OS (left) and PFS (right) among all patients in COMPARZ. D, Demonstrates angiogenesis score by IMDC risk group. E, Demonstrates angiogenesis score by mutation status of PBRM1 and BAP1. TKI, tyrosine kinase inhibitor; OS, overall survival; PFS, progression-free survival; IMDC, International Renal Cell Carcinoma Database Consortium; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models. Sample number per group indicated below each graph. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Close modal

Next, we assessed the impact of somatic mutations on angiogenesis program expression levels. Tumors harboring PBRM1 mutations demonstrated higher angiogenesis gene expression (Wilcoxon test, P = 4.00E−4; Fig. 2E). These findings are consistent with recent in vitro and murine studies, which suggest that PBRM1 inactivation leads to further upregulation of HIF1 via STAT3 (15, 16), as well as a recent report by McDermott and colleagues (17). In line with our earlier findings, we found that patients in cluster 3 had the highest frequency of PBRM1 mutations (54.2%) compared with other clusters (χ2 test, P = 0.009). Patients with BAP1 mutations had lower angiogenesis gene expression (Wilcoxon test, P = 0.01; Fig. 2E). SETD2, KDM5C, TERT promoter, and TP53 mutations were not associated with significant differences in angiogenesis gene expression (data not shown). We do not report the effect of VHL loss on angiogenesis gene score, as data regarding methylation status of VHL were not available. Utilizing The Cancer Genome Atlas (TCGA) KIRC cohort (4), we validated the association of BAP1 mutation status (Wilcoxon test, P = 4.0E−6) with lower angiogenesis gene expression and found that PBRM1-mutant tumors demonstrated a trend toward higher angiogenesis scores; however, this did not reach statistical significance (Wilcoxon test, P = 0.12; Supplementary Fig. S2C and S2D). In summary, PBRM1 and BAP1 mutational status correlated with angiogenesis gene expression and, overall, patients in the Angiohi group demonstrated improved TKI response, OS, and PFS. However, the association between angiogenesis and survival is weakest in cluster 4.

Tumor Immune Infiltration Is Associated with Differences in TKI Efficacy

Recognizing the possible role of the tumor microenvironment (TME) in response to targeted therapy, we next performed pathway analyses to understand the underlying pathways in cluster 4, which exhibited the worst outcomes. We also used gene set enrichment analysis (GSEA) to identify DEGs between clusters 4 and 1–3. Consistent with the aggressive nature of these tumors, we saw enrichment of MYC targets, cell cycle and proliferation programs among genes upregulated in cluster 4. Strikingly, we identified enrichment for several hallmark inflammation signatures, such as inflammatory and IFNγ responses (Fig. 3A; Supplementary Fig. S3A). Using ESTIMATE (18), we found that cluster 4 had the highest immune score, a marker of total immune infiltration, among all clusters (Kruskal–Wallis test, P = 2.20E−16; Fig. 3B). Using IHC, PD-L1 positivity among tumors was found more frequently in cluster 4, where 60% of tumors demonstrated PD-L1 expression on tumor cells, versus clusters 1 to 3, where 66% or more were PD-L1 null (Fisher exact test, P = 1.71E−8; Fig. 3C). These findings suggest that cluster 4 is characterized by an immune-infiltrated and suppressed TME.

Figure 3.

Cluster 4 demonstrates enrichment in inflammatory pathways and macrophage infiltration. A, GSEA of hallmark gene sets comparing cluster 4 vs. clusters 1–3. Enrichment scores are ranked and colored based on the NES and sized by the log10 transformed value of the adjusted P value. B, Comparison of ESTIMATE immune score within each cluster. C, Differences in the proportion of overall PD-L1 tumoral positivity by IHC in each cluster. P value was derived using the Fisher exact test. D, ssGSEA immune deconvolution of cluster 4 vs. clusters 1–3 with the mean infiltration differences noted on the x-axis and specific immune populations on the y-axis. The size of the circles represents the log of the FDR and color represents the directionality of the association. E, Kaplan–Meier analysis demonstrating the impact of macrophage infiltration (Macrophagehi vs. Macrophagelo, based on median score) on OS among all patients in COMPARZ. F, Demonstrates differences in macrophage infiltration by IMDC risk group. NES, normalized enrichment score; OS, overall survival; PFS, progression-free survival; GSEA, gene set enrichment analysis; ssGSEA, single-sample GSEA; IHC, immunohistochemistry; FDR, false discovery rate; HR, hazard ratio; CI, confidence interval; IMDC, International Renal Cell Carcinoma Database Consortium. All HR and CI values for OS were extracted from Cox proportional hazard regression models. Sample number per group is indicated below each graph.

Figure 3.

Cluster 4 demonstrates enrichment in inflammatory pathways and macrophage infiltration. A, GSEA of hallmark gene sets comparing cluster 4 vs. clusters 1–3. Enrichment scores are ranked and colored based on the NES and sized by the log10 transformed value of the adjusted P value. B, Comparison of ESTIMATE immune score within each cluster. C, Differences in the proportion of overall PD-L1 tumoral positivity by IHC in each cluster. P value was derived using the Fisher exact test. D, ssGSEA immune deconvolution of cluster 4 vs. clusters 1–3 with the mean infiltration differences noted on the x-axis and specific immune populations on the y-axis. The size of the circles represents the log of the FDR and color represents the directionality of the association. E, Kaplan–Meier analysis demonstrating the impact of macrophage infiltration (Macrophagehi vs. Macrophagelo, based on median score) on OS among all patients in COMPARZ. F, Demonstrates differences in macrophage infiltration by IMDC risk group. NES, normalized enrichment score; OS, overall survival; PFS, progression-free survival; GSEA, gene set enrichment analysis; ssGSEA, single-sample GSEA; IHC, immunohistochemistry; FDR, false discovery rate; HR, hazard ratio; CI, confidence interval; IMDC, International Renal Cell Carcinoma Database Consortium. All HR and CI values for OS were extracted from Cox proportional hazard regression models. Sample number per group is indicated below each graph.

Close modal

Next, we performed immune deconvolution analyses via single-sample GSEA (ssGSEA; ref. 19) to better characterize the tumor immune microenvironment in cluster 4 compared with clusters 1 to 3. Cluster 4 demonstrated enrichment for many immune populations (Fig. 3D), with macrophages showing the most significant (Wilcoxon test, q = 8.52E−33) infiltration difference and GSEA score [normalized enrichment score (NES) = 2.33, Kolmogorov–Smirnov test, P = 0.0015] compared with clusters 1 to 3 (Supplementary Fig. S3B). DEGs of cluster 4 compared with clusters 1 to 3 also demonstrated enrichment for immune and macrophage-related genes (Supplementary Fig. S3C). We then utilized dual macrophage (CD68+) and PD-L1 IHC staining to orthogonally validate this finding. Cluster 4 demonstrated a higher proportion of PD-L1–positive macrophages compared with clusters 1 to 3 (Fisher exact test, P = 3.50E−7; Supplementary Fig. S3D). Supplementary Table S1 summarizes the immune deconvolution results by cluster for comparison.

We then evaluated the impact of macrophage enrichment on survival within the COMPARZ cohort (regardless of cluster) using ssGSEA. We observed significantly worse OS (HR, 1.54; 95% CI, 1.17–2.03; P = 1.98E−3) among subjects with high macrophage infiltration (Macrophagehi, relative to the median macrophage score; Fig. 3E). Furthermore, patients with PD, as defined by RECIST 1.0 criteria, had significantly higher macrophage infiltration (Kruskal–Wallis test, P = 0.02; Supplementary Fig. S3E). However, there was no significant difference in PFS between high and low macrophage infiltration subject groups, although the same trend was found (HR, 1.22; 95% CI, 0.95–1.56; P = 0.11; Supplementary Fig. S3F). Next, we hypothesized that type 2 macrophages, which are considered to be immunosuppressive (20), are responsible for driving these associations with TKI response. Using CIBERSORT (21), we found that a high type 2 macrophage infiltration (M2hi, relative to the median M2 score) was associated with poor OS (HR, 1.38; 95% CI, 1.06–1.81; P = 0.019) and PFS (HR 1.40; 95% CI, 1.09–1.78; P = 7.90E−3) compared with the M2lo group (Supplementary Fig. S3G). This association was attenuated when we evaluated the impact of both macrophage types (M1hiM2hi, relative to the median M1M2 score) on OS (HR, 1.29; 95% CI, 0.99–1.70; P = 0.063) and PFS (HR, 1.29; 95% CI, 1.01–1.66; P = 0.040; Supplementary Fig. S3H). Of note, somatic mutations in BAP1, PBRM1, and TP53 were not associated with significant differences in macrophage infiltration (data not shown).

Utilizing ESTIMATE immune score (18), total immune infiltration was similar among IMDC risk groups (Kruskal–Wallis test, P = 0.11; Supplementary Fig. S3I). However, the extent of macrophage infiltration by ssGSEA score increased with worsening IMDC risk grouping (Kruskal–Wallis test, P = 0.012; Fig. 3F), with the poor risk group demonstrating the highest macrophage score. In summary, cluster 4 demonstrates an immune-infiltrated and regulated TME that is highly infiltrated with PD-L1–positive macrophages. Furthermore, Macrophagehi infiltration alone is associated with worse OS and decreased objective response in TKI-treated patients.

Flow Cytometry and Immunostaining Analysis of ccRCC Tumors Prospectively Collected at MSKCC

We utilized a cohort of 12 patients with ccRCC whose primary tumor specimens were collected prospectively and analyzed via flow cytometry. Our gating strategy is demonstrated in Supplementary Fig. S4A. All patients went on to receive first-line TKI therapy (5 pazopanib and 7 sunitinib) for metastatic disease. We confirmed a trend toward a higher proportion of activated tumor-associated macrophage (CD45+HLADR+CD3CD19CD14+CD16+) infiltration in TKI nonresponders compared with responders (Wilcoxon one-tailed test, P = 0.052; Fig. 4A and B). In keeping with our genomic analyses, immunofluorescence staining of tumors from nonresponders was more immune-infiltrated overall compared with responders, and in areas with the highest CD45+ cell infiltrates, we noticed a higher density of CD68+ macrophages in nonresponders (Fig. 4C).

Figure 4.

IHC and IF validation of TKI response and macrophage infiltration. A, Macrophage infiltration (CD16+CD14+) by TKI response (responder defined as time to treatment failure of >6 months vs. ≤6 months for nonresponders) in the MSKCC cohort. B, Representative flow cytometry results demonstrating higher macrophage (CD16+CD14+) infiltration in a TKI nonresponder vs. responder in the MSKCC cohort. C, Representative immunofluorescence demonstrating the difference in overall immune (CD45+) and macrophage infiltration (CD68+) in a TKI nonresponder (RCC540) and responder (RCC563) in the MSKCC cohort. D, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage score grouping on OS (left) and PFS (right) among all patients in COMPARZ. TKI, tyrosine kinase inhibitor; IF, immunofluorescence; IHC, immunohistochemistry; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models.

Figure 4.

IHC and IF validation of TKI response and macrophage infiltration. A, Macrophage infiltration (CD16+CD14+) by TKI response (responder defined as time to treatment failure of >6 months vs. ≤6 months for nonresponders) in the MSKCC cohort. B, Representative flow cytometry results demonstrating higher macrophage (CD16+CD14+) infiltration in a TKI nonresponder vs. responder in the MSKCC cohort. C, Representative immunofluorescence demonstrating the difference in overall immune (CD45+) and macrophage infiltration (CD68+) in a TKI nonresponder (RCC540) and responder (RCC563) in the MSKCC cohort. D, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage score grouping on OS (left) and PFS (right) among all patients in COMPARZ. TKI, tyrosine kinase inhibitor; IF, immunofluorescence; IHC, immunohistochemistry; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models.

Close modal

Combined Effects of Angiogenesis and Macrophage Infiltration on Survival Among TKI-Treated Patients

Next, we sought to integrate the notable associations we observed for aspects of angiogenesis (angiogenesis score) and the immune TME (macrophage infiltration) and set out to compare four groups based on angiogenesis and macrophage enrichment scores, categorizing patients individually as either Angiolo/Angiohi and Macrophagelo/Macrophagehi. Strikingly, these groups demonstrated significant differences in median OS (P = 1.74E−5) and PFS (P = 0.0006; Fig. 4D). Patients in the AngioloMacrophagehi group demonstrated the worst outcomes compared with the AngiohiMacrophagelo group, which demonstrated the best survival outcomes, for OS (HR, 3.12; 95% CI, 1.93–5.03; P = 2.91E−6) and PFS (HR, 2.27; 95% CI, 1.51–3.42; P = 8.58E−5; Fig. 4D). These findings suggest that angiogenesis and macrophage infiltration may be important mechanisms defining TKI response. Not surprisingly, cluster 4 is enriched with the Angiolo/hiMacrophagehi groups and cluster 3 demonstrates enrichment in AngiohiMacrophagelo/hi groups (χ2 test, P = 3.37E−44; Supplementary Fig. S4B). Similarly, the poor-risk IMDC group was also enriched with Angiolo/hiMacrophagehi groups (Supplementary Fig. S4C). Furthermore, colony-stimulating factor-1 receptor (CSF1R) expression, involved in the development, survival, and proliferation of monocytes and macrophages (22), is highest among Macrophagehi groups (Supplementary Fig. S4D). Within the non–TKI-treated KIRC TCGA cohort, we demonstrated that angiogenesis score (Supplementary Fig. S4E) is prognostic of cancer-specific survival (CSS); however, macrophage infiltration is not (Supplementary Fig. S4F) associated with CSS.

Integrated Genomic Model Improves Prediction of Clinical Benefit in Patients with TKI-Treated ccRCC

An integrated analysis of our data is shown in Fig. 5A. Cluster 4 is uniquely enriched with inflammasome and myeloid gene expression and cluster 3 is enriched with the angiogenesis expression program and has moderate immune-related gene expression, whereas clusters 1 and 2 appear to have decreased expression of both angiogenesis and immune pathways. We then sought to assess how our genomic findings can improve the performance of the IMDC risk prognostic model to stratify patients receiving first-line TKI therapy. Utilizing standard IMDC risk grouping alone for the COMPARZ biomarker population (n = 409) predicted 2-year OS and PFS with a c-index of 0.63 and 0.60, respectively. Using a stepwise Cox proportional hazard model, we determined the molecular variables that demonstrate prognostic value (P < 0.1) in the presence of IMDC risk grouping and included TP53 and PBRM1 mutational status, angiogenesis score, and macrophage infiltration (Supplementary Table S2). Utilizing molecular variables alone, incorporating angiogenesis, macrophage infiltration, and PBRM1 and TP53 mutational status achieved a c-index of 0.63 for OS. Angiogenesis and macrophage infiltration were the strongest predictors of PFS, with a c-index of 0.61. Notably, integration of molecular and IMDC variables improved the c-index for OS from 0.63 to 0.69 and PFS from 0.60 to 0.65 (Fig. 5B).

Figure 5.

Summary oncoprint highlights the immune infiltration within cluster 4 compared with others and combining genomic markers with IMDC improves survival prediction. A (Top heat map), Angiogenesis and macrophage scores, PD-L1 tumoral expression positivity by IHC, mutational status, the best response to TKI therapy by RECIST 1.0, IMDC and mutational status by cluster. Bottom heat map demonstrates angiogenesis, immune and antigen-presenting machinery, and inflammasome and myeloid gene-expression differences by cluster. APM, antigen-presenting machinery. B, Multivariable model and c-index with the addition of genomic markers for OS (left) and PFS (right). OS, Overall survival; PFS, progression-free survival; IMDC, International Renal Cell Carcinoma Database Consortium; HR, hazard ratio; IHC, immunohistochemistry. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Figure 5.

Summary oncoprint highlights the immune infiltration within cluster 4 compared with others and combining genomic markers with IMDC improves survival prediction. A (Top heat map), Angiogenesis and macrophage scores, PD-L1 tumoral expression positivity by IHC, mutational status, the best response to TKI therapy by RECIST 1.0, IMDC and mutational status by cluster. Bottom heat map demonstrates angiogenesis, immune and antigen-presenting machinery, and inflammasome and myeloid gene-expression differences by cluster. APM, antigen-presenting machinery. B, Multivariable model and c-index with the addition of genomic markers for OS (left) and PFS (right). OS, Overall survival; PFS, progression-free survival; IMDC, International Renal Cell Carcinoma Database Consortium; HR, hazard ratio; IHC, immunohistochemistry. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Close modal

The Effect of Angiogenesis and Macrophage Infiltration on TKI Response May Differ by TKI Received

TKIs for ccRCC are utilized for their potent inhibition of kinases responsible for tumor angiogenesis (VEGFR2). However, TKIs also have significant off-target effects that may affect therapeutic efficacy. Notably, sunitinib is a more potent inhibitor of kinases involved in hematopoiesis (c-Kit, Flt3, and CSF1R) compared with pazopanib (22–25). In vivo and in vitro, sunitinib has been shown to decrease the concentration of myeloid-derived suppressor cells (MDSC) and reverse T-cell suppression by IFNγ in the TME (26). The effects of pazopanib on the immune TME are less well understood and may differ, based on differences in target kinome. Utilizing ssGSEA immune deconvolution analysis, we observed significant differences in the interaction between TKI response and TME profiles when comparing patients treated with pazopanib versus sunitinib. Sunitinib-treated responders, defined as a CR or PR, had increased angiogenesis and cytotoxic T-cell and CD8+ T-cell infiltration compared with those who had PD. By contrast, pazopanib-treated responders had lower macrophage infiltration compared with those who had PD, whereas angiogenesis score did not seem to discriminate TKI response (Fig. 6A). The combination of angiogenesis and macrophage infiltration significantly discriminates survival differences among pazopanib-treated patients, with AngioloMacrophagehi demonstrating the worst OS (HR, 1.86; 95% CI, 1.00–3.46; P = 0.05) and PFS (HR, 6.01; 95% CI, 2.73–13.20; P = 8.00E−6) compared with the AngiohiMacrophagelo group. However, among patients treated with sunitinib, the combination did not significantly discriminate OS (P = 0.08) or PFS (P = 0.07; Fig. 6B and C). We tested angiogenesis alone by drug type and found that the Angiohi group was associated with improved OS (HR, 1.54; 95% CI, 1.07–2.25; P = 0.022) and PFS (HR 1.56; 95% CI, 1.10–2.22; P = 0.012) in the sunitinib cohort only (Supplementary Fig. S5A and S5B). On the contrary, the Macrophagehi group had significantly worse OS (HR, 2.62; 95% CI, 1.71–4.00; P = 8.56E−6) and PFS (HR 1.85; 95% CI, 1.30–2.63; P = 6.74E−4) compared with the Macrophagelo group in the pazopanib cohort only, and survival did not differ between the sunitinib Macrophage groups (Supplementary Fig. S5C and S5D). In summary, the different components of the TME may affect response to the specific TKI used for treatment.

Figure 6.

Immune infiltration differences and TKI response differ by specific type of TKI received. A, Differences in angiogenesis and immune infiltration by drug (sunitinib vs. pazopanib) and response to therapy (CR/PR vs. PD) where the x-axis demonstrates specific immune cell populations. Differences in infiltration within each drug category are represented by the size (log10 transformed nominal P value) and color (difference in mean ssGSEA score) of the circles; asterisks represent significant differences in infiltration between drug categories (P < 0.05). B, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage scores (Angiogenesishi/loMacrophagehi/lo, relative to the median) on OS among patients treated with sunitinib (left) and pazopanib (right). C, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage scores (Angiogenesishi/loMacrophagehi/lo, relative to the median) on PFS among patients treated with sunitinib (left) and pazopanib (right). TKI, tyrosine kinase inhibitor; ssGSEA, single-sample gene set enrichment analysis; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Figure 6.

Immune infiltration differences and TKI response differ by specific type of TKI received. A, Differences in angiogenesis and immune infiltration by drug (sunitinib vs. pazopanib) and response to therapy (CR/PR vs. PD) where the x-axis demonstrates specific immune cell populations. Differences in infiltration within each drug category are represented by the size (log10 transformed nominal P value) and color (difference in mean ssGSEA score) of the circles; asterisks represent significant differences in infiltration between drug categories (P < 0.05). B, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage scores (Angiogenesishi/loMacrophagehi/lo, relative to the median) on OS among patients treated with sunitinib (left) and pazopanib (right). C, Kaplan–Meier analyses demonstrating the impact of angiogenesis and macrophage scores (Angiogenesishi/loMacrophagehi/lo, relative to the median) on PFS among patients treated with sunitinib (left) and pazopanib (right). TKI, tyrosine kinase inhibitor; ssGSEA, single-sample gene set enrichment analysis; OS, overall survival; PFS, progression-free survival; HR, hazard ratio; CI, confidence interval. All HR and CI values for PFS and OS were extracted from Cox proportional hazard regression models. RECIST 1.0 objective response is categorized as PD, progressive disease; CR, complete response; PR, partial response; or SD, stable disease.

Close modal

We performed a comprehensive molecular analysis integrating somatic mutations, RNA, and protein expression with clinical outcomes to evaluate predictors of TKI efficacy in the largest cohort of metastatic RCC reported to date. Archival tumor specimens were collected prospectively in patients receiving pazopanib versus sunitinib on the COMPARZ phase III trial. We identified four distinct molecular subgroups that differed significantly in response and survival. Detailed characterization of these clusters emphasized the central role of the TME for the outcome with TKI therapy and identified angiogenesis and macrophage infiltration as critical determinants of TKI response. Independent of IMDC risk category, we found a superior outcome for patients with higher angiogenesis scores (HR for PFS and OS 0.68 and 0.68, respectively). Although these associations may merely constitute a prognostic signal, we also detected association with objective response (Fig. 2A), a finding validated in the Beuselinck validation cohort (Fig. 2B), suggesting that angiogenesis score may indeed be associated with VEGFR-directed TKI response. Varying effects on the TME, particularly the upregulation of angiogenesis, may also underlie the different clinical impact of mutation status for PBRM1 and BAP1, two of the most commonly altered genes in this disease after VHL and cumulatively altered in >50% of patients with metastatic ccRCC. Loss-of-function mutations in the two genes correlated inversely with outcomes in patients treated on COMPARZ, confirming findings in other TKI-treated data sets (5, 14). We noted upregulation and suppression of angiogenesis observed with loss-of-function mutations in PBRM1 and BAP1, respectively, providing a plausible explanation for the different effects observed in clinical behaviors.

Critically, these findings are consistent with a recent report by McDermott and colleagues (17), in which they evaluated the randomized phase II Immotion150 study of atezolizumab (anti–PD-L1) alone or combined with bevacizumab (anti-VEGF) versus sunitinib in 305 patients with treatment-naïve metastatic RCC. Patients with high angiogenesis gene score demonstrated a PFS advantage in the angio-high arm (HR, 0.31; 95% CI, 0.19–0.55; P < 0.001). Similar to our analyses on COMPARZ, the authors also noted that PBRM1-mutated tumors also had higher angiogenic gene expression. The IMDC prognostic model has been multiply validated and is broadly applied in standard practice and clinical trial design (9). Although it encompasses some “host” factors that indirectly reflect systemic inflammatory effects of the cancer (e.g., anemia, neutrophilia, and thrombocytosis), it is agnostic to any molecular features of the individual cancer that may underlie such changes. With an increasing understanding of tumor biology and our growing ability to decipher the details of such on a molecular level, one should argue that integration of such information is the logical next step to improve our ability to prognosticate patients, possibly to guide the rational choice of agents.

In our analysis, cluster 4 was associated with the worst survival and enriched in the IMDC poor-risk group (46%). This cluster was characterized by a high frequency of BAP1 and TP53 mutations, moderate angiogenesis, high immune infiltration, and notably higher proportion of cases with PD-L1 expression on tumor cells by IHC when compared with clusters 1 to 3. BAP1 mutations have been associated with poor outcomes among patients with ccRCC (27, 28). Much of the above suggests notable overlap between cluster 4 as defined here and an “inflamed subtype” of RCC proposed in a recent report by Wang and colleagues (29), who used an elegant approach across various xenograft models and via transcriptomic analyses defined a hyperinfiltrated molecular variant of RCC, enriched for aggressive phenotypes, including BAP1-deficient tumors.

Further, cluster 4 is similar to the ccrcc4 group identified by Beuselinick and colleagues (6), which was associated with sarcomatoid differentiation, upregulation of MYC and cellular immune pathways, and poor response to sunitinib. Cluster 3 demonstrated the best survival and comprised about one third of the IMDC favorable and intermediate risk groups. This cluster was characterized by enrichment in PBRM1 mutations, high angiogenesis, and moderate immune infiltration. This cluster appears to be similar to the ccA molecular subtype reported by Brannon and colleagues (11) and the ccrcc2 cluster reported by Beuselinck and colleagues (6), both of which are characterized by high angiogenesis and improved survival outcomes.

Independent of clusters, we investigated differences in the angiogenesis program and immune TME across all subjects in COMPARZ. We demonstrated that enrichment in angiogenesis gene expression was associated with improved TKI response, OS, and PFS. These findings are consistent with results from Beuselinck and colleagues (13) where high angiogenesis (VEGFA, VEGFR1, VEGFR2, VEGFR2, and HIF2A) was associated with improved response to sunitinib. Within the immune TME, the strongest signal was seen for tumor-infiltrating macrophages. Due to a lack of available gene signatures for MDSCs, we were unable to assess this population using ssGSEA immune deconvolution (30). In the overall cohort, the Macrophagehi group demonstrated worse PFS and was enriched among the IMDC poor-risk group. As a resistance mechanism, tumors may express the chemokine colony-stimulating factor-1 (CSF1) to help recruit peripheral monocytes that eventually differentiate into tumor-associated macrophages (TAM; ref. 31). TAMs also produce VEGF and other angiogenic proteins that may sustain angiogenesis and promote a state of immunosuppression that may promote TKI resistance (31–33).

With angiogenesis score and macrophage infiltration in the TME showing the strongest associations with survival outcomes, we sought to better understand the combined effect. Our multivariate Cox regression analysis demonstrated that loss-of-function TP53 and PBRM1 mutations, angiogenesis, and macrophage infiltration were independently prognostic for OS and PFS after adjusting for IMDC classifications. Utilizing a combination of these molecular markers, we improved our ability to stratify OS and PFS (Fig. 5B). Specifically, the incorporation of angiogenesis and macrophage infiltration improved the baseline c-index for OS from 0.63 with IMDC risk stratification alone to 0.69 after the addition of molecular markers. This is in comparison with a recent report demonstrating a c-index of 0.63 after the addition of ccA and ccB molecular subtypes (11) to a baseline IMDC model with a c-index of 0.60 (12). Importantly, we also found different associations between treatment and survival when we considered angiogenesis and macrophage infiltration.

Although the primary analysis of the COMPARZ trial suggested similar efficacy for both TKIs (7), we demonstrated distinctly different interactions with the TME for each agent. High macrophage infiltration was a strong adverse prognostic factor in pazopanib-treated patients. In contrast, we found that macrophage infiltration among sunitinib-treated patients did not differentiate survival outcomes. This is consistent with the known roles of sunitinib in reducing accumulation of immunosuppressive MDSCs and T regulatory cells within primary tumors and thus enhancing the infiltration and cytotoxicity of CD8+ and CD4+ tumor-infiltrating lymphocytes (26, 34, 35). These findings were also validated in the Beuselinck cohort, where sunitinib-treated patients demonstrated differences in objective response by angiogenesis alone and not macrophage infiltration. In contrast, reduced macrophage and Th2 gene signatures correlate with improved survival in pazopanib-treated patients. The mechanism of action of pazopanib is poorly elucidated, but our results prompt questions about how these anti-inflammatory populations might counteract the antitumor functions of pazopanib. These findings warrant further investigation and should motivate separate analyses for other TKIs routinely used in this disease; ultimately, they could inform the rational choice of agent based on baseline analysis of tumor tissue. Our results give biological insights into possible mechanisms by which the TME might affect the efficacy of specific drugs. Furthermore, these findings need to be validated and explored in immunotherapy trials. With several agents being applied in combination trials with checkpoint inhibitors (1), it is important to consider these notable differences in the TME interaction as we interpret trial results and conduct correlative analyses. Further, better understanding of the TME may inform the development of non-TKI based combinations. In our study, patients in the Angiohi/loMacrophagehi groups demonstrated the highest CSF1R expression by RNA sequencing (RNA-seq; Kruskal–Wallis test, P = 2.02E−16). CSF1R inhibitors could be a potential therapeutic avenue to explore in patients with high macrophage infiltration. However, we recognize that expression of CSF1R does not guarantee response to these inhibitors and that there are other potential targets such as FLT3 and c-MYC that may be responsible for some of the changes seen in the tumor immune microenvironment and need to be explored further (25).

Our study is not without limitations. First, our Beuselinck and Memorial Sloan Kettering Cancer Center (MSKCC) validation cohorts were small; however, despite these small numbers, we still see an initial signal consistent with the major findings of our study. Furthermore, flow-cytometry studies can be performed only prospectively. We recognize that these findings are based on a detailed characterization of a single region of the primary tumor, and ccRCC demonstrates significant intratumor heterogeneity (36). Serie and colleagues (37) recently demonstrated, among 111 patients, that ccA/ccB molecular subtypes differed between the primary and matched metastatic tumors in 43% of cases. Therefore, future studies need to assess regional differences in angiogenesis and macrophage infiltration in primary tumors and between matched primary–metastatic tumor pairs. Finally, because both arms in the COMPARZ analysis received anti-VEGF therapy, we cannot definitively conclude that our TME subgroups represent the optimal predictive biomarkers.

In conclusion, our angiogenesis and macrophage infiltration grouping demonstrates prognostic value among patients treated with TKI therapy. We speculate that patients in the AngiohiMacrophagelo group may benefit most from pazopanib therapy alone, whereas the AngiohiMacrophagehi/lo group may benefit from a combination of sunitinib- and macrophage-directed immunotherapy (e.g., CSF1R inhibition). Furthermore, the AngioloMacrophagehi group may benefit from targeted CSF1R inhibition (22). In our study, the favorable risk group was most enriched in cluster 3 (33%), which has the highest angiogenesis and lowest macrophage infiltration, and least enriched with cluster 4 (25%). These findings are supported by the recently presented results comparing ipilimumab plus nivolumab versus sunitinib (CheckMate 214), which demonstrated improved outcomes among the favorable IMDC risk group with sunitinib (38). Our findings highlight the importance of angiogenesis and macrophage infiltration in response to TKI therapy, which should be exploited further to personalize treatments, develop novel therapies, and improve outcomes among patients with metastatic RCC.

Clinical Cohorts

COMPARZ Cohort.

We analyzed data from patients who were enrolled in the COMPARZ clinical trial, which was conducted in accordance with the Declaration of Helsinki. All patients provided written informed consent for participation in the clinical trial (7). Formalin-fixed paraffin-embedded (FFPE) tumor blocks were collected at baseline from 486 patients, of whom 453 patients had sufficient tissue available and provided consent for analysis: 221 in the pazopanib arm and 232 in the sunitinib arm. Response to treatment was assessed using RECIST 1.0 (10).

ccRCC (KIRC) TCGA Validation Cohort.

Data were abstracted from 500 patients with ccRCC profiled by TCGA (4). Of these, 424 (84.8%) had clinical data including CSS, 465 (93%) RNA-seq, and 499 (99.8%) mutational data for VHL, PBRM1, BAP1, KDM5C, and SETD2. Systemic therapy treatment and response were not reliably documented and therefore not used in our analyses. RNA-seq and mutation data were downloaded from the NIH Genomic Data Commons (https://gdc.cancer.gov/).

Beuselinck Validation Cohort.

This cohort was composed of 121 patients undergoing nephrectomy from 1994 to 2011 and were subsequently treated with standard doses of sunitinib after the development of metastatic ccRCC (6). Of the 121 patients, 29 (24%) had prior cytokine immunotherapy (IL2 or interferon). Overall, 53 patients representing 53 tumor tissue samples had transcriptomic data (HuGene 1.0ST Affymetrix array) available for review. Gene-expression data (E-MTAB-3267) was downloaded from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/). Response to treatment was assessed using RECIST 1.0 (10).

MSKCC Validation Cohort.

Tumor and adjacent normal tissue were collected from 49 patients with ccRCC at the time of nephrectomy for flow cytometry and mirror FFPE blocks made for immunofluorescence. Twelve (24.5%) patients went on to receive TKI therapy (7 sunitinib/5 pazopanib). TKI response was determined by the time to treatment failure (TTF), defined as progression of disease, recurrence death from disease, or start of second-line therapy. Those with a TTF of ≤6 months were considered nonresponders and those with a TTF >6 months were considered responders.

IHC for Clinical Samples (COMPARZ)

PD-L1 expression quantification in the FFPE samples from the COMPARZ trial was performed using IHC per previously validated protocols and has been described in detail by Choueiri and colleagues (39). Briefly, patients were categorized as PD-L1/B7H1-positive when any tumor cell positivity was detected [H-score (HS) > 0]. In all cases showing any possible staining for PD-L1, a dual-color PD-L1/CD68+ stain was performed on adjacent sections to differentiate PD-L1 expression by tumor cells from that by TAM. The number of TAMs expressing PD-L1 was noted separately and semiquantitatively graded as absent, rare, moderate, or numerous. The TAM PD-L1 staining was not included in the final PD-L1 HS. Patients were categorized as TAM PD-L1–positive when patients had absent, rare, moderate, or numerous TAM PD-L1 staining.

Next-Generation Sequencing (COMPARZ)

All patients had previously given consent for sequencing purposes. DNA from the primary tumors and matched normal tissue was extracted using DNA easy kit (Qiagen) according to a standard protocol and subjected to analysis. Germline mutations were ruled out by analysis of adjacent nontumoral tissue or normal germline for every sample. A minimum of 40 ng of DNA was required for IMPACT sequencing. Samples from 377 (92%) patients were extracted with adequate DNA yields and sequenced using the MSK-IMPACT assay (40), a hybridization capture–based next-generation sequencing assay for targeted deep sequencing (approximately 500×) of all exons and selected introns of 410 oncogenes, tumor suppressor genes, and members of pathways deemed actionable by targeted therapies. Details on the MSK-IMPACT panel of targeted genes are provided in Supplementary Table S3.

Microarray Analyses (COMPARZ)

Four hundred nine patients had microarray and clinical data available for analyses. RNA from the primary tumors was extracted by AltheaDx in 2013 using the Qiagen RNAeasy FFPE kit with a modified deparaffinization step. Gene-expression profiles were derived via Affymetrix GeneChip HTA 2.0 (Affymetrix).

Computational Analysis of Microarray and RNA-seq Data

For the KIRC TCGA cohort only, RNA-seq gene level count values were computed by the summarizeOverlaps function from the R package “GenomicAlignments” (41) with UCSC KnownGene (42) in hg19 as the base gene model. The Union counting mode was used and only mapped paired reads were considered. FPKM (fragments per kilobase million) values were then computed from gene-level counts by using fpkm function from the R package “DESeq2” (43).

RNA microarray data from the COMPARZ trial were normalized to log2 value. Probes without corresponding gene symbol found were excluded from further analysis. For genes matched with multiple probes, the probe with maximum median absolute deviation (MAD) was chosen for representing the expression of the gene. The log2 normalized expression values were used for cNMF clustering and immune deconvolution analyses. For DEG and individual gene analyses, the quantile normalization values were used. QC of microarray did not demonstrate any clusters or similarity.

Gene Signatures.

A previously published angiogenesis signature (44) was utilized to measure an overall angiogenesis score. Two distinct computational methods, ssGSEA (19) and CiberSort (21), were chosen for immune deconvolution analyses. ssGSEA takes the sample gene-expression values as the input and computes an overexpression measure for the given gene list of immune cell type relative to all other genes in the transcriptome. On the other hand, CiberSort also takes gene-expression values as the input but uses a gene-expression signature matrix of particular immune cell types instead to compute the infiltration level of each immune cell type. The LM22 immune cell signature, which is validated and published along with CiberSort, was used. Marker genes of immune cell types for ssGSEA were obtained from Bindea and colleagues (45) and Senbabaoglu and colleagues (46). Infiltration levels for different immune cell types and angiogenesis scores were quantified using the ssGSEA implementation R package “gsva” (19, 47). ESTIMATE was used to calculate an immune score, which is the estimate of immune cells in tumor tissue and is calculated through the “estimate” R package (18) based on given gene-expression profile in FPKM or normalized log2 transformed values.

Clustering Analysis.

cNMF analysis was performed through the R “NMF” package (version 0.20.6; ref. 48). The top 1,500 genes with highest MAD were included for cNMF analysis. The number of clusters was evaluated between 2 and 6, with each being replicated 30 times for consensus clustering. After considering the Cophenetic coefficient, silhouette distribution, residual sum of squares, and reported clustering in the TCGA KIRC (4) and Beuselinck (6) studies, a k = 4 clustering was determined.

Analysis of DEG.

The R package “limma” (version 3.29.0) was used for DEG (49). Limma powers differential expression analyses for RNA-seq and microarray studies. Limma returned empirical Bayes moderated-t P values and adjusted P values (Q-value) to correct for multiple comparisons testing using the Benjamini–Hochberg method to control the false discovery rate (FDR). Genes with an FDR less than 1% (Q < 0.01) and fold change greater than 20% were subjected to IPA analysis.

IPA.

IPA (Ingenuity) was used for running Canonical Pathway and Diseases and Biofunctions analyses over the genes differentially expressed between clusters 4 and clusters 1 to 3.

GSEA.

DEG analysis results were used in GSEA (50) analyses against the MSigDB Hallmark gene sets (http://software.broadinstitute.org/gsea/msigdb).

Flow Cytometry for Clinical Samples (MSKCC)

ccRCC patient tissue specimens were prepared by mechanical disruption using a razor blade followed by treatment with 280 U/mL Collagenase Type 3 (Worthington Biochemical) and 4 μg/mL DNase I (Sigma) at 37°C for 1 hour with periodic vortexing. Digested tissues were passed through 70-μm filters. Resulting cells were resuspended in 44% Percoll–66% Percoll gradient (Sigma) and centrifuged for 30 minutes at 1,900 × g with no brake. Mononuclear cells were collected and immediately stained for flow cytometry analysis following Fc blocking (BioLegend) and live/dead staining (TONBO Biosciences). The antibodies used for flow cytometry were as follows: HLA-DR (L243, BioLegend #307618), CD14 (HCD14; BioLegend #325608), CD45 (2D1; eBioscience 11-9459-42), CD16 (3G8; BioLegend 302008), CD56 (HCD56; BioLegend 318318), and CD3 (7D6; Invitrogen MHCD0317). TAMs were identified as CD45+CD3LinHLA-DR+CD14+CD16+. Data are reported as a percentage of all cells.

Immunofluorescence for Clinical Samples (MSKCC)

Unstained pathologic slides of 6 renal tumors from previously untreated patients who underwent either radical or partial nephrectomy for sporadic, resectable ccRCC were obtained and reviewed by a genitourinary pathologist prior to immunofluorescence (IF). Paraffin-embedded tissue sections underwent heat-induced antigen retrieval in 10 µmol/L sodium citrate buffer. Standard IF was performed with primary antibodies against human CD68 (Thermo Scientific; catalog #ms-397-p) and CD45 (CST; clone #D9M8I). Species-specific secondary antibodies (all Abcam) were used. A Leica upright confocal microscope was used to capture images. Due to the small number of cases analyzed, we report the overall pattern of staining of tumors from nonresponders (n = 5) compared with responders (n = 1), and in areas with the highest CD45+ cell infiltrates, we reported a high or low density of CD68+ macrophages.

Statistical Analysis

Differences in patient characteristics and genomic status were tested using Fisher exact test and χ2 test (SAS 9.4, SAS Institute Inc.) for categorical variables and the Wilcoxon rank sum test [R package “ggpubr” version 0.1.5, Alboukadel Kassambara (2017). ggpubr: “ggplot2” Based Publication Ready Plots. R package version 0.1.5. https://CRAN.R-project.org/package=ggpubr] for continuous variables between sample groups. The Kruskal–Wallis test (R package “ggpubr” & SAS 9.4) was used for comparisons between three or more groups of continuous variables. Survival time was calculated from the first TKI treatment. Patients who were lost to follow-up or alive at the time of the study were treated as censored events. Survival curves were calculated according to the Kaplan–Meier method, and differences between curves were assessed using the log-rank test. The hazard ratio (HR) estimates and 95% confidence intervals (CI) were determined by the Cox proportional hazards regression modeling (SAS 9.4). To find clinical criteria related to PFS or OS, univariate modeling was performed on all the pathologic and clinical covariates. Covariates showing a significant association to prognosis (log-rank P < 0.10) at the univariate level were selected to be analyzed in multivariate models, after the exclusion of redundant covariates (see multivariate analysis; SAS 9.4, SAS Institute Inc.). All statistical analyses were post hoc; statistical significance was set at P < 0.05.

Multivariate Analysis

OS and PFS were calculated from the date of start of TKI therapy to the date of death or last follow-up up to 2 years. Patients alive at 2 years were censored with OS and PFS at 2 years. Variables that demonstrated prognostic value (P < 0.1) were carried forward. Clinicopathologic and genomic factors identified were then introduced into the multivariate analysis via backward stepwise Cox proportional hazard model (Supplementary Table S2). The discrimination of the nomogram for 2-year OS and PFS was measured by the concordance index (C-index) and observed survival probabilities respectively. Intermediate and poor risk IMDC groups were compared with the favorable group. For macrophage and angiogenesis gene-expression scores, a z-score was calculated using log-transformed raw read counts. Z-score was analyzed as a continuous variable (51). Mutational status for PBRM1, BAP1, and TP53 was incorporated with the wild-type as the reference. C-index was calculated for the baseline clinical (IMDC alone), genomic (genomic), and clinical plus genomic (IMDC + genomic) models.

Abbreviations

Manuscript abbreviations are summarized in Supplementary Table S4.

Data Availability

Access to the Beuselinck and KIRC TCGA cohort clinicopathologic, genomic, and transcriptome data is detailed in our methods section. Trial data sets for this trial will be available at https://www.clinicalstudydatarequest.com/Default.aspx after the medicine and indication have received EU and US regulatory approval on or after January 2014 and after the manuscript describing the trial results has been accepted for publication or published and no later than 18 months after the Study Completion Date (estimated June 30, 2018).

M.H. Voss reports receiving commercial research grants from BMS and Genentech, is a consultant/advisory board member for Alexion Pharmaceuticals, Calithera Biosciences, Corvus Pharmaceuticals, Exelixis, Eisai, Natera, Novartis, and Pfizer, and has received other remuneration from Eisai, Takeda, and Novartis. N. Riaz reports receiving commercial research grants from BMS and Pfizer and honoraria from the speakers bureau of Illumina. Y. Cheng owns stock in Novartis Pharm Inc. M.O. Li is a consultant/advisory board member for Leap Therapeutics. T.A. Chan reports receiving commercial research grants from Illumina, BMS, Eisai, AstraZeneca, and An2H, has ownership interest (including stock, patents, etc.) in Gritstone Oncology, and is a consultant/advisory board member for Gritstone Oncology, Illumina, and Bristol Myers. R.J. Motzer reports receiving commercial research grants from Bristol-Myers Squibb, Pfizer, Novartis, Genentech/Roche, and Eisai and is a consultant/advisory board member for Genentech, Pfizer, Merck, Novartis, and Incyte. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.A. Hakimi, M.H. Voss, F. Kuo, P. Patel, A. Reising, M.O. Li, T.A. Chan, R.J. Motzer

Development of methodology: F. Kuo, M. Liu, Y.-B. Chen, P. Patel, T.A. Chan

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.H. Voss, A. Sanchez, M. Liu, B.G. Nixon, L. Vuong, Y.-B. Chen, V. Reuter, P. Patel, R.J. Motzer

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A.A. Hakimi, M.H. Voss, F. Kuo, A. Sanchez, M. Liu, B.G. Nixon, L. Vuong, I. Ostrovnaya, N. Riaz, Y. Cheng, M. Marker, A. Reising, T.A. Chan, R.J. Motzer

Writing, review, and/or revision of the manuscript: A.A. Hakimi, M.H. Voss, F. Kuo, A. Sanchez, L. Vuong, Y.-B. Chen, V. Reuter, N. Riaz, A. Reising, T.A. Chan, R.J. Motzer

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): F. Kuo, P. Patel, T.A. Chan

Study supervision: A.A. Hakimi, P. Patel, M.O. Li, T.A. Chan, R.J. Motzer

This study was supported by the Ruth L. Kirschstein Research Service Award T32CA082088 (A. Sanchez), Weiss Family Kidney Research Fund (A.A. Hakimi), Novartis (R.J. Motzer and M.H. Voss), the Sidney Kimmel Center for Prostate and Urologic Cancers and the NIH/National Cancer Institute to Memorial Sloan Kettering Cancer Center through the Cancer Center Support Grant, award number P30 CA008748 (A.A. Hakimi, R.J. Motzer, and M.H. Voss). A.A. Hakimi also acknowledges the Parker Institute for Cancer Immunotherapy.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Choueiri
TK
,
Motzer
RJ
. 
Systemic therapy for metastatic renal-cell carcinoma
.
N Engl J Med
2017
;
376
:
354
66
.
2.
Gore
ME
,
Larkin
JM
. 
Challenges and opportunities for converting renal cell carcinoma into a chronic disease with targeted therapies
.
Br J Cancer
2011
;
104
:
399
406
.
3.
Brugarolas
J
. 
Molecular genetics of clear-cell renal cell carcinoma
.
J Clin Oncol
2014
;
32
:
1968
76
.
4.
Comprehensive molecular characterization of clear cell renal cell carcinoma
.
Nature
2013
;
499
:
43
9
.
5.
Hsieh
JJ
,
Chen
D
,
Wang
PI
,
Marker
M
,
Redzematovic
A
,
Chen
YB
, et al
Genomic biomarkers of a randomized trial comparing first-line everolimus and sunitinib in patients with metastatic renal cell carcinoma
.
Eur Urol
2017
;
71
:
405
14
.
6.
Beuselinck
B
,
Job
S
,
Becht
E
,
Karadimou
A
,
Verkarre
V
,
Couchy
G
, et al
Molecular subtypes of clear cell renal cell carcinoma are associated with sunitinib response in the metastatic setting
.
Clin Cancer Res
2015
;
21
:
1329
39
.
7.
Motzer
RJ
,
Hutson
TE
,
Cella
D
,
Reeves
J
,
Hawkins
R
,
Guo
J
, et al
Pazopanib versus sunitinib in metastatic renal-cell carcinoma
.
N Engl J Med
2013
;
369
:
722
31
.
8.
Ruiz-Morales
JM
,
Swierkowski
M
,
Wells
JC
,
Fraccon
AP
,
Pasini
F
,
Donskov
F
, et al
First-line sunitinib versus pazopanib in metastatic renal cell carcinoma: results from the international metastatic renal cell carcinoma database consortium
.
Eur J Cancer
2016
;
65
:
102
8
.
9.
Heng
DY
,
Xie
W
,
Regan
MM
,
Warren
MA
,
Golshayan
AR
,
Sahi
C
, et al
Prognostic factors for overall survival in patients with metastatic renal cell carcinoma treated with vascular endothelial growth factor-targeted agents: results from a large, multicenter study
.
J Clin Oncol
2009
;
27
:
5794
9
.
10.
Therasse
P
,
Arbuck
SG
,
Eisenhauer
EA
,
Wanders
J
,
Kaplan
RS
,
Rubinstein
L
, et al
New guidelines to evaluate the response to treatment in solid tumors. European Organization for Research and Treatment of Cancer, National Cancer Institute of the United States, National Cancer Institute of Canada
.
J Natl Cancer Inst
2000
;
92
:
205
16
.
11.
Brooks
SA
,
Brannon
AR
,
Parker
JS
,
Fisher
JC
,
Sen
O
,
Kattan
MW
, et al
ClearCode34: a prognostic risk predictor for localized clear cell renal cell carcinoma
.
Eur Urol
2014
;
66
:
77
84
.
12.
de Velasco
G
,
Culhane
AC
,
Fay
AP
,
Hakimi
AA
,
Voss
MH
,
Tannir
NM
, et al
Molecular subtypes improve prognostic value of international metastatic renal cell carcinoma database consortium prognostic model
.
Oncologist
2017
;
22
:
286
92
.
13.
Beuselinck
B
,
Verbiest
A
,
Couchy
G
,
Job
S
,
de Reynies
A
,
Meiller
C
, et al
Pro-angiogenic gene expression is associated with better outcome on sunitinib in metastatic clear-cell renal cell carcinoma
.
Acta Oncol
2018
;
57
:
498
508
.
14.
Carlo
MI
,
Manley
B
. 
Genomic alterations and outcomes with VEGF-Targeted therapy in patients with clear cell renal cell carcinoma
.
Kidney Cancer
2017
:
49
56
.
15.
Nargund
AM
,
Pham
CG
,
Dong
Y
,
Wang
PI
,
Osmangeyoglu
HU
,
Xie
Y
, et al
The SWI/SNF protein PBRM1 restrains VHL-loss-driven clear cell renal cell carcinoma
.
Cell Rep
2017
;
18
:
2893
906
.
16.
Gao
W
,
Li
W
,
Xiao
T
,
Liu
XS
,
Kaelin
WG
 Jr
. 
Inactivation of the PBRM1 tumor suppressor gene amplifies the HIF-response in VHL-/- clear cell renal carcinoma
.
PNAS
2017
;
114
:
1027
32
.
17.
McDermott
DF
,
Huseni
MA
,
Atkins
MB
,
Motzer
RJ
,
Rini
BI
,
Escudier
B
, et al
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
.
18.
Yoshihara
K
,
Shahmoradgoli
M
,
Martinez
E
,
Vegesna
R
,
Kim
H
,
Torres-Garcia
W
, et al
Inferring tumour purity and stromal and immune cell admixture from expression data
.
Nat Commun
2013
;
4
:
2612
.
19.
Barbie
DA
,
Tamayo
P
,
Boehm
JS
,
Kim
SY
,
Moody
SE
,
Dunn
IF
, et al
Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1
.
Nature
2009
;
462
:
108
12
.
20.
Franklin
RA
,
Li
MO
. 
Ontogeny of tumor-associated macrophages and its implication in cancer regulation
.
Trends Cancer
2016
;
2
:
20
34
.
21.
Newman
AM
,
Liu
CL
,
Green
MR
,
Gentles
AJ
,
Feng
W
,
Xu
Y
, et al
Robust enumeration of cell subsets from tissue expression profiles
.
Nat Methods
2015
;
12
:
453
7
.
22.
Cannarile
MA
,
Weisser
M
,
Jacob
W
,
Jegg
AM
,
Ries
CH
,
Ruttinger
D
. 
Colony-stimulating factor 1 receptor (CSF1R) inhibitors in cancer therapy
.
J Immunother Cancer
2017
;
5
:
53
.
23.
Kumar
R
,
Crouthamel
MC
,
Rominger
DH
,
Gontarek
RR
,
Tummino
PJ
,
Levin
RA
, et al
Myelosuppression and kinase selectivity of multikinase angiogenesis inhibitors
.
Br J Cancer
2009
;
101
:
1717
23
.
24.
Papaetis
GS
,
Syrigos
KN
. 
Sunitinib: a multitargeted receptor tyrosine kinase inhibitor in the era of molecular cancer therapies
.
BioDrugs
2009
;
23
:
377
89
.
25.
Davis
MI
,
Hunt
JP
,
Herrgard
S
,
Ciceri
P
,
Wodicka
LM
,
Pallares
G
, et al
Comprehensive analysis of kinase inhibitor selectivity
.
Nat Biotechnol
2011
;
29
:
1046
51
.
26.
Finke
J
,
Ko
J
,
Rini
B
,
Rayman
P
,
Ireland
J
,
Cohen
P
. 
MDSC as a mechanism of tumor escape from sunitinib mediated anti-angiogenic therapy
.
Int Immunopharmacol
2011
;
11
:
856
61
.
27.
Bi
M
,
Zhao
S
,
Said
JW
,
Adeniran
AJ
,
Xie
Z
,
Nawaf
CB
, et al
Genomic characterization of sarcomatoid transformation in clear cell renal cell carcinoma
.
PNAS
2016
;
113
:
2170
5
.
28.
Joseph
RW
,
Kapur
P
,
Serie
DJ
,
Parasramka
M
,
Ho
TH
,
Cheville
JC
, et al
Clear cell renal cell carcinoma subtypes identified by BAP1 and PBRM1 expression
.
J Urol
2016
;
195
:
180
7
.
29.
Wang
T
,
Lu
R
,
Kapur
P
,
Jaiswal
BS
,
Hannan
R
,
Zhang
Z
, et al
An empirical approach leveraging tumorgrafts to dissect the tumor microenvironment in renal cell carcinoma identifies missing link to prognostic inflammatory factors
.
Cancer Discov
2018
;
8
:
1142
55
.
30.
Veglia
F
,
Perego
M
,
Gabrilovich
D
. 
Myeloid-derived suppressor cells coming of age
.
Nat Immunol
2018
;
19
:
108
19
.
31.
Aparicio
LM
,
Fernandez
IP
,
Cassinello
J
. 
Tyrosine kinase inhibitors reprogramming immunity in renal cell carcinoma: rethinking cancer immunotherapy
.
Clin Transl Oncol
2017
;
19
:
1175
82
.
32.
Toge
H
,
Inagaki
T
,
Kojimoto
Y
,
Shinka
T
,
Hara
I
. 
Angiogenesis in renal cell carcinoma: the role of tumor-associated macrophages
.
Int J Urol
2009
;
16
:
801
7
.
33.
Gordon
SR
,
Maute
RL
,
Dulken
BW
,
Hutter
G
,
George
BM
,
McCracken
MN
, et al
PD-1 expression by tumour-associated macrophages inhibits phagocytosis and tumour immunity
.
Nature
2017
;
545
:
495
9
.
34.
Guislain
A
,
Gadiot
J
,
Kaiser
A
,
Jordanova
ES
,
Broeks
A
,
Sanders
J
, et al
Sunitinib pretreatment improves tumor-infiltrating lymphocyte expansion by reduction in intratumoral content of myeloid-derived suppressor cells in human renal cell carcinoma
.
Cancer Immunol Immunother
2015
;
64
:
1241
50
.
35.
Ozao-Choy
J
,
Ma
G
,
Kao
J
,
Wang
GX
,
Meseck
M
,
Sung
M
, et al
The novel role of tyrosine kinase inhibitor in the reversal of immune suppression and modulation of tumor microenvironment for immune-based cancer therapies
.
Cancer Res
2009
;
69
:
2514
22
.
36.
Gerlinger
M
,
Rowan
AJ
,
Horswell
S
,
Larkin
J
,
Endesfelder
D
,
Gronroos
E
, et al
Intratumor heterogeneity and branched evolution revealed by multiregion sequencing
.
N Engl J Med
2012
;
366
:
883
92
.
37.
Serie
DJ
,
Joseph
RW
,
Cheville
JC
,
Ho
TH
,
Parasramka
M
,
Hilton
T
, et al
Clear cell type A and B molecular subtypes in metastatic clear cell renal cell carcinoma: tumor heterogeneity and aggressiveness
.
Eur Urol
2017
;
71
:
979
85
.
38.
Escudier
B
,
Tannir
NM
,
McDermott
DF
,
Frontera
OA
,
Melichar
B
,
Plimack
ER
, et al
CheckMate 214: Efficacy and safety of nivolumab + ipilimumab (N+I) v sunitinib (S) for treatment-naïve advanced or metastatic renal cell carcinoma (mRCC), including IMDC risk and PD-L1 expression subgroups
.
Ann Oncol
2017
;
28
:
621
2
.
39.
Choueiri
TK
,
Figueroa
DJ
,
Fay
AP
,
Signoretti
S
,
Liu
Y
,
Gagnon
R
, et al
Correlation of PD-L1 tumor expression and treatment outcomes in patients with renal cell carcinoma receiving sunitinib or pazopanib: results from COMPARZ, a randomized controlled trial
.
Clin Cancer Res
2015
;
21
:
1071
7
.
40.
Cheng
DT
,
Mitchell
TN
,
Zehir
A
,
Shah
RH
,
Benayed
R
,
Syed
A
, et al
Memorial Sloan Kettering-Integrated Mutation Profiling of Actionable Cancer Targets (MSK-IMPACT): A hybridization capture-based next-generation sequencing clinical assay for solid tumor molecular oncology
.
J Mol Diagnost
2015
;
17
:
251
64
.
41.
Lawrence
M
,
Huber
W
,
Pages
H
,
Aboyoun
P
,
Carlson
M
,
Gentleman
R
, et al
Software for computing and annotating genomic ranges
.
PLoS Comput Biol
2013
;
9
:
e1003118
.
42.
Rosenbloom
KR
,
Armstrong
J
,
Barber
GP
,
Casper
J
,
Clawson
H
,
Diekhans
M
, et al
The UCSC Genome Browser database: 2015 update
.
Nucleic Acids Res
2015
;
43
:
D670
81
.
43.
Love
MI
,
Huber
W
,
Anders
S
. 
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
44.
Masiero
M
,
Simoes
FC
,
Han
HD
,
Snell
C
,
Peterkin
T
,
Bridges
E
, et al
A core human primary tumor angiogenesis signature identifies the endothelial orphan receptor ELTD1 as a key regulator of angiogenesis
.
Cancer Cell
2013
;
24
:
229
41
.
45.
Bindea
G
,
Mlecnik
B
,
Tosolini
M
,
Kirilovsky
A
,
Waldner
M
,
Obenauf
AC
, et al
Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer
.
Immunity
2013
;
39
:
782
95
.
46.
Senbabaoglu
Y
,
Gejman
RS
,
Winer
AG
,
Liu
M
,
Van Allen
EM
,
de Velasco
G
, et al
Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures
.
Genome Biol
2016
;
17
:
231
.
47.
Hanzelmann
S
,
Castelo
R
,
Guinney
J
. 
GSVA: gene set variation analysis for microarray and RNA-seq data
.
BMC Bioinformatics
2013
;
14
:
7
.
48.
Gaujoux
R
,
Seoighe
C
. 
A flexible R package for nonnegative matrix factorization
.
BMC Bioinformatics
2010
;
11
:
367
.
49.
Ritchie
ME
,
Phipson
B
,
Wu
D
,
Hu
Y
,
Law
CW
,
Shi
W
, et al
limma powers differential expression analyses for RNA-sequencing and microarray studies
.
Nucleic Acids Res
2015
;
43
:
e47
.
50.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
PNAS
2005
;
102
:
15545
50
.
51.
Malone
BM
,
Tan
F
,
Bridges
SM
,
Peng
Z
. 
Comparison of four ChIP-Seq analytical algorithms using rice endosperm H3K27 trimethylation profiling data
.
PLoS One
2011
;
6
:
e25260
.