Skip to main content

Decoding the multicellular ecosystem of vena caval tumor thrombus in clear cell renal cell carcinoma by single-cell RNA sequencing

Abstract

Background

Vascular invasion with tumor thrombus frequently occurs in advanced renal cell carcinoma (RCC). Thrombectomy is one of the most challenging surgeries with high rate of perioperative morbidity and mortality. However, the mechanisms driving tumor thrombus formation are poorly understood which is required for designing effective therapy for eliminating tumor thrombus.

Results

We perform single-cell RNA sequencing analysis of 19 surgical tissue specimens from 8 clear cell renal cell carcinoma (ccRCC) patients with tumor thrombus. We observe tumor thrombus has increased tissue resident CD8+ T cells with a progenitor exhausted phenotype compared with the matched primary tumors. Remarkably, macrophages, malignant cells, endothelial cells and myofibroblasts from TTs exhibit enhanced remodeling of the extracellular matrix. The macrophages and malignant cells from primary tumors represent proinflammatory states, but also increase the expression of immunosuppressive markers compared to tumor thrombus. Finally, differential gene expression and interaction analyses reveal that tumor-stroma interplay reshapes the extracellular matrix in tumor thrombus associated with poor survival.

Conclusions

Our comprehensive picture of the ecosystem of ccRCC with tumor thrombus provides deeper insights into the mechanisms of tumor thrombus formation, which may aid in the design of effective neoadjuvant therapy to promote downstaging of tumor thrombus and decrease the perioperative morbidity and mortality of thrombectomy.

Background

There were more than 431,280 new kidney cancer cases diagnosed in 2020 worldwide [1], of which the most common histological subtype is clear cell renal cell carcinoma (ccRCC), accounting for approximately 70–80% of all renal cell carcinoma (RCC) cases [2,3,4,5]. One unique clinical aspect of RCC is that it can invade through the renal vein into the inferior vena cava (IVC), and even grow up to the right cardiac chambers [6]. Venous tumor thrombus (TT) in the renal vein or inferior vena cava was reported in approximately 15% of RCC patients [7]. Reese et al. showed the tumor thrombus level does not necessarily affect disease-specific survival, but has a major impact on the complexity of surgery with a significant risk of hemorrhage, unstable hemodynamics and death [8,9,10,11]. Neoadjuvant systemic therapy may potentially decrease the TT burden and thus improve the safety and feasibility of tumor thrombectomy potentially improving the curative potential of surgical resection [12,13,14]. However, the underlying mechanism of TT formation is poorly understood, which prevents the design of effective neoadjuvant therapy for downstaging tumor thrombus and decreasing the perioperative morbidity and mortality of thrombectomy.

Bulk genomic studies, including ours, have revealed that most TTs with the propensity of rapid growth [15], harbored limited additional genomic alterations compared with matched primary tumors (PTs) [16,17,18]. The lack of fixation of new driver events in TTs may be due to their rapid extension and/or limited selective pressure in the intravascular space. Notably, we previously showed that the cases with TTs had distinct transcriptomic profiles compared to those without TTs by unsupervised clustering analysis of bulk RNA-seq data [16]. But a limited number of differentially expressed genes were identified between PTs and paired TTs, while multiple pathways related to the immune response and extracellular matrix and structure were significantly enriched in TTs [16]. Moreover, a higher proportion of macrophages in PTs of patients with TTs than in those of patients without TTs was observed by CIBERSORT analysis [16]. Thus, not only cancer cells but also tumor-infiltrating cells may participate in the process of tumor thrombus. And the cell type-specific gene expression patterns will promote the understanding of the intratumoral heterogeneity and the biology of tumor thrombus, and help discover more effective targets for tumor thrombus therapy.

In this study, we dissected the tumor ecosystem in 19 surgical tissue specimens from 8 ccRCC patients (discovery cohort) by single-cell RNA sequencing (scRNA-seq). Moreover, we included the bulk RNA-seq data of 60 tumor thrombus and paired primary tumor specimens from 30 ccRCC patients with TTs [16], the CheckMate 025 cohort [19, 20], the Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma (TCGA KIRC) and the IMmotion 150 cohort [21] as validation cohorts. We aimed to investigate the multicellular heterogeneity of tumor cell communities of tumor thrombus to provide valuable biological and clinical insights into this disease.

Results

scRNA-seq profiling of the tumor ecosystem in tumor thrombus and matched primary tumors

To elucidate the tumor ecosystem in venous tumor thrombus, we collected 19 surgical tissue specimens (discovery cohort) from treatment-naïve ccRCC patients (ccRCC, n = 8) for scRNA-seq. Tumor thrombus (n = 8), paired primary tumors (n = 8) and adjacent renal tissues (n = 3) were referred to as TTs, PTs and ARTs, respectively (Fig. 1a and Additional file 2: Table S1). Moreover, we included 4 different cohorts with bulk RNA-seq data for further validation (Fig. 1b). Approximately 1 billion unique transcripts were obtained from 140,805 cells: 62,035 cells from TTs, 58,695 cells from PTs, and 20,075 cells from ARTs. After filtering, a total of 120,780 cells were used for further analysis (Additional file 1: Fig. S1a-c). We identified and visualized 9 major cell types according to the expression of canonical gene markers using the Uniform Manifold Approximation and Projection (UMAP) method [22] (Fig. 1c-f and Additional file 1: Fig. S1a-c): epithelial cells, T cells, natural killer (NK) cells, myeloid cells, myofibroblasts, endothelial cells (ECs), B cells, plasma cells and mast cells. For the 16 PT and TT samples, all these cell subtypes were shared among patients and between PTs and TTs with no consistent change, albeit at different proportions (Fig. 1g). For the ART samples, in addition to the immune and stromal cells, we obtained a total of 6809 normal epithelial cells and further identified as the proximal tubule cells and the loop of Henle cells based on the classic markers, as described previously [22] (Fig. 1f-g and Additional file 1: Fig. S1d). In summary, our results revealed substantial heterogeneity of ecosystem compositions among the ART, PT and TT samples.

Fig. 1
figure 1

scRNA-seq profiling of the kidney primary tumors and tumor thrombus microenvironments. a, b Schematic diagram of the study strategy. The discovery cohort is shown in (a), and validation cohorts 1, 2, 3 and 4 are shown in (b). Validation cohort 1, 60 tumor thrombus and paired primary tumor specimens from 30 ccRCC patients. Validation cohort 2, 120 samples with advanced clear cell renal cell carcinoma with nivolumab treatment from CheckMate 025 cohort. Validation cohort 3, 531 samples of TCGA KIRC cohort. Validation cohorts 4, 82 samples with metastatic renal cell carcinoma with sunitinib treatment from IMmotion 150 cohort. c UMAP plot showing the annotation and color codes for cell types in the ccRCC ecosystem. Epithelial cell including normal epithelial cell (the proximal tubule cell and the loop of Henle cell) in ARTs and malignant epithelial cell in PTs and TTs. d, e The UMAP plot showing cell origins by color, patient origin (d), and ART, PT or TT origin (e). f Heatmap showing the expression of marker genes in the indicated cell types. g Histogram indicating the proportions of cells in tissues of each patient. N, adjacent renal tissues. C, primary tumors. T, tumor thrombus

T/NK cell clustering and state analysis identifies tumor thrombus are enriched with tissue resident CD8+ T cells in a progenitor exhausted state compared to primary tumors

Considering that recent studies have observed dramatic reduction of tumor thrombus in ccRCC patients after ICB therapy [23,24,25,26], we speculated that tumor-infiltrating immune cells may play important roles in the process. Therefore, we first performed unsupervised clustering of T and NK cells and obtained 12 clusters across ARTs, PTs and TTs, including four subtypes of CD4+ T cells (CD4-C1 to CD4-C4), five subclusters of CD8+ T cells (CD8-C1 to CD8-C5), two NK subclusters (NK1 and NK2) and one NKT cluster (Fig. 2a-d, Additional file 1: Fig. S2a-d and Additional file 3: Table S2).

Fig. 2
figure 2

scRNA-seq revealed heterogeneity in T/NK cells in primary tumors and tumor thrombus. a UMAP plot showing the subtypes of T/NK cells derived from ART, PT and TT samples. Each cluster is color-coded according to cell type. b UMAP plot illustrating T/NK cells clustered and color-coded according to each patient. c Bar plot illustrating the fraction of T/NK subgroups in ARTs, PTs and TTs. d Heatmap indicating the expression of selected gene sets in T/NK subtypes, including naïve, resident, inhibitory, cytokine, costimulatory, transcription factor (TF), and cell type. e Box plots illustrating the average proportion of CD4-C4, CD8-C1, and CD8-C3 subtypes among ARTs, PTs and TTs. p values were determined by a two-sided Wilcoxon test. f Violin and box plots showing the signature score distribution for progenitor exhausted and terminally exhausted CD8+ T cells within CD8-C1 to CD8-C4 subsets. p values were determined by a two-sided Wilcoxon test. ****, p < 0.0001. g Heatmap showing the fold change in the expression of the progenitor exhausted gene set of CD8-C1 to CD8-C4 cells in PTs compared with TTs. h Box plots representing the percentage of CD8-C1 and the scores of progenitor exhausted signature and CD8-C1 signature in paired PT and TT bulk RNA-seq samples. The CD8-C1 cell fraction per sample as inferred by CIBERSORTx. p values were determined by a two-sided Wilcoxon test

For CD4+ T cells, we identified naïve (CD4-C1; CCR7+), helper (CD4-C2; IL2+, CD4-C3; IL7+), and suppressive regulatory Treg (CD4-C4; FOXP3+) CD4+ T cells (Fig. 2d and Additional file 1: Fig. S2d). To determine the differences among ARTs, PTs and TTs, we calculated the percentage of each cluster in ARTs, PTs and TTs, and found that the relative percentage of naïve CD4+ T cells (CD4-C1) in PTs was reduced compared with those in ARTs and TTs (Additional file 1: Fig. S2e), and the proportion of suppressive Tregs (CD4-C4) was higher (Fig. 2e), indicating there might be in a more immune suppressive state in PTs.

For CD8+ T cells, we identified tissue resident (CD8-C1; CD69+), terminal exhausted (CD8-C2 and CD8-C3; PDCD1+ and TOX+), cycling (CD8-C4; MKI67+) and mucosal-associated invariant T (MAIT) (CD8-C5; IL7R+ and CCR4) CD8+ T cells (Fig. 2d and Additional file 1: Fig. S2c-d) according to the expression of canonical markers [27,28,29,30,31]. Similar with the proportion of CD4-C1 cells (Additional file 1: Fig. S2e), the abundance of CD8-C1 cells was also increased in TTs compared with PTs, but much less in both PTs and TTs than in ARTs (Fig. 2e). Although a considerable amount of tissue resident CD8+ T cells has been found in normal kidney tissues [32,33,34], the mechanisms for it are still largely unknown. Consistent with our finding, another study about hepatocellular carcinoma (HCC) also revealed the tissue resident CD8+ T cells were more abundant in adjacent non-tumor tissues than in tumors [35]. Furthermore, as observed previously in ccRCC [36], the proportion of exhausted CD8+ T cells (CD8-C3) was increased in tumor tissues, but with no significant difference between PTs and TTs (Fig. 2e and Additional file 1: Fig. S2e), indicating that CD8+ T cells are exhausted in both PTs and TTs. In pseudotime analysis, we removed MAIT cells, due to their differing TCR characteristics [35]. The continuous developmental trajectory of CD8-C1 to C4 represented a binary branched structure: CD8-C1 was the root, CD8-C3 and C4 were the end states of the two branches, and CD8-C2 was in a transitioning state (Additional file 1: Fig. S2f), which were very similar between PTs and TTs. Pathway analysis from the trajectory indicated T cell activation and lymphocyte differentiation was upregulated at an earlier stage, and cell cycle, ATP metabolism and hypoxia related pathways were highly enriched in the end state of CD8+ T cells (Additional file 1: Fig. S2g). These results suggested the CD8+ T cells shared the same transition trajectory between PTs and TTs, along with dysfunctional and metabolic disorders, and finally turned to be exhausted in both PTs and TTs.

There are evidences showed that certain CD8+ T cells in the progenitor exhausted state could enhance the efficacy of immune checkpoint blockade (ICB) therapy in melanoma and kidney cancer [37, 38]. Meanwhile, considering a striking regression of TTs with ICB therapy in some ccRCC patients [23,24,25,26], we asked whether there was certain subset of CD8+ T cells we recovered may resemble the progenitor exhausted phenotype. Then, we scored the CD8+ T cell subclusters for progenitor and terminally exhausted gene signatures according to the published studies [37, 38]. We noticed that the progenitor exhausted signature was significantly enriched in CD8-C1 subcluster, and the terminally exhausted signature was obviously lower than that of other subsets (Fig. 2f). Moreover, most of the progenitor exhausted signature genes showed upregulated expression in CD8-C1 cells in TT samples (Fig. 2g), which further demonstrated the progenitor exhausted phenotype was enriched in TTs. Furthermore, we examined a larger cohort of patients with paired PT and TT bulk transcriptomes and observed significantly increased proportions of the CD8-C1 subset in TTs compared with PTs, and the progenitor exhausted and CD8-C1 signature scores were both upregulated in TTs (Fig. 2h). Having identified CD8-C1 cells (the tissue resident CD8+ T cells) and its signature were enriched more in TTs than in PTs, we next examined whether CD8-C1 cells we found will show a better response to anti-PD-1 therapy by using the pretherapy bulk RNA-seq data from a larger cohort of ccRCC patients (CheckMate 025) treated with nivolumab (anti-PD-1) [19, 20]. As expected, the high levels of the CD8-C1 signature were associated with improved progression free survival (PFS) (Additional file 1: Fig. S2h), which was also supported by another study published recently that they observed a strong enrichment of the tissue resident CD8+ T cell cluster in one ccRCC patient with complete response to ICB [31]. Collectively, these results suggested that a high signature of tissue resident CD8+ T cells in both primary ccRCC and tumor thrombus may be the hint of a good response to ICB therapy.

Macrophages show enhanced ECM remodeling activity in tumor thrombus, and immune-reactive perturbations in primary tumors

In addition to adaptive immunity, the innate immune cells might be a first barrier for rapid extension of TTs in the intravascular space, which may have important impacts on tumor thrombus formation and resistance to ICB therapy. Thus, we analyzed the transcriptomes of 16,267 myeloid cells in PTs (n = 7689), TTs (n = 7129) and ARTs (n = 1449). Myeloid cells exhibited remarkable heterogeneity and were categorized into 14 clusters. Based on the expression of canonical markers, we annotated these clusters into six subtypes for macrophages (Macro1-Macro6), two for monocytes (Mono1-Mono2), four for DCs (DC1-DC4), one for neutrophils and one for cycling cells (Fig. 3a-c, Additional file 1: Fig. S3a-d and Additional file 4: Table S3). Macrophages were enriched in PTs and TTs compared with ARTs, except Macro4 cluster (Fig. 3b, d and Additional file 1: Fig. S3e). Neutrophils were remarkably depleted in both PTs and TTs compared with ARTs (Fig. 3d), implying tumor cells of PTs and TTs may escape elimination by circulating neutrophils. To better understand the roles of these myeloid-related populations, we first examined the immune checkpoint- and evasion-related gene expression levels. Almost all these genes were highly expressed in the myeloid subpopulations except neutrophils (Fig. 3e). While CD274 (PD-L1) and PDCD1LG2 (PD-L2), both ligands for PD-1 signaling mediating immune checkpoint in T cells, were detected sparsely across all myeloid subpopulations (Fig. 3e). These results suggested that neutrophils might be more effective in confronting cancer cells, and other macrophages, monocytes and DCs might already be ‘educated’ by PT or TT tumor cells.

Fig. 3
figure 3

Detailed characterization of myeloid cells in primary tumors and tumor thrombus. a UMAP plot showing the subtypes of myeloid-derived cells derived from ART, PT and TT samples. Each cluster is color-coded according to cell type. b Bar plot illustrating the average proportion of each myeloid subtype among ARTs, PTs and TTs. c Heatmap showing the expression of marker genes in each subtype of myeloid-derived cells. d Box plots illustrating the fraction of the subclusters of Macro1–3 and neutrophils in ARTs, PTs and TTs. p values were determined by a two-sided Wilcoxon test. e Dot plot showing the percentage of cells and expression level in each cell type expressing immune checkpoint and evasion-related genes. f Violin plots representing the expression levels of M1, M2 and TAM marker genes in macrophage subtypes. p values were determined by a two-sided Wilcoxon test. g Heatmap of GSEA scores indicating the pathways significantly enriched in PTs or TTs. INF-γ, interferon- γ. INF-α, interferon- α. Antigen Proc. and Pres., antigen presentation and processing. Extra. structure organization, extracellular structure organization. EMT, epithelial-mesenchymal transition. h Violin and box plots comparing the expression distributions of immune checkpoint and evasion-related genes and ECM remodeling and angiogenesis-related genes between PTs and TTs. p values were determined by a two-sided Wilcoxon test. ***, p < 0.001

Macrophages are usually classified into two canonical subtypes, proinflammatory M1 and anti-inflammatory M2 [39,40,41]. However, we could not clearly distinguish M1 and M2 macrophages by known marker genes such as CD86, TLR4 (M1) and CD163 and MSR1 (M2), as they were all expressed in Macro1-Macro3 and Macro5-Macro6 but poorly expressed in Macro4 (Fig. 3f). While the tumor-associated macrophage (TAM) markers, as well as CD68 and HLA-DRA as both M1/M2 and tumor-associated macrophage (TAM) markers were highly expressed in each macrophage subtype (Fig. 3f), indicating they were all TAMs. Intriguingly, we noticed a significant enrichment of gene expression signatures in the proinflammatory phenotype, such as response to interferon, and antigen presentation pathways in TAMs of PTs compared to that of TTs (Fig. 3g and Additional file 1: Fig. S3f). Surprisingly, these TAMs in PTs also upregulated expression of immune checkpoint and evasion markers compared to TTs (Fig. 3h). Thus, TAMs in PTs exhibited both proinflammatory and immune suppressive phenotypes, potentially contributing to tumor cells immune escape. Whereas, in TTs, we observed increased signature scores associated with tumor progression-related pathways, including epithelial-mesenchymal transition (EMT), extracellular structure organization and angiogenesis (Fig. 3g-h and Additional file 1: Fig. S3f). These findings were consistent with previous reports that macrophages act as shapers of the tumor ECM to facilitate cell migration and angiogenesis [42, 43]. Collectively, the complexity of TAMs demonstrated distinct dysfunctional states between PTs and TTs, with a more immunosuppressive phenotype in PTs, and enhanced ECM remodeling in TTs.

Malignant cells exhibited increased ECM remodeling in tumor thrombus, and a more immunosuppressive phenotype in primary tumors

To better understand cellular programs active in cancer cells that may function together with the immune cells, we next sought to identify the expression patterns of the malignant cells between PTs and TTs. The tumor cells were confirmed by the detection of copy number variations (CNVs), inferred by scRNA-seq data, using myofibroblasts and endothelial cells as references, according to a published study [44]. Seven of 8 patients showed similar CNV patterns and subclonal structures between TTs and PTs, except patient P06 (Fig. 4a-b and Additional file 1: Fig. S4a), supporting the previous studies of limited additional genomic alterations in TTs [16, 18]. Then, we evaluated the transcriptional difference of malignant cells between PTs and TTs. In total, we obtained 15,012 malignant epithelial cells from PTs and 12,078 cells from TT specimens. Since the proximal tubule cells has been identified as the origin of ccRCC [22], we then compared the hypoxia and angiogenesis signatures, which were two of classic biological states contributed to carcinogenesis, across the malignant cells and the proximal tubule cells. Consistent with previous studies, both PT and TT malignant cells had significantly higher hypoxia and angiogenesis signatures than cells from ARTs (Additional file 1: Fig. S4b). Then, focusing on the malignant cells, we totally obtained 4 clusters, and the proportion in individual patients was different in each subcluster, further supporting a high degree of intertumoral heterogeneity (Fig. 4c-e and Additional file 1: Fig. S4c-f). To better understand the transcriptional heterogeneity of them, we scored the four clusters by GSVA Hallmark analysis, and observed dramatic transcriptional program differences among each cluster (Fig. 4f). In detail, cluster 1 was the most abundant and characterized with increased hypoxia, glycolysis, oxidative phosphorylation (OXPHOS), along with concomitant immune activation and metabolic pathways (Fig. 4f). The coexistence of glycolysis and OXPHOS metabolic states in ccRCC has been reported in another study [38] in which the hybrid phenotype in RCC contributes to metabolic plasticity, allowing cancer cells to adapt to various microenvironment [45,46,47]. Furthermore, interestingly, similar with the TAMs, we observed an enrichment of genes involved in interferon-γ (IFN-γ) response and antigen presenting pathways with upregulated expression of immune checkpoint and evasion genes in PTs compared to TTs (Fig. 4g-h, Additional file 1: Fig. S4g-h and Additional file 5: Table S4). It demonstrated their capacity to promote immune escape and repress the anti-tumor immune response. While, the TT tumor cells were activated in ECM remodeling and response to metal ion related pathways (Fig. 4g and Additional file 1: Fig. S4g-h). Consistent with these results, a significant higher proportion of CD274+ malignant cells in PTs, and more COL4A1+ tumor cells in TTs were observed (Fig. 4i), indicating the higher immunosuppressive state in PTs and the dynamic remodeling of ECM in TTs. Furthermore, we found the classic ferroptosis suppressor gene GPX4 showed upregulated expression in TTs (Additional file 1: Fig. S4h), implying inducing ferroptosis through regulation of GPX4 for tumor therapy might be more effective in TTs than in PTs in ccRCC patients. In summary, at single-cell resolution, we revealed the malignant cells in PTs demonstrated a stronger immunosuppressive phenotype, while tumor cells in TTs showed a higher ECM reprogramming state.

Fig. 4
figure 4

Identification and characterization of malignant cells in primary tumors and tumor thrombus. a, b Heatmaps showing large-scale CNVs for individual malignant cells (rows) from P02 (left) (a) and P06 (right) (b). Endothelial cells and myofibroblasts were treated as references (top), and malignant cells were observed in PTs and TTs (bottom). The color bar red indicates amplifications, and blue indicates deletions. c UMAP plot representing the subtypes of malignant cells from PT and TT samples. d UMAP plot showing malignant cells derived from PTs and TTs, colored by tissue origin. e UMAP plot illustrating color-coded tumor cells, according to each patient. f Differentially expressed pathways were scored per cell by GSVA among 4 malignant cell subtypes. g GSEA of the interferon-γ (IFN-γ) signaling pathway and antigen presentation and processing (Antigen Proc. and Pres.) enrichment scores in all malignant cells in PTs compared with in TTs; collagen containing extracellular matrix (Collagen conta. Extra. matrix) and response to metal ion enrichment scores in all tumor cells in TTs compared with in PTs. h Heatmap depicting the differential expression fold change of PTs compared with TTs for immune checkpoint and evasion-related genes. i Box plots showing the fraction of CD274+ and COL4A1+ tumor cells in PT and TT single-cell samples. p values were determined by a two-sided Wilcoxon test

Endothelial cells represent upregulation of ECM remodeling related pathways in tumor thrombus, and highly-abundant CCL4 + and NDUFA4L2 + endothelial cells are associated with poor survival

To further evaluate endothelial cells roles that may facilitate tumor thrombus growth and extension in the intravascular space, we focused on endothelial cells (ECs), and obtained 16,221 ECs from PTs (n = 4808), TTs (n = 7915) and ARTs (n = 3498). These cells were subclustered into six different clusters (Fig. 5a-c and Additional file 1: Fig. S5a-c), and each of them was identified based on the classically expressed genes in ECs [48, 49], including glomerular like ECs (Endo1: SOST+), cancer-related ECs (Endo2: NDUFA4L2+), arterial ECs (Endo3: GJA5+), ACKR1+ ECs (Endo4: ACKR1+), the tip cells (Endo6: CXCR4+), and one undefined cluster (Endo5) (Fig. 5a-c, Additional file 1: Fig. S5a-c and Additional file 6: Table S5). Furthermore, we found the ratio of Endo2 in TTs was higher than that in ARTs and PTs (Fig. 5d and Additional file 1: Fig. S5d), suggesting that cancer related vasculature formation may facilitate the rapid extension of TTs in the intravascular space. In support, we observed significantly increased percentages of Endo2 subset in TTs through analyzing the bulk transcriptomes from a large cohort of patients with 30 paired PTs and TTs (Fig. 5e). Moreover, the high levels of Endo2 gene signature were significantly associated with poor survival in TCGA KIRC cohort (Fig. 5f and Additional file 1: Fig. S5e). Interestingly, like TAMs and malignant cells, we found significantly increased expression signatures associated with interferon response and antigen binding pathways in ECs of PTs (Fig. 5g and Additional file 1: Fig. S5f-g), but along with variation in genes expression related to immunosuppression among the subsets of ECs between PTs and TTs (Additional file 1: Fig. S5h). At the same time, we noticed in Endo2 subset, the expression level of most of genes related to immunosuppression was upregulated in PTs compared with TTs (Additional file 1: Fig. S5h), suggesting these cells in PTs might be in a more immunosuppressive state. Furthermore, we found that the ECs had enhanced ECM remodeling and cell proliferation activities in TTs (Fig. 5g and Additional file 1: Fig. S5f-g). Thus, these data indicated that in addition to TAMs and cancer cells in TTs, endothelial cells in TTs may modify the ECM to facilitate tumor thrombus growth, and the molecular therapies targeting the cells with Endo2 signature might be able to specifically downstage the vena cava thrombus.

Fig. 5
figure 5

Detailed characterizations of endothelial cells and myofibroblasts in ARTs, PTs and TTs. a UMAP plot showing the subtypes of endothelial cells derived from ART, PT and TT samples. Each cluster is color-coded according to cell type. b Bar plot illustrating the average proportion of each endothelial subtype among ARTs, PTs and TTs. c Heatmap showing the expression of marker genes in each subtype of endothelial cells. d Box plots illustrating the fraction of Endo2 subgroup in ARTs, PTs and TTs. p values were determined by a two-sided Wilcoxon test. e Box plots showing the percentage of Endo2 in paired PT and TT bulk RNA-seq samples, inferred by CIBERSORTx. p values were determined by a two-sided Wilcoxon test. f Kaplan-Meier plot showing that KIRC patients in the TCGA dataset with high expression of Endo2 signature had shorter progression free survival (PFS). g Heatmap of GSEA scores indicating the pathways significantly enriched in PTs or TTs. Extra. structure organization, extracellular structure organization. h UMAP plot showing the subtypes of myofibroblasts derived from the ART, PT and TT samples. Each cluster is color-coded according to cell type. i Bar plot illustrating the average proportion of each myofibroblast subtype among ARTs, PTs and TTs. j Heatmap showing the expression of marker genes in each subtype of myofibroblasts. k Box plot illustrating the average proportion of Myo2 subcluster among ARTs, PTs and TTs. p values were determined by a two-sided Wilcoxon test. l Box plot representing the proportions of subclusters of Myo2 in myofibroblasts between normal (N) and tumor (T) samples in the TCGA cohort. p values were determined by a two-sided Wilcoxon test. m Kaplan-Meier plot showing that KIRC patients in the TCGA dataset with high expression of CD248 had worse progression free survival (PFS). n Kaplan-Meier plot showing that patients from the sunitinib arm of the IMmotion150 cohort with a high Myo2 signature score had improved progression free survival. o Heatmap of GSEA enrichment scores for indicating gene sets significantly enriched in PTs or TTs. Collagen fibril organ, collagen fibril organization. Extra. structure organization, extracellular structure organization. EMT, epithelial- mesenchymal transition. p Violin and box plots demonstrating the expression levels of CD248, COL1A1, COL4A1, and FN1 between PTs and TTs. p values were determined by a two-sided Wilcoxon test. ***, p < 0.001

Myofibroblasts in tumor thrombus enhance ECM remodeling with increased extracellular matrix production compared with primary tumors

Fibroblasts have been shown to be a major contributor to ECM remodeling in tumor progression in previously studies. Thus, we further investigated fibroblasts in ARTs, PTs and TTs. Almost all fibroblasts were positive for α-SMA, a conventional marker of myofibroblasts (Fig. 1f, Additional file 1: Fig. S1c), which is consistent with the findings in a previous ccRCC report [36]. We found five distinct subtypes by reclustering 10,726 myofibroblasts, according to the published studies [36, 50] (Fig. 5h-j, Additional file 1: Fig. S5i-k and Additional file 6: Table S5). In detail, our findings identified a state of stress with highly expressed heat shock proteins in Myo1, a highly vascular pericyte-like phenotype in Myo2, ACTG2+ myofibroblasts association in Myo3, antigen presentation phenotype in Myo4, and the characteristics of cancer-associated myofibroblasts in Myo5 [50]. Moreover, we found that vascular pericyte-like Myo2 was enriched in TTs and PTs compared to ARTs (Fig. 5k and Additional file 1: Fig. S5l), suggesting that Myo2 may help to remodel tumor vessels toward a mature phenotype in tumor thrombus, since pericytes are necessary for vessel maturation [51, 52]. Furthermore, we validated this result in the TCGA KIRC cohort. Indeed, the percentage of Myo2 increased significantly in the tumor tissues of KIRC patients compared to normal tissues (Fig. 5l). In addition, the marker gene of Myo2, CD248, was remarkably upregulated in tumors compared to the adjacent tissues (Additional file 1: Fig. S5m), and significantly associated with poor survival in KIRC patients (Fig. 5m and Additional file 1: Fig. S5n). Antiangiogenic tyrosine kinase inhibitors (TKIs) have been broadly applied in ccRCC patients, and pericytes play crucial roles in integrity of the tumor microvasculature [52,53,54,55]. Therefore, we further evaluated whether the proportion of Myo2 could predict the response to TKIs. We found in the IMmotion150 cohort [21], patients with sunitinib treatment had an improved PFS with high score of Myo2 signature (Fig. 5n), indicating that Myo2 might be a good predictor of response to TKIs in ccRCC patients.

Next, we examined whether there were also transcriptional differences in myofibroblasts between TTs and PTs, and found all the myofibroblast subclusters in PTs were enriched in the TNFα and UV response pathways, while almost all myofibroblasts in TTs increased the signatures associated with angiogenesis and ECM remodeling, including extracellular matrix organization, collagen fibril organization and EMT pathways (Fig. 5o-p and Additional file 1: Fig. S5o). Therefore, myofibroblasts, like TAMs, malignant cells and ECs, may also play promoting roles in ECM remodeling of TTs, thereby possibly representing a potential target for ccRCC patients with TTs.

Tumor-stroma interplay reshaped the extracellular matrix in tumor thrombus associated with poor survival

Given the observations of activated ECM remodeling as common trends in macrophages, malignant cells, ECs and myofibroblasts of TTs compared to PTs, we hypothesized that the different cell populations participated in complex crosstalk. To identify possible non-cell-autonomous effects, we used CellPhoneDB to identify putative signaling between different cell populations via known receptor-ligand pairs (Additional file 7: Table S6). Notably, the interactions in multiple ligand-receptor pairs within the ECM remodeling pathway were stronger in TTs than in ARTs and PTs, and myofibroblasts were predicted to have the strongest signals of significant interactions with tumor cells and ECs, regardless of subgroups of them (Fig. 6a and Additional file 1: Fig. S6a). We additionally validated this finding in the validation cohort with 30 paired PT and TT bulk RNA-seq data. Consistent with the upregulated expression of ECM components and their receptors in stromal cells and other types of cells in TTs, the ECM assembly signature was robustly more strongly correlated with the CIBERSORTx-estimated myofibroblast fraction in TTs compared with in PTs (Fig. 6b and Additional file 1: Fig. S6b), supporting that the ECM remodeling related pathways were more activated in TTs than that in PTs.

Fig. 6
figure 6

Tumor-stroma interplay reshaped the extracellular matrix in tumor thrombus associated with poor survival. a Dot plot showing inferred interactions between epithelial cells (proximal tubule cells from ARTs and malignant cells from PTs and TTs) and macrophages, endothelial cells, and myofibroblasts. Circle size indicates the significance of the interaction, and circle color indicates the mean expression of ligand and receptor genes. The red letters represent ligands, and the black letters represent receptors. b Scatterplots demonstrating the extracellular matrix assembly (left) and EMT (right) signature score versus total myofibroblast cell fraction in bulk RNA-seq of PT-TT paired samples in the validation cohort. Cell-type fractions were inferred using CIBERSORTx. The Pearson coefficient (R) and associated p value are reported for each correlation. c Scatterplot depicting inferred cell-to-cell interactions between tumor cells and myofibroblasts by combining scRNA-seq data and bulk RNA-seq data. d Box plot showing the calculated correlation between extracellular matrix remodeling-related genes and myofibroblast relative abundance involved in PTs and TTs. p values were determined by a two-sided Wilcoxon test. e Scatterplot representing inferred cell-to-cell interactions between myofibroblasts and endothelial cells by combining scRNA-seq data and bulk RNA-seq data. f Box plot showing the calculated correlation between extracellular matrix remodeling-related genes and endothelial cell relative abundance involved in PTs and TTs. p values were determined by a two-sided Wilcoxon test. g Signature scores for extracellular matrix assembly (left) and collagen-containing extracellular matrix (right) in paired PT-TT bulk RNA-seq samples in the validation cohort. p values were determined by a two-sided Wilcoxon test. h Progression free survival for the TCGA KIRC cohort based on high tumor-stroma and stroma-tumor interaction signatures versus low signature expression

Given the strong interactions identified between myofibroblasts and malignant cells and endothelial cells, we first examined whether tumor cells governed stromal myofibroblasts recruitment. We combined the single-cell transcriptomes with bulk RNA-seq according to previous studies [56, 57]. For example, we searched for genes increased in tumor cells rather than in myofibroblasts in single-cell data that were correlated with myofibroblast abundance. In detail, we found extensive interactions between myofibroblasts and tumor cells. In detail, tumor cells upregulated COL20A1, COL28A1 and TGFB1 expression to recruit more myofibroblasts in TTs compared with in PTs (Fig. 6c). Furthermore, myofibroblasts released collagen-related genes such as COL6A2, COL1A2 and COL4A2 to promote extracellular matrix assembly and have higher correlation with the abundance of myofibroblasts in TTs (Fig. 6c). The correlation between the relative abundance of myofibroblasts and extracellular matrix remodeling-related genes up-regulated in tumor cells was significantly higher in TTs than in PTs (Fig. 6d and Additional file 1: Fig. S6c). These data suggest that tumor cells might recruit more myofibroblasts in TTs by secreting ECM-related components. Additionally, we also found that myofibroblasts conferred enhanced recruitment of ECs in TTs than that in PTs, which might potentially facilitate angiogenesis (Fig. 6e-f and Additional file 1: Fig. S6d). Furthermore, the expression signature scores of ECM assembly, collagen-containing ECM and extracellular structure organization were all at a higher level in TTs than that in PTs in the validation cohort of 30 paired PT and TT ccRCC patients (Fig. 6g and Additional file 1: Fig. S6e). Thus, we showed that the stronger tumor-stroma interplay in the tumor thrombus was accompanied by the upregulated expression of a specific set of ECM related genes and pathways in both tumor cells and stromal cells.

To further investigate the association between tumor-stroma interplay and patient prognosis, we defined two interaction signature scores using significant L-R pairs to calculate the strengths of the ECM remodeling-related interactions in the TCGA KIRC cohort. The tumor-stroma interaction score was calculated between tumor cells and myofibroblasts, endothelial cells, and macrophages, while, the stroma-tumor interaction score was generated between myofibroblasts and tumor cells, endothelial cells, and macrophages. Notably, we found a stronger tumor-stroma interplay was associated with poor survival (Fig. 6h and Additional file 1: Fig. S6f). In summary, tumor-stroma interplay reshaped the extracellular matrix in tumor thrombus associated with poor survival.

Discussion

Surgical treatment of locally advanced ccRCC with TTs, especially with those that extend to the inferior vena cava, remains a clinical challenge. Neoadjuvant targeted therapies, such as VEGFR or mTOR inhibitors, have been shown to reduce tumor thrombus levels and potentially make thrombectomy more feasible. However, the percentage of benefited patients showing thrombus regression from these treatments were extremely variable (from 25 to 67%) in different clinical studies [12]. Specifically, neoadjuvant immunotherapy was only effectively reducing the tumor thrombus size in some patients [23,24,25,26]. Thus, systematically deciphering the molecular differences between tumor thrombus and primary tumors is urgently needed for the development of new predictive markers and more effective therapies for ccRCC with TTs. Here, we presented a comprehensive single-cell transcriptomic atlas characterizing the heterogeneity of tumor cells, immune cells, and stromal cells in primary tumors and paired tumor thrombus. Cell type-specific biological features associated with the rapid extension of tumor thrombus in the intravascular space were also identified, which are potentially useful in designing more effective therapies for patients with tumor thrombus.

Our results showed that macrophages, malignant cells, endothelial cells and myofibroblasts in TTs exhibited enhanced remodeling of the extracellular matrix pathways compared to matched primary cancer cells, and the tumor-stroma interplay reshaped the extracellular matrix in tumor thrombus and was associated with poor survival. These results suggest that malignant cells of TTs are not independently suspended in the blood vessels, but require to form an ECM network embedded with immune cells and stroma cells to promote tumor thrombus growth. Although the cancer cells of TTs have intravasated into vein, and might have more opportunity to spread to lung, there is no direct experimental or clinical evidence that ccRCC patients with tumor thrombi are more likely to develop pulmonary metastasis. Previous studies demonstrated the ECM is important at initial stages of metastasis, where interactions between tumor cells and extracellular matrix induce an invasive phenotype [58, 59]. Moreover, TT grade has been reported as an independent predictor for metastatic outcome irrespective of the grade of the PTs [60], which suggests the invasion ability of tumor cells in TTs is very important for metastatic potential. Our study demonstrated that ECM remodeling is crucial for both metastasis and tumor thrombus growth. Thus, targeting ECM remodeling factors associating with TT formation might represent effective counter-measurements in ccRCC patients with TTs.

The high efficacy of ICB therapies in some ccRCC patients with tumor thrombus suggest the need to identify what kind of features or markers are likely to benefit from ICBs, based on a better understanding of the immune microenvironment between the primary tumor and tumor thrombus. Our data revealed that TAMs did not polarize to either M1-like or M2-like subsets, but showed distinct molecular and metabolic phenotypes with no subset diversity between PTs and TTs. Interestingly, TAMs from PTs showed both pro-inflammatory and anti-inflammatory phenotypes compared with TTs, suggesting TAMs in PTs may potentially facilitate the escape of tumor cells from immune surveillance. Importantly, more abundant subcluster of CD8+ T cells with the progenitor exhausted state in TTs was further validation in a bulk RNA-seq cohort. These cells have similarities to those described as mediating ICB response in melanoma and advanced RCC patients [38, 61]. Our findings may provide a mechanistic explanation as to why some patients with TTs respond to anti-PD-1 immunotherapy effectively. Moreover, we detected an association between pre-therapy levels of tissue resident CD8+ T cell signature and better response to anti-PD-1 therapy, suggesting that anti-PD-1 immunotherapy may obtain a greater benefit in reducing both primary tumor and tumor thrombus with higher proportions of tissue resident CD8+ T cells. Furthermore, these findings also suggest a potential value of tissue resident CD8+ T cell signature as a significant indicator for ICB response in ccRCC patients. Thus, strategies enhancing tissue resident CD8+ T cells may have good application prospects in immunotherapy for ccRCC patients with TTs. Moreover, it is necessary for further clinical and mechanistic investigation regarding the associations between tissue resident CD8+ T cells with ICB response in ccRCC patients with TTs given the small sample size of our single-cell and validation cohorts.

Conclusions

In summary, this study provides evidence of phenotypic heterogeneity between primary tumors and tumor thrombus, in terms of tumor cells, immune cells, and stromal cells. Our data can be a valuable resource, facilitating a deeper understanding of the mechanisms associated with tumor thrombus and assisting in developing more effective neoadjuvant molecular therapies and biomarkers for advanced ccRCC patients with TT.

Methods

Clinical information and sample collection

All samples were obtained from Peking University Third Hospital, Beijing, China. Patient characteristics and clinical information are shown in Additional file 2 Table S1.

Single-cell suspension preparation and sequencing

Single-cell suspensions for single-cell RNA-seq were obtained by mechanical and enzymatic dissociation. According to the manufacturer’s protocol, the Single Cell 3′ Library and Gel Bead Kit V3.1 (10X Genomics, 1,000,075) and Chromium Single Cell B Chip Kit (10x Genomics, 1,000,074) were used to prepare barcoded scRNA-seq. The libraries were finally sequenced using an Illumina NovaSeq6000 sequencer with a sequencing depth of at least 100,000 reads per cell with a paired-end 150 bp (PE150) reading strategy (performed by Capital Bio Technology, Beijing).

scRNA-seq data processing

CellRanger (version 4.0.0) coupled with human reference version GRCh38 was used to process the 10X single-cell RNA-seq raw data for each sample, following the DoubletFinder R package (version 1.2.2) [62] to computationally infer and remove doublets in each sample individually, with default parameters. After removal of doublets, we employed the Seurat R package (version 3.2.2) [63] to analyze the output-filtered gene expression matrices. In brief, low-quality cells were removed if they met the following conditions: (i) > 5000 or < 200 genes and (ii) > 50% UMIs derived from the mitochondrial genome according published papers [22, 48, 64,65,66,67]. Then, according to the published studies [50, 65, 68,69,70], the filtered expression matrix was then normalized with the function “NormalizeData”, followed by the identification of 2000 genes of high cell-to-cell variation by using the function “FindVariableFeatures”. For multi-sample integration, we employed the function “FindIntegrationAnchors” to obtain “anchors” across individual samples. By inputting the anchors into the function “IntegrateData” and regressing out the influence of library size, percentage of mitochondria genes and cell-cycle scores, we created a “batch-corrected” expression matrix of all cells on the 2000 highly-variable genes. Based on the batch-corrected data, we performed Principal Component Analysis (PCA) with top 2000 variable features by using the function “RunPCA”. Cells were then clustered using the functions “FindNeighbors” and “FindClusters” with the first 50 principal components (PCs). Finally, we conducted nonlinear dimensional reduction for data visualization. In brief, UMAP was performed on the top 50 PCs by using the function “RunUMAP”.

Cell type annotation and cluster marker identification

The function “FindAllMarkers” function was used to find markers for each of the identified clusters and annotated on the basis of the expression of canonical markers of particular cell types and the annotation reference created by SingleR (version 1.2.4) [71] based on the published single-cell transcriptome data of kidney cancer [22].

Differential gene expression analysis

The parameter MAST in the function “FindAllMarkers” was used to perform differential gene expression analysis. A gene was considered significantly different with adjusted p < 0.05 [72].

CNV estimation, identification of malignant cells and malignant cell subset analysis

To infer CNV patterns from the scRNA-seq data, we used an approach described on the website tutorial (https://github.com/broadinstitute/inferCNV). The identification of malignant cells was performed based on previous reports [44, 73]. Endothelial cells and myofibroblasts were considered references. And 30% of them (as spike-in cells) were randomly selected together with epithelial cells for CNV inference and hierarchical clustering in each patient. We considered the epithelial cells that clustered together with spike-in control cells to be ‘CNV-low’ cells, whereas the remaining cells were considered ‘CNV-high’, as malignant cells for further analysis.

Diversity score

We calculated the diversity score of a tumor based on the gene expression profiles of malignant cells according to a published literature from Ma et al. [74]. The PCA analysis was performed to project the original expression profiles of all malignant cells to reducing the dimensionality of such datasets and to derive PCs, which could increase interpretability but at the same time minimize information loss. Then we calculated tumor diversity score to measure the degree of intratumoral heterogeneity based on the PCs within tumor by referring to the diversity score algorithms.

Gene set enrichment analysis

We applied gene set enrichment analysis (GSEA) to analyze different pathways in different subclusters on the hallmark pathways and GO terms documented in the molecular signature database [75,76,77]. We applied the GSVA R package (version 1.38.0) to estimate pathway activity scores for single cell [78]. The differential activities of pathways were calculated using the limma R package (version 3.46.0) [79].

Defining cell state scores

We used cell scores to evaluate the degree to which individual cell expressed a certain predefined expression gene set. To define progenitor exhausted and terminally exhausted phenotypes, we used the function “AddModuleScore” to calculate cell scores by using well-defined progenitor exhausted and terminally exhausted CD8+ T cell signatures [38, 61].

Pseudotime trajectory analysis

To characterize the potential process of CD8+ T functional changes and determine the potential lineage differentiation, we applied the Monocle2 (version 2.16.0) R package [80] excluding MAIT cells. The gene-cell matrix of UMI counts was provided as the input to Monocle, and then, the newCellDataSet function was employed to create a CellDataSet with the parameter expressionFamily = negbinomial.size. We then performed the differentialGeneTest function to identify significantly different genes over time.

Investigating genes correlated with certain cellular abundances

We estimated the abundance of each cell type according to the average expression levels of the top 50 cell-specific marker genes with average expression > 2 in each cell type in the bulk expression profiles. In addition, we manually excluded genes were obviously expressed in many other cell types [57]. Then, we calculated correlations between each gene and the abundance of certain cell types by using the function “corr” in the validation data. Furthermore, we identified different expression genes (DEGs) between two indicated cell types from single-cell profiles with log 2-transformed expression ratios > 1.5 or < − 1.5 and compared the relative correlation of the abundance of myofibroblasts or ECs associated with extracellular matrix remodeling-related genes between PTs and TTs.

Quantification and statistical analysis

We performed all the statistical analyses using R software (version 4.0.0). For comparison of the signature scores or CIBERSORTx-inferred immune and nonimmune fractions between different cell groups and bulk RNA-Seq sample groups, a two-sided Wilcoxon test was used. Detailed descriptions of the statistical tests performed for individual analysis are provided in the Figure legends and Methods.

Availability of data and materials

The raw single cell RNA sequencing data has been deposited in the Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences. The accession number for the sequencing data reported in this paper is HRA000963 at https://ngdc.cncb.ac.cn/gsa-human [81]. For the validation cohort with PT and paired TT, raw bulk RNA-seq data were obtained from the database of NCBI Sequence Read Archive (SRA) under the accession code PRJNA596338 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA596338) [82]. For the CheckMate 025 cohort, normalized bulk RNA-seq and clinical data were obtained from published results [19, 20]. For TCGA-KIRC cohort, preprocessed bulk RNA-seq data were downloaded from UCSC Xena (https://xena.ucsc.edu/). For the IMmotion150 cohort, the datasets are available from the European Genome-Phenome Archive (EGA) repository, the accession number is EGAS00001002928 (https://ega-archive.org/studies/EGAS00001002928) [83]. All code used to analyze data are available on GitHub, at https://github.com/zhangqi234/RCC-tumor-thrombus-scRNA.

References

  1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49.

    PubMed  PubMed Central  Google Scholar 

  2. Lopez-Beltran A, Carrasco JC, Cheng L, Scarpelli M, Kirkali Z, Montironi R. 2009 update on the classification of renal epithelial tumors in adults. Int J Urol. 2009;16:432–43.

    PubMed  Google Scholar 

  3. Leibovich BC, Lohse CM, Crispen PL, Boorjian SA, Thompson RH, Blute ML, et al. Histological subtype is an independent predictor of outcome for patients with renal cell carcinoma. J Urol. 2010;183:1309–15.

    PubMed  Google Scholar 

  4. Teloken PE, Thompson RH, Tickoo SK, Cronin A, Savage C, Reuter VE, et al. Prognostic impact of histological subtype on surgically treated localized renal cell carcinoma. J Urol. 2009;182:2132–6.

    PubMed  PubMed Central  Google Scholar 

  5. Kim JK, Kim TK, Ahn HJ, Kim CS, Kim KR, Cho KS. Differentiation of subtypes of renal cell carcinoma on helical CT scans. AJR Am J Roentgenol. 2002;178:1499–506.

    PubMed  Google Scholar 

  6. Noguchi K, Hori D, Nomura Y, Tanaka H. Renal cell carcinoma with tumor-thrombus extension into the right ventricle. Ann Vasc Dis. 2012;5:376–80.

    PubMed  PubMed Central  Google Scholar 

  7. Psutka SP, Leibovich BC. Management of inferior vena cava tumor thrombus in locally advanced renal cell carcinoma. Ther Adv Urol. 2015;7:216–29.

    PubMed  PubMed Central  Google Scholar 

  8. Reese AC, Whitson JM, Meng MV. Natural history of untreated renal cell carcinoma with venous tumor thrombus. Urol Oncol. 2013;31:1305–9.

    PubMed  Google Scholar 

  9. Skinner DG, Pfister RF, Colvin R. Extension of renal cell carcinoma into the vena cava: the rationale for aggressive surgical management. J Urol. 1972;107:711–6.

    CAS  PubMed  Google Scholar 

  10. Skinner DG, Pritchett TR, Lieskovsky G, Boyd SD, Stiles QR. Vena caval involvement by renal cell carcinoma. Surgical resection provides meaningful long-term survival. Ann Surg. 1989;210:387–392; discussion 392-384.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Berczi A, Flasko T, Szerafin T, Thomas B, Bacso Z, Berczi C. Surgical management and outcome of renal cell carcinoma with inferior vena cava tumor thrombus. Urol Int. 2017;99:267–71.

    PubMed  Google Scholar 

  12. Dey S, Peabody HN, Noyes SL, Lane BR. Neoadjuvant targeted molecular therapy before renal surgery. Urol Clin North Am. 2017;44:289–303.

    PubMed  Google Scholar 

  13. Borregales LD, Adibi M, Thomas AZ, Wood CG, Karam JA. The role of neoadjuvant therapy in the management of locally advanced renal cell carcinoma. Ther Adv Urol. 2016;8:130–41.

    CAS  PubMed  Google Scholar 

  14. Cost NG, Delacroix SE Jr, Sleeper JP, Smith PJ, Youssef RF, Chapin BF, et al. The impact of targeted molecular therapies on the level of renal cell carcinoma vena caval tumor thrombus. Eur Urol. 2011;59:912–8.

    CAS  PubMed  Google Scholar 

  15. Woodruff DY, Van Veldhuizen P, Muehlebach G, Johnson P, Williamson T, Holzbeierlein JM. The perioperative management of an inferior vena caval tumor thrombus in patients with renal cell carcinoma. Urol Oncol. 2013;31:517–21.

    PubMed  Google Scholar 

  16. Wang XM, Lu Y, Song YM, Dong J, Li RY, Wang GL, et al. Integrative genomic study of Chinese clear cell renal cell carcinoma reveals features associated with thrombus. Nat Commun. 2020;11:739.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Warsow G, Hubschmann D, Kleinheinz K, Nientiedt C, Heller M, Van Coile L, et al. Genomic features of renal cell carcinoma with venous tumor thrombus. Sci Rep. 2018;8:7477.

    PubMed  PubMed Central  Google Scholar 

  18. Turajlic S, Xu H, Litchfield K, Rowan A, Chambers T, Lopez JI, et al. Tracking cancer evolution reveals constrained routes to metastases: TRACERx renal. Cell. 2018;173:581–594 e512.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med. 2020;26:909–18.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Motzer RJ, Escudier B, McDermott DF, George S, Hammers HJ, Srinivas S, et al. Nivolumab versus Everolimus in advanced renal-cell carcinoma. N Engl J Med. 2015;373:1803–13.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. McDermott DF, Huseni MA, Atkins MB, Motzer RJ, Rini BI, Escudier B, et al. Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Nat Med. 2018;24:749–57.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Young MD, Mitchell TJ, Vieira Braga FA, Tran MGB, Stewart BJ, Ferdinand JR, et al. Single-cell transcriptomes from human kidneys reveal the cellular identity of renal tumors. Science. 2018;361:594–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Labbate C, Hatogai K, Werntz R, Stadler WM, Steinberg GD, Eggener S, et al. Complete response of renal cell carcinoma vena cava tumor thrombus to neoadjuvant immunotherapy. J Immunother Cancer. 2019;7:66.

    PubMed  PubMed Central  Google Scholar 

  24. Shepherd ARH, Joshi R, Tan CP, Cohen P, Brook NR. Inferior vena cava thrombectomy following complete response to nivolumab/ipilimumab for metastatic renal cell carcinoma. ANZ J Surg. 2020;90:1517–9.

    PubMed  Google Scholar 

  25. Bhat A, Nahar B, Venkatramani V, Banerjee I, Kryvenko ON, Parekh DJ. Metastatic renal cell carcinoma with level IV thrombus: contemporary management with complete response to neoadjuvant targeted therapy. Case Rep Urol. 2019;2019:7102504.

    PubMed  PubMed Central  Google Scholar 

  26. Berends J, Gourley E, Kaushik D. Robust response to nivolumab in patient with renal cell carcinoma inferior vena cava tumour thrombus. BMJ Case Rep. 2019;12(4):e227030.

    PubMed  PubMed Central  Google Scholar 

  27. Thommen DS, Schumacher TN. T cell dysfunction in cancer. Cancer Cell. 2018;33:547–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Khan O, Giles JR, McDonald S, Manne S, Ngiow SF, Patel KP, et al. TOX transcriptionally and epigenetically programs CD8(+) T cell exhaustion. Nature. 2019;571:211–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Ge G, Peng D, Xu Z, Guan B, Xin Z, He Q, et al. Restoration of 5-hydroxymethylcytosine by ascorbate blocks kidney tumour growth. EMBO Rep. 2018;19:e45401.

    PubMed  PubMed Central  Google Scholar 

  30. Wingender G, Kronenberg M. OMIP-030: characterization of human T cell subsets via surface markers. Cytometry A. 2015;87:1067–9.

    PubMed  Google Scholar 

  31. Krishna C, DiNatale RG, Kuo F, Srivastava RM, Vuong L, Chowell D, et al. Single-cell sequencing links multiregional immune landscapes and tissue-resident T cells in ccRCC to tumor topology and therapy efficacy. Cancer Cell. 2021;39:662–677 e666.

    CAS  Google Scholar 

  32. van der Putten C, Remmerswaal EBM, Terpstra ML, van der Bom ND, Kers J, Ten Berge IJM, et al. CD8 and CD4 T cell populations in human kidneys. Cells. 2021;10(2):288.

    PubMed  PubMed Central  Google Scholar 

  33. Park JG, Na M, Kim MG, Park SH, Lee HJ, Kim DK, et al. Immune cell composition in normal human kidneys. Sci Rep. 2020;10:15678.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Turner JE, Becker M, Mittrucker HW, Panzer U. Tissue-resident lymphocytes in the kidney. J Am Soc Nephrol. 2018;29:389–99.

    CAS  PubMed  Google Scholar 

  35. Sun Y, Wu L, Zhong Y, Zhou K, Hou Y, Wang Z, et al. Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. Cell. 2021;184:404–421 e416.

    CAS  PubMed  Google Scholar 

  36. Hu J, Chen Z, Bao L, Zhou L, Hou Y, Liu L, et al. Single-cell transcriptome analysis reveals intratumoral heterogeneity in ccRCC, which results in different clinical outcomes. Mol Ther. 2020;28:1658–72.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. Miller BC, Sen DR, Al Abosy R, Bi K, Virkud YV, LaFleur MW, et al. Subsets of exhausted CD8(+) T cells differentially mediate tumor control and respond to checkpoint blockade. Nat Immunol. 2019;20:326–36.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Bi K, He MX, Bakouny Z, Kanodia A, Napolitano S, Wu J, et al. Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma. Cancer Cell. 2021;39:649–661 e645.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Engblom C, Pfirschke C, Pittet MJ. The role of myeloid cells in cancer therapies. Nat Rev Cancer. 2016;16:447–62.

    CAS  PubMed  Google Scholar 

  40. Okabe Y, Medzhitov R. Tissue biology perspective on macrophages. Nat Immunol. 2016;17:9–17.

    CAS  PubMed  Google Scholar 

  41. Murdoch C, Muthana M, Coffelt SB, Lewis CE. The role of myeloid cells in the promotion of tumour angiogenesis. Nat Rev Cancer. 2008;8:618–31.

    CAS  PubMed  Google Scholar 

  42. Deligne C, Midwood KS. Macrophages and extracellular matrix in breast cancer: partners in crime or protective allies? Front Oncol. 2021;11:620773.

    PubMed  PubMed Central  Google Scholar 

  43. Kovaleva OV, Samoilova DV, Shitova MS, Gratchev A. Tumor associated macrophages in kidney cancer. Anal Cell Pathol (Amst). 2016;2016:9307549.

    Google Scholar 

  44. Maynard A, McCoach CE, Rotow JK, Harris L, Haderk F, Kerr DL, et al. Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing. Cell. 2020;182:1232–1251 e1222.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Wettersten HI, Aboud OA, Lara PN Jr, Weiss RH. Metabolic reprogramming in clear cell renal cell carcinoma. Nat Rev Nephrol. 2017;13:410–9.

    CAS  PubMed  Google Scholar 

  46. Yu L, Lu M, Jia D, Ma J, Ben-Jacob E, Levine H, et al. Modeling the genetic regulation of cancer metabolism: interplay between glycolysis and oxidative phosphorylation. Cancer Res. 2017;77:1564–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Jia D, Lu M, Jung KH, Park JH, Yu L, Onuchic JN, et al. Elucidating cancer metabolic plasticity by coupling gene regulation with metabolic pathways. Proc Natl Acad Sci U S A. 2019;116:3909–18.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. Menon R, Otto EA, Hoover P, Eddy S, Mariani L, Godfrey B, et al. Single cell transcriptomics identifies focal segmental glomerulosclerosis remission endothelial biomarker. JCI Insight. 2020;5(6):e133267.

    PubMed Central  Google Scholar 

  49. Goveia J, Rohlenova K, Taverna F, Treps L, Conradi LC, Pircher A, et al. An integrated gene expression landscape profiling approach to identify lung tumor endothelial cell heterogeneity and angiogenic candidates. Cancer Cell. 2020;37:421.

    CAS  PubMed  Google Scholar 

  50. Xing X, Yang F, Huang Q, Guo H, Li J, Qiu M, et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv. 2021;7(5):eabd9738.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. Benjamin LE, Hemo I, Keshet E. A plasticity window for blood vessel remodelling is defined by pericyte coverage of the preformed endothelial network and is regulated by PDGF-B and VEGF. Development. 1998;125:1591–8.

    CAS  PubMed  Google Scholar 

  52. Raza A, Franklin MJ, Dudek AZ. Pericytes and vessel maturation during tumor angiogenesis and metastasis. Am J Hematol. 2010;85:593–8.

    CAS  PubMed  Google Scholar 

  53. Escudier B, Eisen T, Stadler WM, Szczylik C, Oudard S, Siebels M, et al. Sorafenib in advanced clear-cell renal-cell carcinoma. N Engl J Med. 2007;356:125–34.

    CAS  PubMed  Google Scholar 

  54. Cella D, Escudier B, Rini B, Chen C, Bhattacharyya H, Tarazi J, et al. Patient-reported outcomes for axitinib vs sorafenib in metastatic renal cell carcinoma: phase III (AXIS) trial. Br J Cancer. 2013;108:1571–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Sternberg CN, Davis ID, Mardiak J, Szczylik C, Lee E, Wagstaff J, et al. Pazopanib in locally advanced or metastatic renal cell carcinoma: results of a randomized phase III trial. J Clin Oncol. 2010;28:1061–8.

    CAS  PubMed  Google Scholar 

  56. Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352:189–96.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. Jin SZ, Li RY, Chen MY, Yu C, Tang LQ, Liu YM, et al. Single-cell transcriptomic analysis defines the interplay between tumor cells, viral infection, and the microenvironment in nasopharyngeal carcinoma. Cell Res. 2020;30:950–65.

    CAS  PubMed  PubMed Central  Google Scholar 

  58. Weaver VM, Petersen OW, Wang F, Larabell CA, Briand P, Damsky C, et al. Reversion of the malignant phenotype of human breast cells in three-dimensional culture and in vivo by integrin blocking antibodies. J Cell Biol. 1997;137:231–45.

    CAS  PubMed  PubMed Central  Google Scholar 

  59. He X, Lee B, Jiang Y. Cell-ECM interactions in tumor invasion. Adv Exp Med Biol. 2016;936:73–91.

    CAS  PubMed  Google Scholar 

  60. Kim K, Zhou Q, Christie A, Stevens C, Ma Y, Onabolu O, et al. Determinants of renal cell carcinoma invasion and metastatic competence. Nat Commun. 2021;12:5760.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. Sade-Feldman M, Yizhak K, Bjorgaard SL, Ray JP, de Boer CG, Jenkins RW, et al. Defining T cell states associated with response to checkpoint immunotherapy in melanoma. Cell. 2019;176:404.

    CAS  PubMed  PubMed Central  Google Scholar 

  62. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8:329–337 e324.

    CAS  PubMed  PubMed Central  Google Scholar 

  63. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36:411–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. Stewart BJ, Ferdinand JR, Young MD, Mitchell TJ, Loudon KW, Riding AM, et al. Spatiotemporal immune zonation of the human kidney. Science. 2019;365:1461–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. Liao M, Liu Y, Yuan J, Wen Y, Xu G, Zhao J, et al. Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19. Nat Med. 2020;26:842–4.

    CAS  PubMed  Google Scholar 

  66. Borcherding N, Vishwakarma A, Voigt AP, Bellizzi A, Kaplan J, Nepple K, et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol. 2021;4:122.

    CAS  PubMed  PubMed Central  Google Scholar 

  67. He Q, Mok TN, Yun L, He C, Li J, Pan J. Single-cell RNA sequencing analysis of human kidney reveals the presence of ACE2 receptor: a potential pathway of COVID-19 infection. Mol Genet Genomic Med. 2020;8:e1442.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Melms JC, Biermann J, Huang H, Wang Y, Nair A, Tagore S, et al. A molecular single-cell lung atlas of lethal COVID-19. Nature. 2021;595:114–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. Chen Z, Zhou L, Liu L, Hou Y, Xiong M, Yang Y, et al. Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma. Nat Commun. 2020;11:5077.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. Wang L, He T, Liu J, Tai J, Wang B, Chen Z, et al. Pan-cancer analysis reveals tumor-associated macrophage communication in the tumor microenvironment. Exp Hematol Oncol. 2021;10:31.

    CAS  PubMed  PubMed Central  Google Scholar 

  71. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019;20:163–72.

    CAS  PubMed  PubMed Central  Google Scholar 

  72. Finak G, McDavid A, Yajima M, Deng J, Gersuk V, Shalek AK, et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 2015;16:278.

    PubMed  PubMed Central  Google Scholar 

  73. Kim N, Kim HK, Lee K, Hong Y, Cho JH, Choi JW, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun. 2020;11:2285.

    CAS  PubMed  PubMed Central  Google Scholar 

  74. Ma L, Hernandez MO, Zhao Y, Mehta M, Tran B, Kelly M, et al. Tumor cell biodiversity drives microenvironmental reprogramming in liver cancer. Cancer Cell. 2019;36:418–430 e416.

    CAS  PubMed  PubMed Central  Google Scholar 

  75. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  76. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739–40.

    CAS  PubMed  PubMed Central  Google Scholar 

  77. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    PubMed  PubMed Central  Google Scholar 

  79. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.

    PubMed  PubMed Central  Google Scholar 

  80. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32:381–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  81. Shi Y, Zhang Q, Bi H, Lu M, Tan YZ, Zou DJ, Ge LY, Chen ZG, Liu C, Ci WM, Ma LL. Decoding the multicellular ecosystem of vena caval tumor thrombus in clear cell renal cell carcinoma by single-cell RNA sequencing. NGDC. 2022. https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000963.

  82. Wang XM, Lu Y, Song YM, Dong J, Li RY, Wang GL, Wang X, Zhang SD, Dong ZH, Lu M, et al. Integrative genomic study of Chinese clear cell renal cell carcinoma reveals features associated with thrombus. NCBI SRA. 2020. https://www.ncbi.nlm.nih.gov/bioproject/PRJNA596338.

  83. McDermott DF, Huseni MA, Atkins MB, Motzer RJ, Rini BI, Escudier B, Fong L, Joseph RW, Pal SK, Reeves JA, et al. Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Eur Genome Phenome Arch. 2018. https://ega-archive.org/studies/EGAS00001002928.

Download references

Acknowledgements

We are grateful to Prof. Yongliang Zhao, Fai Bai, Yungui Yang and Yang Liang for their thoughtful and valuable comments and discussions.

Review history

The review history is available as Additional file 8.

Peer review information

Anahita Bishop and Stephanie McClelland were the primary editors of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Funding

This work was supported by the CAS Strategic Priority Research Program (XDA16010102 to W.C.), the National Key R&D Program of China (2018YFC2000100 and 2019YFA0110900 to W.C.), the National Natural Science Foundation of China (31800695 to Y.S., 82173061 to W.C., 81772731 to K.C, 82102975 to H.B., and 82070778 to C.L.).

Author information

Affiliations

Authors

Contributions

CL, WC and LM conceived the project. YS, QZ, HB, ML, ZC, LG, YT and DZ performed the experiments. YS and QZ conducted the bioinformatics analyses. YS, QZ and WC wrote the manuscript with help from all the authors. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Cheng Liu, Weimin Ci or Lulin Ma.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Peking University Third Hospital Medical Science Research Ethics Committees (M2017147), and was conducted in accordance with the Declaration of Helsinki. All the patients in this study have provided written informed consents.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1.

Supplementary Figures S1-S6 and corresponding legends.

Additional file 2: Table S1

. Clinical variables for samples and patients.

Additional file 3: Table S2.

Differential expression results for all lymphoid clusters, differential expression results for CD8+ T cell clusters, CD8-C1 related gene signatures, and progenitor/terminally exhausted signature scores for all CD8+ T cells, related to Fig. 2.

Additional file 4: Table S3.

Differential expression results for all myeloid cell clusters, for PT and TT from all TAMs clusters, and published gene signatures for GSEA, related to Fig. 3.

Additional file 5: Table S4.

Differential expression results for all tumor cell clusters, for PT and TT from all tumor clusters, related to Fig. 4.

Additional file 6: Table S5.

Differential expression results for all endothelial clusters and all myofibroblast clusters, and published gene signatures for GSEA, related to Fig. 5.

Additional file 7: Table S6.

Ligand-receptor means and p-values for all cell interactions in ART, PT and TT (from CellPhoneDB v2.0) and created gene signatures, related to Fig. 6.

Additional file 8.

Review history.

Rights and permissions

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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Shi, Y., Zhang, Q., Bi, H. et al. Decoding the multicellular ecosystem of vena caval tumor thrombus in clear cell renal cell carcinoma by single-cell RNA sequencing. Genome Biol 23, 87 (2022). https://doi.org/10.1186/s13059-022-02651-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13059-022-02651-9

Keywords

  • Tumor thrombus
  • Clear cell renal cell carcinoma
  • Single-cell RNA sequencing
  • Tumor heterogeneity
  • Tumor microenvironment