Transcriptomic Profiling of the Tumor Microenvironment Reveals Distinct Subgroups of Clear Cell Renal Cell Cancer: Data from a Randomized Phase III Trial

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.

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.

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 (χ² 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).

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

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 Angio^hi 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 Angio^lo 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 Angio^hi group (relative to the median) did not demonstrate improvement in either OS (P = 0.05) or PFS (P = 0.08) when compared with the Angio^lo 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.

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 (χ² 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 Angio^hi group demonstrated improved TKI response, OS, and PFS. However, the association between angiogenesis and survival is weakest in cluster 4.

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.

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 (Macrophage^hi, 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 (M2^hi, 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 M2^lo group (Supplementary Fig. S3G). This association was attenuated when we evaluated the impact of both macrophage types (M1^hiM2^hi, 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, Macrophage^hi infiltration alone is associated with worse OS and decreased objective response in TKI-treated patients.

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⁺CD3⁻CD19⁻CD14⁺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).

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 Angio^lo/Angio^hi and Macrophage^lo/Macrophage^hi. Strikingly, these groups demonstrated significant differences in median OS (P = 1.74E−5) and PFS (P = 0.0006; Fig. 4D). Patients in the Angio^loMacrophage^hi group demonstrated the worst outcomes compared with the Angio^hiMacrophage^lo 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 Angio^lo/hiMacrophage^hi groups and cluster 3 demonstrates enrichment in Angio^hiMacrophage^lo/hi groups (χ² test, P = 3.37E−44; Supplementary Fig. S4B). Similarly, the poor-risk IMDC group was also enriched with Angio^lo/hiMacrophage^hi 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 Macrophage^hi 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.

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

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 Angio^loMacrophage^hi 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 Angio^hiMacrophage^lo 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 Angio^hi 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 Macrophage^hi 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 Macrophage^lo 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.

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 Macrophage^hi 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 Angio^hi/loMacrophage^hi 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 Angio^hiMacrophage^lo group may benefit most from pazopanib therapy alone, whereas the Angio^hiMacrophage^hi/lo group may benefit from a combination of sunitinib- and macrophage-directed immunotherapy (e.g., CSF1R inhibition). Furthermore, the Angio^loMacrophage^hi 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.

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

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

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

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.

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.

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.

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

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 log₂ 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 log₂ 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.

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 log₂ transformed values.

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.

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

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

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⁺CD3⁻Lin⁻HLA-DR⁺CD14⁺CD16⁺. Data are reported as a percentage of all cells.

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.

Differences in patient characteristics and genomic status were tested using Fisher exact test and χ² 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.

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.

Manuscript abbreviations are summarized in Supplementary Table S4.

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.