Mapping single-cell transcriptomes in the intra-tumoral and associated territories of kidney cancer

Abstract

Tumor behavior is intricately dependent on the oncogenic properties of cancer cells and their multi-cellular interactions. To understand these dependencies within the wider microenvironment, we studied over 270,000 single-cell transcriptomes and 100 microdissected whole exomes from 12 patients with kidney tumors, prior to validation using spatial transcriptomics. Tissues were sampled from multiple regions of the tumor core, the tumor-normal interface, normal surrounding tissues, and peripheral blood. We find that the tissue-type location of CD8⁺ T cell clonotypes largely defines their exhaustion state with intra-tumoral spatial heterogeneity that is not well explained by somatic heterogeneity. De novo mutation calling from single-cell RNA-sequencing data allows us to broadly infer the clonality of stromal cells and lineage-trace myeloid cell development. We report six conserved meta-programs that distinguish tumor cell function, and find an epithelial-mesenchymal transition meta-program highly enriched at the tumor-normal interface that co-localizes with IL1B-expressing macrophages, offering a potential therapeutic target.

Keywords: IL1B; kidney cancer; multi-regional sequencing; pseudocapsule; renal cell carcinoma; single-cell sequencing; tumor microenvironment.

Graphical abstract

Figure 1

Sampling strategy and overall tissue distribution of the major cell types in RCC (A) Sampling strategy for each of 12 patient donors. a, c, d, and e represent four different regions of the tumor core; g, tumor-normal interface; f, perinephric fat; n, normal kidney; b, peripheral blood; h, normal adrenal gland; i, adrenal metastasis; t, thrombus. a1, a2, a3, and a4 represent LCM biopsies in tumor region a; ST, spatial transcriptomics. (B) Overall UMAP of all cells in our study. (C) Heatmap showing top differentially expressed genes (DEGs) in each of the major cell types. (D) UMAP and bar plots showing tissue distribution of the major cell types. Colors in the bar plots correspond to those in the UMAP here and in (B). (E) Dot plot showing tissue distribution of the fine-grained annotated cell types. Proportions are calculated as dividing cell numbers by total cell numbers of a certain major cell compartment. See also Figures S1 and S2; Table S2. Information on scRNA-seq and Visium spatial transcriptomics, related to Figure 1, Table S3. Information on LCM samples and WES, related to Figure 1, Table S4. Top differentially expressed genes in the subclustering analysis of different cell lineages, related to Figure 1

Figure 2

Cross-study comparisons of different cell types Cell types and annotations our study (rows) were compared with annotations or cluster numbers (columns) reported by four previous RCC studies, namely Biet al., Krishna et al., Braun et al., and Borcherding et al. The comparison was based on logistic regression trained models using CellTypist with our data as the reference. Dot size represents the fraction of cells predicted as certain cell types, and color scale represents mean probability of prediction.

Figure 3

CD8⁺ T cell characterization, clonality, exhaustion, and regional enrichment (A) Dot plot showing marker gene expression defines principal CD8⁺ cell types. EM, effector memory; Act, activated; EFF, effector; EX, exhausted. (B) UMAP depicting the pseudotime inference of CD8⁺ cells. (C) Expression of canonical exhaustion markers across cells ordered by pseudotime analysis. All marker genes are statistically significant across pseudotime values (q value = 0; Moran’s I test in Monocle 3). (D) UMAP showing the ten most expanded clones from patient PD43948. Gray dots represent cells outside the ten most expanded clones. (E) Box plot depicting the most expanded clonotypes (>100 cells) across all patients, ordered by their mean pseudotime values, showing the median, interquartile range, and outlier pseudotime values. Statistical analysis by two-sided Wilcoxon rank-sum test. Data are presented as boxes with the median ± first and third quartiles with notches depicting 95% confidence intervals (top panel); bar plots showing the maximum expansion for the most expanded region (middle panel); and percentage of cycling cells (lower panel). (F) The probability of detecting a given TCR clone in peripheral blood as a function of minimal clone size and mean pseudotime value of the clone. (G) Mean pseudotime values based on the categorization of clonotypes according to their principal region of enrichment. Data are presented as median ± first and third quartiles with whiskers depicting 95% confidence intervals. Statistical analysis by Tukey’s test. See also Figure S3.

Figure 4

Somatic mutation calling and the relationship with TCR clonotypic heterogeneity (A) Reconstructed phylogenies from WES of multi-regional LCM biopsies. Each node represents a mutant clone present in one or more of the biopsies. (B) Comparison of WES-derived phylogenies (left) with geographic location (center) and CD8⁺ TCR clonotype expansion (right). Colors reference somatic clones to spatial localization. Each column in the right panel represents a TCR clonotype; those with significant regional enrichment are highlighted in red. a, c, d, and e represent four different regions of the tumor core; g, tumor-normal interface; n, normal kidney; b, peripheral blood; h, normal adrenal gland. (C) Scatterplot of Mantel correlation between tree distances. x axis represents the correlation coefficient between WES-derived clones and TCR clonotype distances. y axis represents the correlation coefficient between spatial localization and TCR clonotype distances. (D) Dot plot showing inferred TCR groups with a high probability of sharing antigen specificity. The exhaustion levels represent pseudotime values. (E) Benchmarking results for scRNA-seq-derived calls against WES data for each patient donor. See also Figure S4.

Figure 5

Myeloid cell characterization, regional enrichment, and evolution (A) UMAP re-presentation of all myeloid cells, their annotation, and their regional contribution. Mono, monocyte; TR Mac, tissue-resident macrophage; TAM, tumor-associated macrophage. (B) The relative enrichment of different myeloid cell subsets across different regions sampled. (C) Dot plot depicting top DEGs for macrophage clusters. (D) Heatmap showing mean scaled scores for macrophage subsets by macrophage function of M1/M2 polarization and suppressive, angiogenesis, and phagocytosis activity. (E) Heatmap showing the results of pathway enrichments of macrophage subsets using gene set variation analysis. (F) UMAP with superimposed RNA velocity analysis of the monocyte and macrophage subsets, with zoomed-in windows highlighting possible directional flows from monocytes to macrophages. (G) Neighbor-joining tree depicting the relationship of different monocyte and macrophage clusters, utilizing the somatic mutations called from scRNA-seq data. The numbers of supporting votes in bootstrapping (100 times) are labeled. See also Figure S5.

Figure 6

RCC cell expression programs, regional enrichment, and prognosis (A) Heatmap showing expression programs derived in a representative patient using NMF. (B) Heatmap depicting shared expression meta-programs across all patients. (C) UMAP representing clusters of tumor cell population. (D) Relative expression scores of meta-programs in each RCC cell cluster (left) and the distributions of cells with different meta-programs in tumor core and tumor-normal interface. (E) Cells from patient donors PD45815 and PD45816, ranked by decreasing EMT score with corresponding PT score and cell location. (F) Spatial mapping of EMT and PT tumor cells in a representative tumor-normal interface sample (PD47171) using cell2location. Estimated abundance (color intensity) is overlaid on a histology image. Scale bars, 1 mm. (G) Box plots showing the EMT and PT scores of TCGA samples in different molecular subtypes, which are m1 (n = 145), m2 (n = 90), m3 (n = 93), and m4 (n = 86). Data are presented as median ± first and third quartiles with whiskers depicting 95% confidence intervals. ^∗∗∗p < 0.001 (two-sided Wilcoxon rank-sum test). (H) Dot plot showing gene scores of previously defined signatures in tumor cell populations. See also Figure S6 and Table S5.

Figure 7

Cellular interactions in the ccRCC microenvironment (A) Heatmap depicting the potential regulation of genes expressed by the EMT meta-program and ligands expressed by macrophages. (B) Spatial mapping of EMT tumor cells, PT tumor cells, and TR Mac.2 in Visium data for representative tumor-normal interface (PD45816) and tumor core (PD47171) samples using cell2location. Estimated abundance for cell types (color intensity) across locations (dots) is overlaid on histology images. Scale bars, 1 mm. (C) Abundance correlation between EMT tumor cells and TR Mac.2 in tumor cores and tumor-normal interfaces from multiple tumors. See also Figure S7.