Participant characteristics

Between August 24, 2021 and May 12, 2023, 36 participants were assessed for eligibility, of whom 28 were enrolled into this phase 2 trial (CASTLE-01, NCT05010668) (Fig. 1a,b). All eligible participants were histologically diagnosed with locally advanced or metastatic ICC and had radiological evidence of progressive disease at enrollment. All enrolled participants had to have at least two distinct lesions, with one intrahepatic lesion amenable to cryoablation and the remaining lesion(s), either intrahepatic or extrahepatic, designated as the target lesion(s) for response assessment (Fig. 1c). Baseline characteristics for the study population are shown in Table 1 and Supplementary Table 2. The average age of the male (61%) and female (39%) population was 62.6 ± 6.8 years. The Eastern Cooperative Oncology Group (ECOG) scores were 0 for 17 participants (61%), 1 for eight participants (29%) and 2 for the remaining three participants (11%). Five participants (18%) had hepatitis B infection and three (11%) had liver cirrhosis (all graded as Child–Pugh class A). A total of 22 participants (79%) had tumor–node–metastasis (TNM) stage IV disease and two (7%) had stage IIIB disease, while four participants (14%) had postoperative recurrence.

Fig. 1: Trial design and participant outcomes.

a, Overall trial design. qd, once daily; i.v., intravenous injection; p.o, oral administration; q3w, once every 3 weeks. b, Trial screening, enrollment, treatment and discontinuation profile. Participants were considered to have completed treatment if they finished 24 months on lenva + sint. c, Lesion selection for cryoablation versus response evaluation. d, Duration of treatment exposure and efficacy assessments. PD, progressed disease. e, Best percentage change in the sum of the longest diameter of the target lesion from baseline (tumor burden). Each bar represents one participant. The ORR and DCR are listed. *CR was confirmed by the reduction of all pathological lymph nodes to a short-axis diameter < 10 mm, which were the only target lesions in this participant. f, Changes in investigator-assessed target tumor burden over time. Tumor responses were measured in accordance with RESICT version 1.1 every 6 weeks after the initiation of immunotherapy. Each line represents one participant. g, Kaplan–Meier survival curves for PFS and OS. m, median. h, Representative images of three participants with significantly prolonged survival (PR for P02; CR for P03 and P09) at baseline and best response. For d–f, 28 participants were included.

Source data

Table 1 Baseline characteristics of participantsClinical efficacy

All 28 (100%) participants had lesions that were evaluable for response, which excluded cryoablation-treated areas (Fig. 1c). As described in the protocol, the primary endpoint was confirmed ORR; secondary efficacy endpoints included disease control rate (DCR), duration of response (DOR), PFS and OS according to Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1. Efficacy data are listed in Supplementary Table 3. The investigator-assessed ORR was 75% (95% confidence interval (CI): 59–91%), thereby achieving the prespecified primary endpoint. Of these participants, two (7.1%) achieved a confirmed complete response (CR), 19 (67.9%) achieved a confirmed partial response (PR) and seven (25%) achieved stable disease (SD). The DCR was 100% (95% CI: 100–100%). Figure 1 shows the duration of objective responses (Fig. 1d), overall change in disease burden (Fig. 1e) and change in target tumor burden over time (Fig. 1f). Detailed responses in individual participants are shown in Supplementary Table 4. The median time to response (TTR) was 3.6 months (95% CI: 2.6–6.2) and the median DOR was 17.0 months (95% CI: 7.9–not reached (NR)). As of the data cutoff on March 2, 2024, PFS events had occurred in 15 participants (54%), among whom three were still alive. With a median follow-up time of 16.9 months (95% CI: 14.6–NR), the median PFS and OS were 16.8 months (95% CI: 11.5–NR) and 25.4 months (95% CI: 13.3–NR), respectively (Fig. 1g). Typical images of the 28 participants at baseline and best response after treatment are provided in Fig. 1h and Extended Data Fig. 1.

Safety

To evaluate the safety endpoint, treatment-related adverse events (TRAEs) were collected. TRAEs of any grade occurred in all 28 (100%) participants, with the most common events being fatigue (68%), hypertension (65%), hypoalbuminemia (57%) and hypothyroidism (57%) (Supplementary Table 5). Grade (G)3 TRAEs occurred in ten (36%) participants and included hypertension (three, 11%), hyperbilirubinemia (three, 11%) and decreased white blood cell count (two, 7%). Other G3 TRAEs included fatigue, increased alanine aminotransferase, palmar–plantar erythrodysesthesia syndrome (PPES), increased γ-glutamyl transferase and postoperative hemorrhage (each in one participant). Hyperbilirubinemia, a G4 TRAE, occurred in one (4%) participant and led to treatment discontinuation. Most other TRAEs were reversible by a reduction in lenvatinib (lenva), dose interruption or administration of antihypertensive drugs. Immune-mediated AEs (irAEs) are summarized in Supplementary Table 5 and were mostly G1–G2 (G3–G4 irAEs occurred in four (14%) participants). Systemic corticosteroid use for irAEs was reported in three (11%) participants, attributed to G3 fatigue, G3 PPES and G2 PPES (Supplementary Table 5). No treatment-related deaths occurred.

Dynamic landscape of the tumor immune microenvironment during CASTLE therapy

Exploratory post hoc analyses of baseline and on-treatment liver biopsies and peripheral blood mononuclear cells (PBMCs) were performed to capture dynamic changes in cellular composition, functional status and the localization of cell populations during CASTLE therapy. In total, 58 liver tumor biopsies and 51 PBMC samples were collected from 26 participants at three time points: at baseline, 2 weeks after cryoablation and after two cycles of lenva and sintilimab (lenva + sint) (Fig. 2a and Supplementary Table 6). Response to CASTLE therapy did not correlate with the actionable alterations, including IDH1 mutations, MET and ERBB2 amplifications, tumor mutational burden (TMB) and PDL1 expression, but correlated with PFS (19.37 months for CR and PR versus 9.77 months for SD; P = 0.0088) and OS (25.9 months for CR and PR versus 11.5 months for SD; P = 0.0012) (Fig. 2b–e). Bulk RNA transcriptomics of liver biopsies (n = 22) revealed a distinct pattern of immune cell infiltration at baseline (Fig. 2f). Among these, 11 (50%) participants were clustered on the basis of a ‘hot’ tumor phenotype characterized by increased infiltration of immune cells (especially tumor-infiltrating T cells), upregulation of inhibitory checkpoints, activation of interferon (IFN) responses and increased production of proinflammatory chemokines. The remaining 11 (50%) participants comprised a ‘cold’ phenotype cluster characterized by scarce or absent T cell infiltration and defects in the IFN pathways. Surprisingly, these tumor phenotypes did not significantly correlate with treatment response. Compared to cold-phenotype participants, those with the hot phenotype did not have a prolonged PFS (16.8 months for hot versus 11.5 months for cold; P = 0.84) or OS (25.4 months for hot versus 11.5 months for cold; P = 0.17) (Fig. 2g,h). Longitudinal analysis before and after treatment showed that the CASTLE regimen had the potential to turn cold tumors into inflamed, hot tumors (Fig. 2f).

Fig. 2: Overview of the genomic and transcriptomic landscape of the CASTLE-01 cohort.

a, Experimental design of multimodel sequencing, comprising data from tumor tissues and PBMCs collected from the CASTLE-01 cohort of participants with ICC at three time points: before treatment (pre), after cryoablation (cryo) and after lenva + sint (detailed sample collection in Supplementary Table 6). PD, progressed disease. b, Mutational landscape in the cohort identified by WES obtained from 25 baseline samples. c, MSI status, TMB and PDL1 expression (as determined by IHC) from 28 participants at baseline. d, Kaplan–Meier survival curves for PFS and OS in the cohort, analyzed by the log-rank test. e, Distribution of treatment responses (CR, PR and SD) by status of biomarkers PDL1 (IHC), TMB, MET, ERBB2, IDH1 and BRAF in pretreatment samples, analyzed by two-sided Fisher’s exact test. Participant counts are displayed within the grid. f, Heat map of gene set variation analysis scores of representative pathways and expression of selected genes across 46 bulk RNA-seq samples (pre, n = 22; cryo, n = 13; lenva + sint, n = 11) collected from 23 participants at three time points. Hierarchical clustering separated the TME of pretreatment samples into hot and cold phenotypes. Paired, PDL1 (IHC) high-expression and high-TMB samples are labeled. g, Kaplan–Meier survival curves of PFS and OS in the cohort, grouped by TME phenotype, analyzed by the log-rank test. h, Distribution of treatment responses by TME phenotype in pretreatment samples, analyzed by two-sided Fisher’s exact test. Participant counts are displayed within the grid. i, UMAP plots of six major cell classes captured across 45 scRNA-seq samples: CD8+ T cells, CD4+ T cells, myeloid cells, EPs, mesenchymal cells and ECs. Colors indicate 48 different cell types. j, Proportional changes in cell types across 45 scRNA-seq samples compared pairwise across three treatment time points (each time point, n = 15), analyzed by a paired two-sided Wilcoxon test. Colors indicate the cell type assignment. Red dashed lines represent the P-value threshold of 0.05. FC, fold change.

Source data

Further single-cell gene expression profiles of 45 paired samples before treatment, after cryoablation and after lenva + sint treatment were obtained to create a single-cell atlas (Fig. 2a and Supplementary Table 6). After filtering low-quality cells and applying batch correction using the R package Harmony, we obtained single-cell RNA (scRNA) profiles from 267,018 cells, identifying immune and nonimmune cell clusters shared across participants and treatments (Extended Data Fig. 2a,b). The major cell types and subclusters were identified using tissue-type-specific canonical marker genes in the scRNA sequencing (scRNA-seq) data (Fig. 2i and Extended Data Fig. 2c,d). Because all participants achieved tumor control (DCR = 100%) and participants with SD achieved a median PFS of 9.77 months, greatly surpassing the historical response to monotherapy or combination therapy with ICIs, we focused on the cellular changes occurring at three treatment time points without comparing the CR + PR and SD groups (Extended Data Fig. 2a,e,f). We also applied imaging mass cytometry (IMC) to gain a spatial-resolution-based phenotyping of the tumor microenvironment (TME) and to validate the scRNA profile results (Fig. 2a). We acquired images from 88 regions of interest (ROIs) derived from 15 paired samples and based on the expression of 30 proteins (Supplementary Table 6), yielding data on 88,924 cells (Extended Data Fig. 2g–i). First, we observed dramatic increases in CD8+ T cells, CD4+ T cells and B cells and a decrease in malignant epithelial cells (EPs) in both scRNA and IMC data (Extended Data Fig. 2a,e,f,j,k). Longitudinal analysis showed that cryoablation induced recruitment of CD8+ effector T cells (Teffs) to the TME, which was boosted by subsequent treatment with lenva + sint (Fig. 2j). Other changes included increased abundance of CD4+ memory T cells following cryoablation. Treatment with lenva + sint significantly decreased the abundances of stalk and tip endothelial cells (ECs), CD4+ T helper (Th)17 cells and regulatory T (Treg) cells, while increasing the abundances of CD8+PD1hi Teffs, CD4+CXCL13+ T follicular helper (Tfh) cells, RELN+ cancer-associated fibroblasts (CAFs) and postcapillary venule (PCV) ECs (Fig. 2j). Taken together, our results suggested that the CASTLE regimen reshaped the TME, modulating nonimmunoreactive tumors to become responsive to immunotherapy.

Cryoablation-triggered continuous replenishment of intratumoral PD1hi CD8 Teffs from peripheral blood

We stratified subsets of tumor-infiltrated CD8+ T cells on the basis of expression profiles of pertinent markers (Fig. 3a and Extended Data Fig. 3a). In total, seven discrete CD8+ T cell groups were identified, with four displaying increased PD1 expression (subsets T1–T4) and increased immune response to tumor cells (Fig. 3b and Extended Data Fig. 3b). T1 expressed terminal exhaustion T cell (Tex) markers (LAG3, HAVCR2 and CTLA4), T2 displayed increased levels of proliferation markers (STMN1 and MKI67), T3 demonstrated augmented expression of IFN-stimulated gene modules (MX1 and ISG15) and T4 showcased elevated levels of effector molecules (GZMK and EOMES). The subset of CD8+PD1hi Teffs in T4 was identified as the precursor population of T1 CD8+PD1hi Tex, as supported by RNA velocity analysis (Extended Data Fig. 3c–e). Longitudinal analysis showed that cryoablation-induced CD8+PD1hi Teff recruitment to tumors was boosted after lenva + sint treatment (Fig. 3c and Extended Data Fig. 3f), which was confirmed by IMC data (Fig. 3d,e and Extended Data Fig. 3g,h). Additionally, we observed increased progenitor-like and effector phenotypes and decreased Tex phenotype in CD8+PD1hi T cells (Fig. 3f,g). Furthermore, the abundance of clonal CD8+ T cells (clone size ≥ 2) in the CD8+PD1hi Teff cell subset increased after cryoablation and lenva + sint therapy (Extended Data Fig. 3i).

Fig. 3: Cryoablation-triggered continuous replenishment of intratumoral CD8+PD1hi Teffs from peripheral blood.

a, Heat map showing standardized average expression values of cluster-defining genes in CD8+ T cell clusters on population level. b, Spearman correlation map showing associations of the proportions of CD8+ T cell subtypes among 15 pretreatment samples in the scRNA-seq data. c, Box plot of scRNA-seq data showing the intraparticipant dynamic changes in the frequency of CD8+PD1hi T cell subtypes in all CD8+ T cells of paired samples stratified by three time points: pre, cryo and lenva + sint (each time point, n = 12). d, Representative IMC images showing the distribution of CD8+PD1hiGZMK+ Teffs in participant P02 at the three time points. Scale bars, 100 µm and 50 µm. e, Box plot of IMC data showing the dynamic changes in proportion of CD8+PD1hi Teffs in all CD8+ T cells at the inter-ROI level stratified by the three time points (pre, n = 32; cryo, n = 27; lenva + sint, n = 29; obtained from five participants each time point). f, Heat map displaying standardized average expression values in typical genes of T cell-related modules of CD8+PD1hi T cells at the three time points on population level. g, Violin plots showing scores for exhaustion, progenitor, proliferation and effector or memory signatures of CD8+PD1hi T cells at the three time points. h, Intraparticipant clonal expansion levels (left) and clonal transition levels (right) of CD8+PD1hi Teffs in paired samples at the three time points, quantified using the STARTRAC index. i, Alluvial diagram displaying CD8+PD1hi clonotype-mediated cluster associations. A flow or association among the three time points represents all pairs of cells sharing the same TCR clonotype but may transit to different phenotypes during treatment. j,k, Violin plot (j) revealing the change in clone size of persistent CD8+PD1hi Teff clones at the three time points at a TCR level (each time point, n = 40 clones) and UMAP plot (k) displaying the distribution of these clones over time (representing as blue dots). l, UMAP distribution of novel clones of CD8+PD1hi Teffs at the three time points (representing as blue dots). m, Pie charts showing the percentage of novel clones of CD8+PD1hi Teffs among the total preexisting cells (left) and TCR clones (right) at the lenva + sint time point. n, Clone size of novel and preexisting clones of CD8+PD1hi Teffs at the lenva + sint time point at a TCR level. o, Dot plots showing the expression levels of selected genes in novel and preexisting clones at the lenva + sint time point. Dot size indicates the percentage of cells in which gene expression was detected; color indicates the mean expression level. p,q, Frequencies of the TCR nucleotide sequences of novel clones identified in peripheral blood bulk TCR samples at the three time points illustrated by a heat map (p) and a box plot (q) (each time point, n = 24 clones). For a,f,g, data were obtained from 45 scRNA-seq samples at three time points (each time point, n = 15). For h–l, data were obtained from 27 scTCR-seq samples at three time points (each time point, n = 9). For m–o, data were analyzed in nine scTCR-seq samples after lenva + sint treatment. For c,e,g,h,j,q, the center line indicates the median value, bottom and top hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5× the interquartile range. For c,h,j,q, a paired two-sided Wilcoxon test was applied. For e,g,n, an unpaired two-sided Wilcoxon test was applied.

Source data

Next, we discovered that cryoablation may promote CD8+ Teff infiltration in two ways by tracing the T cell receptor (TCR) expansion and migration patterns. First, cryoablation induced expansion of preexisting tumor-infiltrating CD8+ Teffs, as evidenced by increased clonal expansion and transition scores of CD8+PD1hi Teffs (Fig. 3h,i), and increased number and clone size of shared TCRs between baseline and after cryoablation, which was further boosted after treatment with lenva + sint (Fig. 3j,k). Second, we observed a notable emergence of novel clones of tumor-infiltrating CD8+ Teffs after treatment with lenva + sint that were absent at the previous time points (Fig. 3l). These novel clones accounted for 45.6% of cell numbers and 66.8% of the TCR clones of T4 CD8+PD1hi Teffs in tumors after treatment with lenva + sint (Fig. 3m). They also exhibited smaller clone size and higher expression of progenitor, proliferation and migration markers (Fig. 3n,o) compared to preexisting CD8+ Teffs, suggesting that they recently infiltrated into the TME with a progenitor-like, effector phenotype. Tracking of the emerging novel CD8+PD1hi Teff clones in the TME after lenva + sint treatment in the peripheral bulk TCR-seq data revealed that, among the 24 TCRs detectable in peripheral blood, 23 (95.8%) emerged in circulation after cryoablation (Fig. 3p,q), indicating that the newly emerged TCRs infiltrated the TME after lenva + sint treatment.

We also identified seven distinct phenotypic clusters of CD4+ T cells in the scRNA-seq data (Extended Data Fig. 3j,k), among which the T8 cluster exhibited high expression of PD1 and CXCL13, characteristic of Tfh cells. T8, defined as CD4+CXCL13+ Tfh cells, displayed high immune response to tumor cells (Extended Data Fig. 3k) and proportionally increased after lenva + sint therapy (Extended Data Fig. 3l). Positive correlations at a sample level of CD4+CXCL13+ Tfh cells with CD8+PD1hi Teffs, B cells and plasma cells were identified (Extended Data Fig. 3m). Examining the immune infiltration based on the hematoxylin and eosin (H&E) stainings also confirmed that both intratumoral and stromal lymphocytes exhibited a significant increase following lenva + sint treatment (Extended Data Fig. 3n–p). Furthermore, after treatment, the formation of tertiary lymphoid structures (TLSs) was observed in P02, who achieved PR (Extended Data Fig. 3q). The signature of TLS was also upregulated after CASTLE treatment in the bulk RNA-seq data (Extended Data Fig. 3r). These findings collectively indicate that the CASTLE regimen had the potential to boost the CD8+ and CD4+ T cells infiltration.

Cryoablation-increased tumor immunogenicity

Cluster analysis of isolated tumor cells revealed a total of 13 distinct EP subgroups (Extended Data Fig. 4a). All subgroups were confirmed to exhibit malignancy using the inferCNV program (Extended Data Fig. 4b). There was upregulation of antigen presentation and the IFNγ and IFNα modules after cryoablation, accompanied by downregulation of the transforming growth factor (TGF)β module after lenva + sint therapy (Fig. 4a and Extended Data Fig. 4c–e). IMC displayed greater abundance of HLAI+ EPs at the two later time points (Fig. 4b,c and Extended Data Fig. 4f). Given the participant specificity of EP9–EP13, we sought to focus on the malignant subgroups EP1–EP8 (Extended Data Fig. 4g). Cryoablation or lenva + sint therapy did not alter the proportions of these subgroups within tumors (Extended Data Fig. 4h). Using the bulk transcriptome data of our cohort, we identified genes related to poor PFS or OS (cutoffs: hazard ratio (HR) > 1 and P < 0.05). Interestingly, projection of the two negative gene sets onto the uniform manifold approximation and projection (UMAP) embedding revealed significant enrichment in the EP1 cluster (Fig. 4d,e and Extended Data Fig. 4i,j). Compared to the other malignant EP subgroups, EP1 showed upregulation in glycolysis pathways, alongside downregulation in immune activation-related pathways (Fig. 4f). In the context of developmental lineages, EP1 served as the progenitor of other tumor cells (Extended Data Fig. 4k). A notable surge in immune-related gene expression within EP1 subsequent to interventions with cryoablation and lenva + sint therapy was also observed (Fig. 4g and Extended Data Fig. 4l,m). Additionally, immune-related modules including IFNγ, IFNα and antigen presentation were markedly upregulated in the EP1 subcluster (Extended Data Fig. 4n,o). Using the top differentially expressed genes in the EP1 subgroup, EP1 scores of pretreatment bulk transcriptomic samples were calculated (Supplementary Table 7). Notably, bulk RNA samples with higher EP1 scores exhibited elevated expression of immune inhibitory signaling pathways and greater infiltration of immune inhibitory cellular components (Fig. 4h and Extended Data Fig. 4p,q). Conversely, bulk RNA samples with lower EP1 scores demonstrated abundant expression of immune-activated pathways and TME components. Additionally, immune responses after cryoablation also differed between EP1low and EP1high groups (Extended Data Fig. 4r). At the single-cell transcriptomic level, a significant negative correlation was observed between the proportions of EP1 and CD4+CXCL13+ Tfh and CD8+PD1hi Teff cells (Extended Data Fig. 4s,t). Thus, the baseline EP1 score demonstrated robust predictive value for PFS and OS after CASTLE treatment (Fig. 4i).

Fig. 4: Cryoablation increased tumor immunogenicity and boosted DC activation.

a, Violin plots showing the signature scores of Kyoto Encyclopedia of Genes and Genomes (KEGG) antigen processing and presenting, Hallmark IFNγ response, IFNα response and TGFβ signaling in EPs at three time points. b, Representative IMC images showing the distribution of HLAI+ EPs of participant P03 at the three time points. Scale bars, 100 µm and 50 µm. c, Box plot of IMC data showing the dynamic changes in proportion of HLA+EPs in all EPs at the inter-ROI level stratified by the three time points from five participants. d,e, UMAP plots of EP clusters colored by gene signatures identified in bulk RNA data related to poor PFS (d) and poor OS (e). f, Comparative gene set enrichment analysis of pretreatment gene sets in the EP1 cluster and other malignant EP clusters in 15 pretreatment scRNA-seq samples. NES, normalized enrichment score. P values were determined using a one-tailed permutation test. g, Heat map displaying standardized average expression values in typical genes of immune-related modules in EP1 cells at the three time points on population level. h, Heat map showing abundance of pathways and TME components per sample in pretreatment bulk RNA-seq samples stratified by median EP1 signature scores into high-EP1 and low-EP1 signature groups. Significantly differentially expressed pathways are indicated (unpaired two-sided t-test, *P < 0.05). i, Kaplan–Meier survival curves for PFS (left) and OS (right), grouped by EP1 signature scores of pretreatment bulk RNA-seq samples, analyzed by the log-rank test. j, Volcano plots showing the genes with differential expression between the time points before treatment and after cryoablation in all DCs in the scRNA-seq data (each time point n = 15). The x axis and y axis values were calculated by the Seurat method. Blue and green dots indicate genes with significant upregulation at the time points after cryoablation and before treatment, respectively. Significant genes of interest are labeled. k, Dot plots showing the expression levels of selected genes in all DCs at three time points. Dot size indicates the percentage of cells in which the gene was detected. Color indicates the mean expression. l, Numbers of enriched genes (count) of enriched pathways in all DCs at the time points before treatment and after cryoablation (each time point, n = 15) in the scRNA-seq data, analyzed by the hypergeometric test. m, Representative IMC images showing the distribution of CD11c+HLA-DR+IDO1+ APCs of participant P21 at the three time points. Scale bars, 100 µm and 50 µm. n, Quantification of the percentages of APCs in myeloid cells in IMC data at the inter-ROI level stratified by the three time points. o, Box plot of IMC data showing the dynamic changes in IDO1 expression in APCs at the inter-ROI level, stratified by the three time points (pre, n = 17; cryo, n = 20; lenva + sint, n = 23; obtained from five participants each time point); ROIs without IDO1+ APCs were excluded from analysis. p, Interaction score between CD8+PD1+ Teffs (middle) and APCs (neighbor) in IMC at the inter-ROI level, stratified by the three time points. Black bars indicate the mean and 95% CI. q, Quantification of the distance of CD8+PD1+ Teffs to the nearest APC in IMC, stratified by the three time points. For a,d,e,g,k, data were obtained from 45 scRNA-seq samples at three time points (each time point, n = 15). For a,c,j,n–q, an unpaired two-sided Wilcoxon test was applied. For c,n,p,q, ROIs obtained from five participants each time point (pre, n = 32; cryo, n = 27; lenva + sint, n = 29) were analyzed. For a,c,n,o,q, the center line indicates the median value, bottom and top hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5× the interquartile range.

Source data

Cryoablation-boosted dendritic cell activation and antigen presentation

Dendritic cells (DCs) were categorized into three distinct subgroups: DC1 (CD1C), DC2 (CLEC9A) and DC3 (LAMP3) (Extended Data Fig. 5a,b). The lack of change in the proportions of the DC group and subgroups after cryoablation (Extended Data Fig. 5c,d), likely because of limited DC capture numbers, prompted us to explore functional differences in the subgroups. Representing traditional antigen-presenting cells (APCs), DC1 and DC2 exhibited high expression levels of CD74 and major histocompatibility complex class II (MHC II) genes (Extended Data Fig. 5e) and enrichment of antigen presentation pathways (Extended Data Fig. 5f). By contrast, DC3 exhibited elevated expression of costimulatory and tertiary lymphoid markers (CD80, CD86, CD40, CCL19 and CCR7) (Extended Data Fig. 5e) and involvement in T cell immunoregulation and migration (Extended Data Fig. 5f). DCs showed upregulation in the expression of antigen presentation and costimulatory genes after cryoablation (Fig. 4j,k), as well as enrichment in immune processes related to T cell activation and DC differentiation (Fig. 4l and Extended Data Fig. 5g,h). After cryoablation, DC1 and DC2 showed further upregulation of the antigen-presentation-related molecules (Extended Data Fig. 5i,j), while DC3 showed enhanced expression of the costimulatory molecules (Extended Data Fig. 5k). In terms of their interaction with T cells, all three subsets of DCs exhibited an upregulated ability to present antigen to T4 CD8+PD1hi Teffs (Extended Data Fig. 5l). Within the receptor–ligand network, molecules facilitating T cell migration and maturation were also elevated, thus fostering an environment conducive to subsequent immunotherapy (Extended Data Fig. 5o). Consistently, IMC showed a significant increase in the frequency of CD11c+HLA-DR+IDO1+ APCs after cryoablation (Fig. 4m,n, and Extended Data Fig. 5m,n), along with upregulated expression of IDO1 (Fig. 4o), indicative of the IFN response29. Furthermore, spatial interaction between APCs and CD8+PD1hi Teffs was stronger after cryoablation, as evidenced by higher interaction scores and closer distance at spatial resolution (Fig. 4p,q).

While neutrophils may potentially form neutrophil extracellular traps (NETs) after ablation to create an immunosuppressive TME30, our comprehensive analysis of neutrophil changes following cryoablation revealed a significant reduction in neutrophil proportions along with downregulation of neutrophil activation and chemotaxis-related pathways (Extended Data Fig. 5c,p,q). Multiplex immunofluorescence further confirmed decreased NET formation (Extended Data Fig. 5r,s). Importantly, we observed concurrent upregulation of antigen presentation genes and immune-activating pathways in the scRNA-seq data (Extended Data Fig. 5p,q).

Intratumoral lymphocyte influx boosted by targeting tumor vasculature

In the aforementioned experiments, we noticed that many novel CD8+PD1hi Teff clones entered circulation after cryoablation. However, these clones only infiltrated tumors after lenva + sint therapy. We speculated that lenva may normalize the tumor vasculature, thus boosting intratumoral lymphocyte influx31,32. Therefore, we annotated EC subtypes with markers identified from published datasets, resulting in eight EC clusters33,34 (Fig. 5a and Extended Data Fig. 6a). Traditional angiogenic cell phenotypes, such as tip cells and stalk cells, were significantly reduced following lenva treatment, possibly because of their elevated expression of vascular endothelial growth factor receptor (VEGFR; KDR and FLT1)33 (Fig. 5b). Concurrently, a noticeable increase in PCV ECs, also known as high endothelial venules, which expressed high levels of adhesion markers (SELE and SELP), were observed after lenva34 (Fig. 5b). Pseudotime analysis unveiled sequential transcriptional alterations during EC development, with tip and stalk cells at the terminus of the trajectory, indicative of their status in immature neovascularization (Fig. 5c,d and Extended Data Fig. 6b). Pathway analysis indicated that PCV ECs exhibited enrichment in modules associated with proinflammatory effects and lymphocyte adhesion, whereas tip and stalk cells showed enrichment in pathways related to tumor promotion, angiogenesis and extracellular matrix (ECM) regulation (Fig. 5d,e and Extended Data Fig. 6c). The use of lenva upregulated genes related to cytokine ligands, cell adhesion and antigen presentation, while downregulating TGFβ-related genes and lenva-targeting genes (Extended Data Fig. 6d). These observations were verified in pathway analysis (Extended Data Fig. 6e). Specifically, lenva decreased the expression of angiogenesis and ECM modules in stalk and tip cells (Fig. 5f,g and Extended Data Fig. 6f) while increasing immune-related functions, including antigen presentation and cell adhesion, in PCV cells (Fig. 5h,i and Extended Data Fig. 6g). Furthermore, the frequency of CD8+PD1hi Teff and CD4+CXCL13+ Tfh cells was positively correlated with PCV ECs but negatively correlated with stalk and tip cells (Fig. 5j and Extended Data Fig. 6h), suggesting an association between vasculature normalization and lymphocyte influx.

Fig. 5: Boosting intratumoral lymphocyte influx through targeting of the tumor vasculature.

a, Heat map showing standardized average expression values of cluster-defining genes in eight EC clusters on population level. b, Box plot showing dynamic proportional changes in all EC subtypes in paired samples at three time points of scRNA-seq data (each time point, n = 11). Data were analyzed using a paired two-sided Wilcoxon test. c, Pseudotime-ordered analysis of EC clusters EC1–EC4 in scRNA-seq data. Dot colors indicate pseudotime inferred by Monocle2 (top), cell type (middle) and time point (bottom). d, Heat map of the dynamic changes in gene expression along the pseudotime of the EC1–EC4 transition (divided into three phases), showing the distribution of the four EC subtypes over pseudotime. Top: subclusters are labeled by color. Genes showing differential expression during pseudotime (left) and enriched pathways (right) are labeled. e, Heat map showing scaled mean activities of selected modules in the eight EC clusters on population level. GO, Gene Ontology. f–i, Dynamic changes in transcriptional characteristics in EC subclusters at the three time points. Heat maps showing scaled average activities of selected modules in stalk and tip cells (f) and PCV cells (h) on population level. Dot plots showing the expression levels of selected genes in stalk and tip cells (g) and PCV cells (i). Dot size indicates the percentage of cells in which the gene was detected. Color indicates the mean expression. j, Scatter plots showing significant correlations between the proportion of CD8+PD1hi Teffs in CD8+ T cells and the percentages of stalk cells (top), tip cells (middle) and PCV cells (bottom) in all ECs. Data analyzed by the two-sided Pearson correlation test (top) or the two-sided Spearman correlation test (middle and bottom). k, Scatter plots showing significant correlations between the proportion of CD8+PD1hi Teffs in CD8+ T cells of scRNA-seq data and the TVNS in bulk RNA-seq samples (pre, n = 12;