An Empirical Approach Leveraging Tumorgrafts to Dissect the Tumor Microenvironment in Renal Cell Carcinoma Identifies Missing Link to Prognostic Inflammatory Factors

The interplay between a tumor and its microenvironment influences a variety of processes, including invasion, metastasis, and response to therapy. However, much still remains unknown about this interplay and the triggered inflammatory response, which could involve immune cells from both the innate and adaptive immune systems. To understand the intricate interactions between a tumor and its microenvironment, it is helpful to identify and characterize each cell population present in bulk tumor. One school of in silico approaches relies on signature genes and a backend deconvolution method, such as single sample gene set enrichment analysis (ssGSEA), an algorithm that correlates particular cell types with a set of signature genes (1). Indeed, several sets of gene signatures have been produced (2–4). However, these signatures suffer from poor definition of immune-specific genes, which can be affected by the type of cancer these cells infiltrate or, when derived from laser microdissection approaches, by technically challenging sample processing (5). Newer approaches involve single-cell sequencing, but they are still quite expensive and also require extensive sample handling, which can affect gene expression (6).

Therefore, improved approaches are essential for more accurate dissection of the heterogeneous cell populations in tumors. To generate such gene signatures with empirical evidence, it is helpful to obtain the in vivo expression of immune cells infiltrating bulk tumors by dissection of bulk tumor expression profiles. However, there are two major problems associated with current tumor dissection studies. First, in addition to the immune cell and the tumor cell components, normal parenchyma may contaminate the bulk tumor specimen. Thus, earlier two-component dissection models (7, 8) may be insufficient. Second, tumor cell line cultures, which have been used in some studies as surrogates of the tumoral component (9, 10), are subject to significant phenotypic drift from the original tumors (11, 12) and are not optimal.

Here, we undertook a novel approach to dissect the tumor cells and the tumor microenvironment (TME) in renal cancer by leveraging tumorgrafts. Tumorgrafts, or patient-derived xenografts (PDX), represent original tumor cell populations transplanted directly from a patient into a host, typically a mouse. In comparison with cell line–based xenografts, tumorgrafts more faithfully reproduce the histology, gene expression, mutations, and treatment responsiveness of individual human tumors (13–15). Despite recent reports questioning the validity of tumorgrafts (16), the reported limitations affect only certain tumor types subject to genome duplication over long passages. Importantly, in the tumorgrafts, only the human tumor cells expand, and the tumor stroma is soon replaced by the host (13, 15, 17). Thus, inasmuch as human transcripts can be distinguished from the mouse transcripts (through alignment to both reference genomes), tumorgrafts can be used to specifically define the tumor cell transcriptome. By subtracting this transcriptome and the normal tissue transcriptome from the corresponding patient tumor through a three-component dissection analysis, a gene expression profile may be generated corresponding to the immune/stroma component, namely, the TME component.

We applied this approach to renal cell carcinoma (RCC). RCC includes clear-cell RCC (ccRCC), the most common subtype accounting for 70% of all RCCs; papillary RCC (pRCC); and chromophobe RCC (chRCC; ref. 18). RCC is a paradigmatic tumor type characterized by a rich inflammatory stroma (19, 20). Interestingly, several extensively validated prognostic factors in RCCs (21, 22) are also indicators of systemic inflammation, such as thrombocytosis, anemia, and neutrophilia. However, a detailed description of how these extensively validated RCC prognostic factors are linked to the biology of the tumor cells and the TME is lacking. RCC is also one of two tumors traditionally regarded as immune-responsive (23–25). Two other tumors considered immunogenic include melanoma and, recently, lung cancer. Whereas melanoma and lung cancer stand out by their high mutation burdens, RCCs exhibit a modest overall mutational load (26). However, a recent study reported that RCCs have the highest number of indel mutations on a pan-cancer basis (19), which may yield a disproportionately higher amount of neoantigens and elicit robust adaptive immune responses.

Here, we set out to evaluate the TME in RCC using DisHet, a three-component dissection algorithm we developed. Based on DisHet analyses, we obtained an empirically defined TME gene expression signature (eTME). Using eTME, we refined previously published signatures and assigned other eTME genes to specific immune cell types using the BLUEPRINT project data. Based on expression of eTME genes, we defined an inflamed subtype (eTME-IS), which is correlated with systemic manifestations of inflammation such as thrombocytosis and anemia. Patients with eTME-IS have abundant tumor infiltration of regulatory T cells, natural killer (NK) cells, T_H1 cells, neutrophils, macrophages, B cells, and CD8⁺ T cells and are enriched for clear-cell tumors with BAP1 mutations and papillary tumors of type 2 histology. Our analyses also show that tumor cells could actively modulate the immune response in the TME.

We explored whether tumorgrafts could be used to dissect the immune/stroma component's gene expression through subtraction from the corresponding patient tumor. First, we evaluated whether tumorgrafts are purely comprised of tumor cells. We leveraged the UT Southwestern (UTSW) Kidney Cancer Program (KCP) platform and identified 35 patients (29 ccRCC and 6 non-ccRCC) for whom we have carried out RNA sequencing (RNA-seq) for one bulk tumor, one tumorgraft, and one normal tissue sample (Supplementary Table S1). We profiled the expression of genes from the ESTIMATE signature (2), a signature made up of stromal and immune genes. Overall, we found that ESTIMATE gene expression was much lower in the human cell component of tumorgrafts when compared with the bulk tumor (Fig. 1A; note log scale). As a control, we evaluated three previously identified RCC tumor-specific genes, FABP7 (27), PTHLH (28), and CD70 (29). As shown in Fig. 1A, all three transcripts were expressed at high levels in both the bulk tumor and the human component of the tumorgrafts. Next, we identified 38 patients with ccRCC for whom we have carried out exome sequencing for at least one tumor, one tumorgraft, and one normal sample (27 overlapping with the 35 previously mentioned patients). We calculated the average somatic variant allele frequencies (VAF) in both tumors and tumorgraft samples and found that, as we previously showed (15), VAFs generally increased in the tumorgrafts (Fig. 1B). As only the human tumor cells carry these somatic mutations, these data are consistent with the notion that only the tumor cell population expands in mice and that the stroma is replaced by the host.

One premise of our dissection analysis is that the gene expression of tumorgrafts does not shift significantly from those of the original human tumors. We carried out RNA-seq for one patient (XP389) for whom we have obtained primary tumors and have established one tumorgraft and one tumor cell line. We calculated the proportions of genes whose expression deviation in tumorgrafts or cell culture from primary tumors are larger than a series of cutoffs. We calculated this proportion with respect to all genes or to genes whose expression in both tumorgraft and cell culture is larger than that in the primary tumors (such genes are more likely to be the tumor-specific genes). Figure 1C shows that although there are still some differences between tumorgrafts and the tumor cells in primary tumors, the phenotypic drift in tumorgrafts is smaller than that in cell culture. These data are consistent with our previous report that the majority of tumorgrafts cluster with the corresponding patient tumor in unsupervised tumor analyses (13). Ben-David and colleagues (16) recently reported significant clonal evolution in tumorgrafts over passages through analyses of copy-number variation (CNV) data in multiple tumor types. However, it is worth noting that the drift applied almost exclusively to tumors having undergone genome duplication, which is a rare event in RCC. We also show the whole-genome CNVs in the primary tumor and tumorgraft of XP166 (11th passage) and XP243 (1st passage). In most parts of the genomes of both early and late passages of tumorgrafts, the CNVs are in a direction (amplification or loss) consistent with that of the primary tumor (Fig. 1D). This suggests that the genomic instability in RCC tumorgrafts is relatively small, even in late passages. The CNVs in the tumorgraft are generally larger in amplitude than in primary tumors, which is expected, as tumorgrafts have a higher representation of human tumor cells and the CNVs are not masked by nontumor cells. In addition, one other advantage of tumorgrafts is that they allow rapid sample processing and minimize ischemia time. On the other hand, approaches such as laser capture microdissection (LCM) and single-cell analyses would allow direct measurement of cells in the stroma and tumor tissues, but are limited by preanalytic variables such as ischemia time.

Building on the previous analysis, we subtracted from the bulk tumor gene expression, tumor cell and normal parenchyma expression using tumorgraft and adjacent normal tissue RNA-seq data. To meet the analytic needs, we developed DisHet, for Dissecting Heterogeneous bulk tumors, an algorithm based on a Bayesian hierarchical model. The dissection task can be expressed as solving for Ax = b, where A is a matrix of component-specific expression, x is a vector of mixture proportions and b is a vector of bulk tissue expressions. The three components include tumor cells, immune/stroma, and normal parenchyma. Importantly, the immune/stroma component further contains various types of immune cells, activated fibroblasts, and other cell types such as endothelial cells and pericytes. The bulk tumor expression (b) and the normal and tumor components of the A matrix are observed from the bulk tumor, normal tissue, and tumorgraft RNA-seq data, respectively, and the goal is to estimate the third unknown component in A. This concept is shown in Fig. 1E. The bulk sample expression is determined by both the amount and expression of each cellular component, both of which can be estimated by DisHet. The Supplementary File explains the full technical details of DisHet.

To evaluate DisHet, we initially performed simulation analyses. As shown in Supplementary Fig. S1A, the component proportions estimated by DisHet are almost the same as the simulated true proportions. Supplementary Figure S1B shows that DisHet can accurately estimate the average expression of each gene across all patients. Finally, DisHet can accurately estimate the individual variation of each gene by calculating a Spearman correlation between the simulated expression for each gene across all patients and the estimated expression by DisHet in that gene across all individuals. We show that most of these correlations are very close to 1 (Supplementary Fig. S1C), confirming the accuracy of DisHet in capturing the individual expression variation.

We applied DisHet on the real data from the 35 patients with ccRCC for which matched bulk tumor, tumorgraft, and normal tissue RNA-seq were conducted. As a validation, pathologic analyses of sections flanking the analyzed tumor sample for tumor proportions correlate well (Spearman correlation = 0.65, P = 0.0004) with estimations of DisHet (Fig. 2A). Second, there is a strong positive correlation (Spearman correlation = 0.72, P = 2.5 × 10⁻⁵) between the predicted tumor proportions with average VAFs of somatic mutations (Fig. 2B). Finally, ssGSEA was applied to the combined ESTIMATE stromal and immune cell genes. We observed a positive correlation (correlation = 0.88, P < 1 × 10⁻¹⁰) between the ESTIMATE enrichment score and the immune/stroma component proportions predicted by DisHet (Fig. 2C).

In addition to component proportions, we also validated the dissected immune/stroma component expression. We applied ssGSEA on the dissected immune/stroma component expression with B- and T-cell signatures from the Immunome immune cell signatures (4). We found that patients with a higher level of lymphocyte infiltration, according to hematoxylin and eosin (H&E) slides reviewed by our pathologist, did in fact have higher ssGSEA scores for T-cell and B-cell genes (Fig. 2D). We also applied ssGSEA using the 16 Winslow signatures defined in breast cancer (3) with genes in the same signature being highly correlated. We showed that genes in the same signature are indeed highly correlated in the expression of the immune/stroma component (Supplementary Fig. S2A). Lastly, we performed DisHet analyses of 9 RNA-seq datasets from 4 unique patients, each with more than one set of matched bulk tumor, tumorgraft, and normal tissue samples available. We showed that these 9 samples largely clustered by each patient according to their immune/stroma component expression (Supplementary Fig. S2B). Overall, these extensive analyses attest to the accuracy of the dissected gene expression across individuals.

We compared the performance of DisHet with another popular algorithm, DeMix (8), for gene expression deconvolution in a mixture of cell types. However, DeMix considers only two components: tumor and stroma, with the stroma component being the sum of the normal parenchyma and immune/stroma components defined by DisHet. DisHet shows much better dissection accuracy than DeMix (compare Supplementary Fig. S3A and S3B with Fig. 2A and C). DisHet is also much more computationally efficient than DeMix (Supplementary Fig. S3C).

The matched trio of bulk tumor, tumorgraft, and normal tissue RNA-seq data and DisHet analyses enabled us to define immune/stroma-specific genes with empirical evidence. We created one immune/stroma-specific gene list with 2,080 genes that are expressed on average at least 3 times higher in the immune/stroma component than in the tumor, and another with 904 genes that are expressed at least 20 times higher (Supplementary Table S2). The false discovery rates (FDR) are estimated to be 4.7% for the first set and 0.0044% for the second set based on 500 permutations. We termed these immune/stroma-specific signature genes eTME, for empirically defined gene signature of the TME (immune/stroma). We investigated the relationship of eTME genes with previously published signatures. To our surprise, 65% of the Immunome signature genes (4) are not abundantly expressed in the RCC immune/stroma component even according to the less stringent cutoff of 3× (Fig. 3A). In addition, 22% of the ESTIMATE genes and 9% of the Winslow genes are also not abundantly expressed in the RCC immune/stroma component. Although we applied these gene signatures in previous sections for validating the gene expression dissected by DisHet, they are not optimal, and further refinement filtering out nonspecific genes would be fitting. Conversely, among the 904 eTME genes expressed at >20× higher levels, we identified 778 genes without a clear interpretation in the Immunome signatures. Furthermore, even after excluding genes in ESTIMATE and Winslow signatures, there were 610 immune/stroma-specific genes that were novel.

The eTME genes may hold clues to the kidney cancer inflammatory environment and may help explain why renal cancer is highly immunogenic. To investigate the identity of these genes, we carried out gene ontology analyses for the 904 eTME genes. As expected, there was a significant enrichment for immune pathways in these genes (Supplementary Table S3). Dominant enrichment of immune-related pathways is seen even in the subset of 778 and 610 eTME genes that are not captured by the Immunome alone, or Immunome + ESTIMATE + Winslow signatures (Supplementary Table S3). eTME also contains a small fraction of genes related to nonimmune pathways/cell types. For example, fibroblast-related genes such as FAP (30) and POSTN (31) are found in the set of 904 eTME genes. This is also reflected by the extracellular matrix–related gene ontology terms that are weakly enriched in eTME (Supplementary Table S3). To derive a more specific interpretation for each gene, we evaluated the gene expression data of hematopoietic lineages (BLUEPRINT, Expression Atlas; http://www.ebi.ac.uk/gxa/home) and called signature genes by requiring their expression to be at least 2 times larger than any of the other 33 lineages, and conducting a Kruskal–Wallis test on signature gene sets requiring P < 0.001. For these analyses, we used the eTME genes with 3-fold cutoff. In this way, we assigned cell type identities to 293 genes, including 231 genes missed by the previous signatures (Supplementary Table S2). Despite the use of multiple signatures, more than 50% of the eTME genes remain unassigned. Although these genes may be expressed in more than one lineage, they may also carry information about novel cell types or differentiation states.

To validate the eTME gene signatures, we carried out single-cell RNA-seq for five patients with ccRCC without treatments prior to tumor resection. The lymphoid and myeloid cells were separately enriched from these 5 patients by FACS (using CD45 positivity and separated based on scatter properties), generating two pools of cells (n = 482 for lymphoid and n = 487 for myeloid, all patients combined), which were sequenced separately. We first aggregated single-cell gene expression data and normalized each gene by expression of the same gene in the tumor component, to derive a ratio similar to that in Fig. 3A. We plotted these ratios (logged) for “Non-eTME” genes—genes not belonging to the eTME gene signature (3× cutoff), and eTME genes defined with 3× and 20× fold cutoffs. In Fig. 3B, we show that these single-cell/tumor expression ratios for the eTME genes are mostly positive, proving the overall concordance between single-cell data and eTME. The P values of two-way t tests comparing non-eTME with eTME (3×) genes, and non-eTME with e-TME (20×) genes are smaller than 1 × 10⁻¹⁵. We also asked a complementary question, by examining genes that are most highly expressed in the single cells. By limiting to genes whose log_{(single cell/tumor)} ratios are larger than 3.5, 4, and 4.5, 71%, 86%, and 98% of the genes can be recapitulated by eTME (3×).

Next, we determined the immune cell type of each sequenced cell according to expression of eTME signature genes. Figure 3C shows the number of cells assigned to each cell type and the proportions of these cells in the myeloid pool and lymphoid pool. Reassuringly, most designated neutrophils, macrophages, and eosinophils were found in the myeloid pool and most designated T cells, NK cells, and B cells were found in the lymphoid pool. Moreover, most lymphoid cells were designated as T cells, which is consistent with previous reports of counts of lymphoid cells in RCC (32). Admittedly, technical issues may exist that led to some apparent artifacts, such as assignment of some B cells in the myeloid pool. To further validate eTME, we profiled the infiltration levels of several immune cell types in a small set of patient RNA-seq data using ssGSEA on the refined gene signatures that passed the eTME criterion of 3×. For the same patients, we carried out immunohistochemistry (IHC) for the cell-type surface markers. On average, the correlation between predicted infiltration abundances and IHC was 0.69 (Fig. 3D). On balance, these results confirm the accuracy of the eTME gene signatures that we defined with empirical evidence. We have also incorporated the eTME signatures in our DisHet package and have provided them publicly.

The identification of a highly specific RCC eTME gene signature enabled us to evaluate the functional importance of the tumor immune infiltrates. We deployed the eTME genes to assess the largest publicly available cohort of RCC samples with RNA-seq data, 884 RCCs from The Cancer Genome Atlas (TCGA). This cohort includes tumors of different histologies. We conducted principal component analysis (PCA) on these tumors based on eTME gene expression with the more stringent cutoff of 20× and used tumor-specific genes as a control. Notably, the eTME gene expression signature was able to cluster RCCs of the same histology with good accuracy (Supplementary Fig. S4). We contrasted the eTME-based clustering with clustering based on tumor-specific genes and observed a similar clustering structure based on tumor-specific gene expression. However, the clustering was less robust in separating out the pRCC and chRCC histologies than eTME genes. These results show that the TME carries important discriminating features corresponding to the different histologies.

Building on the previous observations, we took a more focused approach to elucidate the heterogeneity of eTME expression profiles in each patient. We focused on ccRCC first and clustered all TCGA patients with ccRCC (n = 529) based on expression of the eTME genes (20×). We observed two clusters of patients, one accounting for about 42% of patients and the other accounting for about 58% of patients (Fig. 4A). To characterize each cluster, we evaluated the infiltrating levels of various types of immune cells in each group of patients by applying the ssGSEA method on cell type–specific signature genes that passed the more relaxed eTME criterion of 3×. We found that the first cluster of patients was enriched for Tregs, NK cells, T_h cells, neutrophils, macrophages, eosinophils, B cells, CD8⁺ T cells, and aDC cells. This cluster was also characterized by high expression levels of the complement system proteins (represented by the C1q signature from the Winslow signatures). The second cluster was, in contrast, enriched for angiogenesis, pDC, and mast cells. Thus, the eTME signature separated tumors into an IS that included an extensive immune infiltrate and a second subtype that was not inflamed (NIS). We found that IS was enriched for mutations in BAP1 (P = 7.7 × 10⁻⁵), a tumor suppressor associated with high-grade cancers (15, 33). In contrast, PBRM1 mutations were more homogeneously distributed across both subtypes, though enriched in NIS.

We extended our observations to two other, albeit significantly smaller, ccRCC cohorts. We trained a support vector machine (SVM) classifier based on the unsupervised clusters and predicted the subtype identities of patients with ccRCC and non–clear-cell RCC (nccRCC) from the UTSW KCP database (181 patients in total, 39% ccRCCs, 24% pRCCs, 15% chRCCs, 22% other) with bulk-tumor RNA-seq data available (Fig. 4B). As in the TCGA ccRCC cohort, in the KCP cohort, eTME-IS patients with ccRCC were associated with enriched BAP1 mutations, although statistical significance was not achieved (P = 0.21, N = 70). Next, we analyzed the exome sequencing and RNA-seq data from another ccRCC cohort (34). In this cohort, we found 9 patients with BAP1 mutations out of all 95 available patients. The SVM model predicted 41 eTME-IS patients, which include 7 BAP1-mutated tumors (P = 0.032). Given the intriguing correlation of BAP1 and eTME-IS, we carried out CD4 and CD8 IHC staining in our recently established RCC mouse models, Six2-Cre;Vhl^F/F;Pbrm1^F/F and Six2-Cre;Vhl^F/F;Bap1^F/+ (35). We observed more CD4 and CD8 T-cell infiltration in Bap1-mutated mice than in Pbrm1-mutated mice, suggesting a causal relationship may exist between BAP1 mutations and the eTME-IS phenotype (Supplementary Fig. S5).

Next, we extended our analyses to nccRCC. We evaluated TCGA patients with pRCC and chRCC (Fig. 4C). The eTME-IS patients accounted for a minority of nccRCC tumors. Interestingly, the pRCC tumors in the eTME-IS group were predominantly type 2 (P = 7.2 × 10⁻⁵), which have a worse prognosis than type 1 (36). We also carried out ssGSEA analysis using the refined eTME signatures in both KCP and TCGA nccRCC cohorts (Fig. 4B and C). Across all cohorts examined, there was enrichment of Tregs, NK cells, T_H1 cells, neutrophils, macrophages, B cells, and CD8⁺ T cells, and the complement system in eTME-IS patients. Therefore, the grouping of patients with ccRCC based on eTME expression profiles revealed a pan-RCC subtype with shared similarities across different RCC histologies.

Recently, Ayers and colleagues (37) reported an 18-gene IFNγ gene signature whose activation predicts clinical response to PD-1 blockade. We analyzed this gene signature in all three patient cohorts (Fig. 4), and we discovered a very strong correlation of the IFNγ activation signature with the eTME-IS patient subtype. In synergy with this, the expression of PD-1 (PDCD1) is higher in eTME-IS patients than in eTME-NIS patients (Fig. 4A). These results reveal the potential of the eTME analysis and eTME-IS/NIS patient subtyping to be translated to clinical practice.

Next, we asked whether the differences in inflammatory phenotype correlated with tumor characteristics on imaging. We focused on the 35 patients with ccRCC analyzed by DisHet. The SVM classifier predicted that 21 of these patients belonged to eTME-IS. An experienced radiologist (I. Pedrosa) examined all MRI/CT images for these patients (32 as available) blindly and observed (Fig. 4D) a significant correlation (P = 0.034) between eTME-IS and infiltrative tumor margins (including partial infiltration). In addition, we also reviewed H&E staining images (P. Kapur), which showed that eTME-IS samples are associated with more necrosis than eTME-NIS samples (Fig. 4E), which also indicates poor prognosis (38).

To determine whether the inflamed TME is driven by the tumor cells, we asked whether human tumors, when transplanted in mice, stimulated an inflammatory response similar to that in the human stroma as determined by tumor infiltration by inflammatory murine cells. We first focused on the 35 patients with ccRCC. Because the tumorgraft hosts, NOD/SCID mice, are inherently deficient in B and T cells, we focused on neutrophils and NK cells. By quantifying the human and mouse immune cell populations in the original human tumors and the mouse-specific gene expression from tumorgrafts, respectively, we found that the mouse neutrophil and NK cell infiltrates were induced to an extent comparable in humans (Fig. 5A). Comparing eTME-IS and eTME-NIS tumors, we showed that eTME-IS tumors have a higher level of neutrophils than eTME-NIS tumors (P = 0.021), and they also elicit higher levels of neutrophils in the mouse (P = 0.031). A similar trend was also observed for NK cells (Fig. 5A, P = 0.036 in human and P = 0.23 in mice).

Although the above analyses showed a trend of concordancy, the characteristics of the immune infiltrates could be significantly different between human disease and murine tumorgraft environment relating to factors such as species and microbiota. Therefore, we used another analysis that more broadly investigates the eTME genes as a whole. We hypothesized that expression of mouse immune/stroma-specific genes in different tumorgrafts established from the same patients should be similar to each other, if tumors indeed drive stromal responses. We analyzed RNA-seq data for a set of 37 tumorgrafts from a total of 14 unique patients from our database, for whom at least two tumorgrafts had been established and analyzed by RNA-seq for each patient (Supplementary Table S1). We carried out PCA of mouse gene expression quantified from the tumorgraft RNA-seq data, for genes that belong to the innate immunity system and satisfied the eTME cutoff of 3×. As can be seen in Fig. 5B, tumorgrafts from the same patients tend to cluster into the same entities. As a negative control, we carried out this analysis for all non-eTME (3× cutoff) mouse transcripts and observed almost no separation between tumorgraft samples of the same patients (Fig. 5B). To further illustrate this analysis from another perspective, we calculated the Euclidean distances between tumorgraft pairs from the same patients as one group, and pairs from different patients as another group, determined by murine expression of eTME genes in innate immunity system, all eTME genes, and all other genes. The density plot in Fig. 5C shows that the distances between tumorgrafts from the same patients are much smaller than those from different patients, in terms of expression of eTME genes in the innate immune system (Kolmogorov–Smirnov test, P = 7.7 × 10⁻¹⁵). In contrast, interpatient and intrapatient distances are almost indistinguishable according to expression of other non-eTME genes (P = 0.98, Fig. 5C). Interestingly, echoing this active role of tumor cells on stromal immune responses, our previous work showed that RCC tumors that induced hypercalcemia in humans also induced hypercalcemia in tumorgraft mice (13).

We next asked whether there was a correlation between the inflammatory subtype (eTME-IS) we identified from expression profiles with systemic manifestations of inflammation and other systemic factors. We extracted clinical lab records from UTSW KCP patients from our database. We found that patients with eTME-IS tumors had higher platelet counts (P = 9.4 × 10⁻⁴), and lower hemoglobin levels (P = 8.1 × 10⁻¹⁰; Fig. 6A), which are established prognostic markers in RCC and may signal systemic inflammation (22). In contrast, neutrophil counts were not significantly different between these two groups of patients (P = 0.71). As a control, creatinine, a marker of renal function (39), showed no significant difference between patients with eTME-IS and eTME-NIS tumors (P = 0.20; Supplementary Fig. S6). Although not all of these parameters are available for the TCGA cohort, we asked whether a correlation existed between the inflamed TME with thrombocytosis and anemia. In all the TCGA patients with RCC combined, higher proportions of eTME-IS patients had elevated platelet counts (P = 0.039) and low hemoglobin levels (P = 0.019) compared with eTME-NIS patients (Fig. 6B). Overall, these data show that tumors with eTME-IS are associated with thrombocytosis and anemia, which may represent systemic manifestations of tumor-induced inflammation.

Terminal cancer, poorly controlled inflammatory states, and wasting syndromes are associated with malnourishment and lower protein (albumin) levels in the blood. We asked whether a correlation existed between the inflammatory microenvironment and low albumin. Interestingly, eTME-IS patients had lower albumin levels than eTME-NIS patients (P = 2.8 × 10⁻⁷; Supplementary Fig. S6). Another indicator of aggressive cancer (40), whose pathogenesis is not always understood, is low circulating sodium. Our results show that eTME-IS patients seem to have lower sodium levels as well (P = 0.014; Supplementary Fig. S6).

Given our findings that the inflamed TME correlates with systemic markers of inflammation, which are validated prognostic factors in RCC (21, 22), we asked whether the inflamed TME correlated with patient survival. Our Kaplan–Meier analyses showed that eTME-IS patients had significantly worse overall survival across all 3 patient cohorts analyzed [Fig. 6C; HR = 1.8 (95% CI = 1.3–2.5) and P = 1.4 × 10⁻⁴ in the TCGA ccRCC cohort; HR = 11.9 (4.1–34.5) and P = 4.6 × 10⁻⁶ in the UTSW cohort; and HR = 2.3 (1.2–4.5) and P = 0.014 in the TCGA nccRCC cohort].

There are several unique advantages of our methodologies over previous efforts to dissect the TME. In our study, we took the novel approach of including matched tumors, tumorgrafts, and normal tissue samples from the same patients to dissect out the different components of the TME. The tumorgrafts provided more accurate transcriptomic profiles of tumor cells grown in vivo than cultured tumor cells and allow the accurate subtraction of the tumor cell component from the tumor bulk. Second, this analysis involving matched trio samples is empowered by DisHet, a novel and efficient Bayesian algorithm. DisHet properly considers the possible existence of all three cellular components in a bulk tumor.

With the DisHet analyses of the matched trio RNA-seq data, we refined, with accurate empirical evidence, previous immune/stroma-specific genes such as Immunome, ESTIMATE, and Winslow and in addition annotated many other immune/stroma-specific genes not previously reported. Our analyses enabled the generation of an empirical TME gene signature, which we refer to as eTME. Importantly, immune-cell quantification methods such as ssGSEA rely on accurate gene signatures. The eTME genes defined by the DisHet analysis could be used to refine these signatures, as is done in our study, and to make the immune cell quantification much more relevant. Moreover, the TME of each cancer type could be very different. Our eTME gene signatures are defined with RCC data and should have very high accuracy in RCC, whereas ESTIMATE identified leukocytes in breast, colorectal, and ovarian cancers. Ideally, such tumor-specific gene signatures should be established for each type of cancer. The NCI Patient-Derived Models Repository (PDMR) repository was recently released and has provided RNA-seq data for more than 600 tumorgraft models for 24 cancer types. Valuable resources such as PDMR can be utilized for achieving this goal.

The advantages of our approach complement insights that can be provided through other techniques such as single-cell sequencing and LCM. It will be more powerful if these approaches and tumorgraft models can be utilized in a synergistic manner to tackle tumor heterogeneity and tumor immunity, as is also exemplified by some of the analyses carried out in this study. Over time, the cost and complexity of single-cell analysis and LCM are expected to decrease, which should make these approaches more affordable and feasible.

Recently, Ben-David and colleagues (16) reported that tumorgrafts do not represent an exact replica of the original tumors in multiple tumor types. Notably, their study shows that RCCs are the second-most stable tumor next to brain tumors. Indeed, our own work (13, 15) has shown that passage has insignificant impact on CNV and gene expression, even in late passage. Although in this work all tumorgrafts except for one were from early passage (≤3), one was from passage 11, and we did not find substantial differences compared with earlier passages. In addition, more generally, the genomic instability reported by Ben-David and colleagues is inferred from gene expression data, rather than directly measured at the DNA level, and focuses on arm-level CNV. Furthermore, as reported by the authors, their findings are largely related to tumors that experience genome duplication, which is rare in RCC.

The eTME signature we identified was applicable across RCC histologies, and we show that it offers greater discriminatory power of histologic subtypes than a signature focusing on tumor-specific genes. These data suggest the existence of unique TME features that accompany the different histologic subtypes. Furthermore, our work shows that eTME gene expression can separate tumors into a highly inflammatory subtype (eTME-IS) and a noninflammatory subtype (eTME-NIS) with a characteristic immune cell infiltration repertoire and stromal pathway activation. This clustering of eTME-IS and eTME-NIS patients and their associated immune infiltration patterns seems to be corroborated by flow cytometric findings by Giraldo and colleagues (41). In their study, the authors identified a ccRCC “immune-regulated” subtype of patients, characterized by inflammation, CD4⁺ICOS⁺ cells with a Treg phenotype, and high risk of disease progression, which likely corresponds to the eTME-IS phenotype. On the other hand, we found higher numbers of mast cells in eTME-NIS tumors, which are associated with a better outcome in keeping with a recent report showing that tumor-infiltrating mast cells are associated with a survival advantage in localized ccRCC (42).

The TME may have either a protumor or antitumor effect. Interestingly, eTME-IS tumors may be more susceptible to immunotherapy. eTME-IS tumors are characterized by an IFNγ signature, which has been associated with response to PD-1 blockade (37). It is also in alignment with the fact that patients with worse prognosis may have a higher chance of response to nivolumab plus ipilimumab, and perhaps even nivolumab alone (43, 44).

Unsupervised clustering analyses revealed that BAP1-mutant ccRCC tumors are correlated with an inflammatory TME (eTME-IS). Interestingly, BAP1 deficiency has been associated with a peritoneal inflammatory response upon exposure to asbestos fibers in mice (45) and with an inflammatory phenotype in primary uveal melanoma (46). Our CD4 and CD8 IHC experiments in the genetically engineered mouse models suggest a causal relationship between BAP1 loss and immune infiltrations in RCC, though more extensive experiments are needed for confirmation.

Our data provide insight into a missing link between tumor subtypes and important systemic prognostic factors, which may be indicators of inflammation. Specifically, anemia, thrombocytosis, and neutrophilia correlate with inflammatory conditions beyond cancer. These markers are associated, for instance, with autoimmune conditions such as rheumatoid arthritis, Crohn disease, and ulcerative colitis, as well as some infections such as granulomatous infections. However, their molecular pathogenesis, and specifically how these manifestations are linked to RCC, is poorly understood. Our study shows that the eTME-IS phenotype is associated with anemia and thrombocytosis. Inasmuch as the inflammatory TME may account for thrombocytosis or anemia, this study sets the foundation for identifying mediators of these systemic responses and, as such, begins to shed light on features in the tumor/TME that may affect prognosis.

Overall, our work, which was enabled by important methodologic advances, provides novel biological insights into the TME based on empirical analyses of renal tumors, which could eventually lead to novel therapeutic interventions. Especially in the modern era of immunotherapy, a better characterization of tumor–microenvironment interactions may hold clues to more informed treatment decisions involving immunotherapy. The mechanism of tumor immunogenicity could be vastly different in different tumor types. For example, in some other tumor types (e.g., melanoma or head and neck cancers), high lymphocyte infiltration is associated with better prognosis. However, our methodology is general and should have wide applicability for understanding the roles of immunogenicity and inflammation in development, progression, metastasis, and treatment of a large variety of cancers.

The study was performed following a protocol approved by the Institutional Review Board of the University of Texas Southwestern Medical Center that involves analyses of kidney cancer (or control) samples donated by patients. The studies were conducted in accordance with the Belmont Report and U.S. Common Rule guidelines. The investigators obtained informed written consent from the subjects regarding whether their tissue samples can be used in this study. We carried out RNA-seq for 403 samples and exome sequencing for 343 samples from 195 patients of the UTSW KCP cohort (Supplementary Table S1). RNA-seq and exome-sequencing analyses were performed from the same samples using protocols that allow simultaneous DNA and RNA extraction (47). Small samples to be extracted (∼50 mg) were flanked by sections used for pathologic imaging analyses. In addition to genomics data, we also queried clinical lab results and follow-up data for the UTSW patients. H&E images were reviewed by a pathologist (P. Kapur) blindly to report percentages of tumor nuclei. Among the 35 tumorgraft models that were used in the primary analysis (Supplementary Table S1), 25 are first passage, 7 are second passage, 2 are third passage, and one is the 11th passage. The cell line derived from the tumor of patient XP389 (Fig. 1C) was established in 2012 in the lab of Dr. James Brugarolas.

Four-micron-thick sections of formalin-fixed, paraffin-embedded tissue were subjected to heat-induced epitope retrieval using CC1 (Ventana Medical Systems, Inc.), a Tris-based buffer at pH 8–8.5, and IHC staining. We used rabbit monoclonal CD3 (Ventana Medical Systems, Inc.), mouse monoclonal CD4 (Leica Biosystems Ltd.), rabbit monoclonal CD8 (Ventana Medical Systems, Inc.), mouse monoclonal CD20 (Ventana Medical Systems, Inc.), mouse monoclonal CD68 (Dako), mouse monoclonal CD56 (Dako), and mouse-monoclonal CD117 (Dako). The slides were scanned for each case using the Aperio ScanScope (Aperio Technologies). After scanning the slides, relevant areas of tumor were selected for analysis. Aperio membrane IHC algorithm was used for each marker, which is based on the spectral differentiation between brown (positive) and blue (counter) staining. Results of Aperio IHC quantification were also confirmed manually. We calculated the proportion of positive cells out of all cells in the sampled areas. The final results were validated by another pathologist, M. Chen.

The DisHet model dissects the bulk tumor RNA-seq data such that patient-to-patient differences in the immune/stroma gene expression profiles are minimized, based on the assumption that the immune/stroma components of different patients share similarities in gene expression. At the same time, the “remaining” patient-to-patient differences can reflect the true individual variation in gene expression of the immune/stroma component. The DisHet algorithm is based on a Bayesian hierarchical model, and details of model setup and estimation procedures can be found in the Supplementary File.

RNA-seq reads that passed quality control were aligned to mouse genome (mm10) and human genome (hg38) using HISAT2. Reads that were primarily aligned to mouse (mm10) were removed. By aligning a human tumor RNA-seq dataset competitively to the combined genome of human and mouse, we estimated the possibility of misalignment of a sequencing read to the wrong genome to be only 0.047%. The remaining reads were then realigned to the human reference genome (hg38). Exon-level read counting was performed using featureCounts. Multi-mapping reads were included in the final count for each exon. For each sample, a constant value of 1 was added to the raw read counts per gene. Expression matrix was quantile normalized. For exome-sequencing data, BWA was used to map Fastq files to the combined reference genome of human (hg19) and mouse (mm10). Reads that had primary alignments in the mouse genome were dropped. Then remaining reads were mapped again to hg19 using BWA. Duplicates were removed by Picard tools and sequencing data of the same sample across lanes were merged. Realignment was carried out around indels using RealignerTargetCreator from the GATK toolkit, and BaseRecalibrator from the GATK toolkit was used to carry out base recalibration. Somatic SNP and indel mutations were called using MuTect2, Varscan, Strelka, and SpeedSeq and combined. CNV was detected from exome-sequencing data using in-house software (publication in preparation).

Resected tumors from 5 patients with ccRCC were rinsed with PBS, dissociated mechanically into tiny pieces, and transferred into a 50 mL conical tube containing 10 mL prewarmed Digest media (10 mL Processing media, 40 μL Dispase stock, 40 μL collagenase stock, and 20 μL DNAse stock, preprepared). Tumor tissues were incubated in digestion media for 10 minutes at 37°C, then vortexed for 10 seconds, and pipetted up and down for 1 minute using pipettes of descending sizes (25, 10, and 5 mL). The cell suspension was filtered using a 70-μm nylon mesh to remove residual cell clumps and centrifuged at 580 × g at 4°C for 6 minutes. CD45-positive cells of different granularity that were viable (Calcein blue high) were sorted into each well of a 96-well plate (1 plate myeloid and 1 plate lymphoid) in the UTSW Flow Cytometry Core. Reverse transcription and whole transcriptome amplification of single cells were performed with SMART-Seq v4 Ultra Low Input RNA Kit (Takara) following the manufacturer's protocol. Amplified cDNA products were then purified with 0.9X Agencourt AmPure XP-beads (Beckman). Following purification, the PCR product of each single cell was analyzed for fragment size distribution with Bioanalyzer (Agilent Technologies) using High Sensitivity DNA reagent (Agilent Technologies) and quantified with Quant-IT picogreen dsDNA assay kit (Thermo Fischer scientific). Single-cell RNA-seq libraries were then constructed using Nextera XT DNA sample preparation kit (Illumina) and Nextera XT index kit (Illumina). Libraries were pooled together, quality checked, and quantified using Bioanalyzer (Agilent Technologies) and Qubit HsDNA kit (Invitrogen). Constructed libraries were then sequenced by an Illumina Hiseq2500 or Hiseq4000 sequencer with 1 × 50 bp reads.

Single-cell RNA-seq data were processed by standard bioinformatics pipelines and assembled into a gene expression matrix, |${E_{ij}}\;( {i = 1 \ldots I,\;j = 1 \ldots J} )$|⁠, of I genes × J cells with log transformation and quantile normalization. We used the following approach to assign cell type to each sequenced cell. First, we centered each gene's expression by |$E_{ij}^c = {E_{ij}} - {\rm{median}}( {{E_{i.}}} )$|⁠. We denote the set of signature genes for each cell type, |$k = 1 \ldots K$|⁠, as |${S_k}$|⁠. The signature enrichments for this cell type across all cells, j = = 1…J, are calculated as |$e_j^k = {\rm{median}}( {E_{ij}^c\;{\rm{|}}i \in {S_k}} )$|⁠. This set of enrichment scores for each cell type is normalized by median and standard deviation: |$e_j^{n,k} = ( {e_j^k - {\rm{median}}( {e_.^k} )} )/sd( {e_.^k} )$|⁠. A sequenced cell |$j$| is assigned cell type |$k$| if and only if |$e_j^{n,k} = \max ( {e_j^{n,.}} )$| and |$e_j^{n,k} \gt {\rm{median}}( {e_.^{n,k}} )$|⁠. We successfully assigned cell-type identities to 663 cells (68.5% of all cells).

All computations and statistical analyses were carried out in the R computing environment. All hypothesis tests performed in this study are two-sided, unless otherwise specified. For comparison with DeMix, we ran the DeMix Linux version with the following parameters: Cit = 64, Matched = 1, Nstart = 2, Nref = 0, and Nhuman = 1 (8). ssGSEA analysis was performed using the R GSVA package by calling the gsva function with parameter method = “ssgsea” and rnaseq = T (48). In this study, the ESTIMATE stromal and immune signatures are combined as one unique signature. Gene ontology analysis is carried out on the GOrilla server (49). Hierarchical clustering of TCGA patients with ccRCC was based on Euclidean distance and the “complete” agglomeration method. We used an SVM to train a classifier to learn the two subtypes generated by unsupervised clustering of TCGA patients with ccRCC. We used the svm function in the e1071 R package with sigmoid kernel. For survival analyses, we used the log-rank test based on dichotomized patient sets by the eTME-IS/NIS subtypes. Patients' overall survival data were visualized by the Kaplan–Meier estimator. For all box plots appearing in this study, box boundaries represent interquartile ranges, whiskers extend to the most extreme data point which is no more than 1.5 times the interquartile range, and the line in the middle of the box represents the median. We used Spearman correlation for all correlations calculated for validation purposes in Fig. 2 and Supplementary Fig. S3.

Data are provided for UTSW patients specifically consenting to placement of their raw genomic data in a protected publicly accessible database. Sequencing data from consented patients are deposited in the European Genome–phenome Archive (EGA), with accession numbers EGAS00001002786 and EGAS00001000926 (also shown in Supplementary Table S1).

The DisHet algorithm is freely available at https://cran.r-project.org/web/packages/DisHet and http://lce.biohpc.swmed.edu/dishet. The eTME gene signatures for immune cells we defined in this study are also incorporated into both the R package and the webserver.

No potential conflicts of interest were disclosed.

Conception and design: T. Wang, Z. Zhang, M. Chen, J. Brugarolas

Development of methodology: T. Wang, J.J. Luke, X. Wang, M. Chen, A. Filatenkov, J. Torrealba, J. Brugarolas

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P. Kapur, B.S. Jaiswal, R. Hannan, Z. Zhang, I. Pedrosa, Q. Yousuf, A. Joyce, M.S. Kim, M. Chen, A. Filatenkov, J. Torrealba, E. Stawiski, Z. Modrusan, S. Durinck, S. Seshagiri, J. Brugarolas

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): T. Wang, R. Lu, P. Kapur, B.S. Jaiswal, R. Hannan, Z. Zhang, I. Pedrosa, J.J. Luke, H. Zhang, L.D. Goldstein, X. Wang, M. Chen, A. Filatenkov, X. Luo, J. He, S. Durinck, J. Brugarolas

Writing, review, and/or revision of the manuscript: T. Wang, B.S. Jaiswal, R. Hannan, Z. Zhang, I. Pedrosa, J.J. Luke, X. Wang, M. Chen, J. Torrealba, S. Seshagiri, J. Brugarolas

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Hannan, Q. Yousuf, Y.-F. Gu, T. McKenzie, M.S. Kim, D. Luo, O. Onabolu, C. Stevens, M. Chen, W. Guo, J. He, E. Stawiski, J. Brugarolas

Study supervision: T. Wang, J. Brugarolas

Other (performed management, curation, and preprocessing of patient data from the UTSW Kidney Cancer Program database): Z. Xie

We acknowledge the UT Southwestern Tissue Shared Resource for collecting patient tissue samples. We acknowledge the UTSW Flow Cytometry Core for providing help with the single-cell sequencing experiments. We acknowledge the UTSW Kidney Spore Histology core for providing help with all IHC experiments. We acknowledge Eva Moresco and Jessie Norris for their helpful comments on the writing of this article. This study was supported by the NIH (R03 ES026397-01 to T. Wang; R01 GM110050-01A1 to M. Chen; R01 CA175754 to J. Brugarolas; SPORE 1P50CA19651601 to P. Kapur, J. Brugarolas, T. Wang, and I. Pedrosa; R01 CA154475 to I. Pedrosa; CCSG 5P30CA142543 to T. Wang), Center for Translational Medicine of UT Southwestern (SPG2016-018 to T. Wang), Cancer Prevention and Research Institute of Texas (CPRITRP150596 to M.S. Kim, Z. Zhang, R. Lu, D. Luo), and American Cancer Society (RSG-16-001-01-CCE to R. Hannan). This work was also partially supported by fundraising efforts orchestrated by the KCP Patient Council and the Kidney Cancer Coalition.

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.