Skip to main content

Spatial transcriptomics analysis of neoadjuvant cabozantinib and nivolumab in advanced hepatocellular carcinoma identifies independent mechanisms of resistance and recurrence



Novel immunotherapy combination therapies have improved outcomes for patients with hepatocellular carcinoma (HCC), but responses are limited to a subset of patients. Little is known about the inter- and intra-tumor heterogeneity in cellular signaling networks within the HCC tumor microenvironment (TME) that underlie responses to modern systemic therapy.


We applied spatial transcriptomics (ST) profiling to characterize the tumor microenvironment in HCC resection specimens from a prospective clinical trial of neoadjuvant cabozantinib, a multi-tyrosine kinase inhibitor that primarily blocks VEGF, and nivolumab, a PD-1 inhibitor in which 5 out of 15 patients were found to have a pathologic response at the time of resection.


ST profiling demonstrated that the TME of responding tumors was enriched for immune cells and cancer-associated fibroblasts (CAF) with pro-inflammatory signaling relative to the non-responders. The enriched cancer-immune interactions in responding tumors are characterized by activation of the PAX5 module, a known regulator of B cell maturation, which colocalized with spots with increased B cell marker expression suggesting strong activity of these cells. HCC-CAF interactions were also enriched in the responding tumors and were associated with extracellular matrix (ECM) remodeling as there was high activation of FOS and JUN in CAFs adjacent to the tumor. The ECM remodeling is consistent with proliferative fibrosis in association with immune-mediated tumor regression. Among the patients with major pathologic responses, a single patient experienced early HCC recurrence. ST analysis of this clinical outlier demonstrated marked tumor heterogeneity, with a distinctive immune-poor tumor region that resembles the non-responding TME across patients and was characterized by HCC-CAF interactions and expression of cancer stem cell markers, potentially mediating early tumor immune escape and recurrence in this patient.


These data show that responses to modern systemic therapy in HCC are associated with distinctive molecular and cellular landscapes and provide new targets to enhance and prolong responses to systemic therapy in HCC.


Hepatocellular carcinoma (HCC) is one of the most common causes of cancer-associated deaths globally and the most rapidly rising cause of cancer death in the USA [1,2,3]. Immune checkpoint inhibitors (ICIs) that target programmed cell death protein-1 (PD1) have modest clinical activity as monotherapy in HCC but may be more effective in combination with other therapeutic agents, including anti-angiogenic therapies. The combination of bevacizumab (an anti-VEGF antibody) plus atezolizumab (an ICI targeting the PD1 axis) was recently established as a preferred first-line standard of care for patients with unresectable HCC, and clinical efficacy has also been reported with multiple other anti-angiogenic/immune checkpoint inhibitor (ICI) combinations (cabozantinib + atezolizumab and apatinib + camrelizumab) [4,5,6,7,8,9,10]. While anti-angiogenic plus ICI combinations have demonstrated clinical benefit, a significant proportion of patients does not respond to such therapies, and biomarkers to predict response and determine which patients will truly benefit from the treatment are currently not available [11, 12].

Window of opportunity studies provide an unparalleled opportunity to elucidate the mechanisms of response and resistance to systemic therapy, yielding abundant tissue for deep interrogation of the tumor microenvironment (TME) in patients receiving neoadjuvant systemic therapy. We recently reported the results of a window of opportunity clinical trial of angiogenic/ICI therapy in which HCC patients were treated with 8 weeks of cabozantinib and nivolumab (CABO/NIVO) followed by attempted surgical resection (NCT03299946) [10]. In this clinical study, 5 out of 15 patients achieved pathologic responses, with outstanding disease-free survival noted among patients with a pathologic response. Spatial proteomics profiling of the HCC surgical specimens identified an enriched immune effector infiltrate, and reduced immunosuppressive macrophages, among patients with pathologic response to CABO/NIVO [10]. Nevertheless, molecular pathways that underlie these immune interactions as well as the tumor intrinsic mechanisms of resistance to ICIs in HCC are not clear. In this current study, we performed a high-dimensional, unbiased profiling of patients enrolled in our window of opportunity clinical trial of cabozantinib and nivolumab in HCC to uncover the intrinsic molecular and cellular tumor features of the tumor microenvironment that underlie the differential responses to systemic therapy.

The recent development of technologies that provide spatially resolved gene expression data introduced powerful methods to profile the TME and understand how tumor intrinsic features are associated with the distribution of other crucial cell types for tumor development and response to therapies [13]. Using spatial transcriptomics (ST), it is possible to examine the samples’ molecular and cellular compositions and interactions among the different cellular components as it maintains the tissue architecture [14]. Thus, ST provides spatial visualization of immune cell distribution in relation to the cancer cells and the correlation with the molecular profile that drive the infiltration of certain immune types. Overall, it is an exceptional experimental approach to understand the responses and resistance to immunotherapies as the molecular mechanisms driving or inhibiting effective anti-tumor response can now be examined within the cellular context.

In this study, we used ST to investigate the tumor-related features that explain the response or resistance to CABO/NIVO. The ST profiling of responder and non-responder samples identified transcriptional signatures that are associated with immune activation and metabolism, respectively. The intercellular interaction analysis shows that in responders these cancer-immune communications lead to B cell activation, while in the areas of cancer-cancer associated fibroblasts (CAF) interactions there is activation of extracellular matrix (ECM) remodeling genes. The molecular and cellular analysis with ST indicates that the mechanisms of response depend on the infiltration of immune cells into the HCC TME combined with the ability of cancer cells to trigger the immune response. On the other hand, non-responder samples are poor in immune cells because of unsuccessful antigen presentation. The ST analysis also identified remarkable HCC heterogeneity in one sample with a tumor cluster resembling the responders’ samples and another cluster resembling a non-responder tumor that is potentially associated with disease recurrence due to the expression of cancer stem cell (CSC) features not observed in other samples from our cohort. Our study demonstrates the power of ST analysis to understand the responses to therapies in the context of TME composition and molecular features and identifies independent mechanisms of resistance and recurrence in HCC.


Patients and sample acquisition

Clinical outcomes and detailed study methods for our clinical trial of CABO/NIVO were recently described by our group [10]. The trial was prospectively registered at as NTC03299946 (, and the institutional review board of Johns Hopkins University approved the protocol. All enrolled patients provided written informed consent. Briefly, we enrolled 15 patients with potentially resectable HCC on a single-arm, open-label, phase 1 clinical trial of neoadjuvant cabozantinib and nivolumab. The enrollment period was from 2018–05-14 to 2019–12-12. The primary outcome was the safety and feasibility of this combination therapy. The secondary outcomes are the percentage of patients eligible for resection, percentage of participants with complete pathologic response, percentage of participants with major pathologic responses, objective response rate, median overall survival (OS), and disease-free survival. Patients were enrolled at the Liver Cancer Multidisciplinary Clinic at the Johns Hopkins Sidney Kimmel Comprehensive Cancer Center in Baltimore. The key eligibility criteria included borderline resectable or locally advanced HCC, age greater than or equal to 18 years, an Eastern Cooperative Oncology Group performance status score of 0 or 1, and preserved liver function with a Child–Pugh score of A. The enrolled patient population included patients with high-risk tumor features that historically predict poor outcomes with upfront surgical resection such as multinodular disease, portal vein invasion, or large tumors. Cabozantinib and nivolumab were provided by Exelixis and Bristol Myers Squibb, respectively.

Of the 15 patients enrolled in the clinical trial, 12 patients achieved successful R0 resection. Of these 12 patients, 5 achieved a major or complete pathological response (defined as 90% or greater tumor necrosis in the surgical resection specimen) and were classified as pathologic responders. The remaining 7 patients were considered non-responders [10]. Tumor samples from all surgical resections were immediately submerged in the optimal cutting temperature (OCT) compound and immersed in liquid nitrogen for quick freezing. All OCT frozen samples were stored at − 80 °C until use.

Spatial transcriptomics data generation

To prepare the ST slides, the samples preserved in OCT were sectioned, stained with hematoxylin and eosin, and examined by an experienced pathologist (R.A.). The areas for ST analysis were chosen based on tumor viability, presence of stroma with immune infiltration, and cancer-associated fibroblasts when possible. For each sample, a 6 × 6 mm region with those characteristics was selected for sectioning and mounting onto the ST slides.

The ST data was generated using the commercial platform Visium (10× Genomics). Briefly, from each surgical sample, a 5-µm section was placed in the designated area at the Visium slide and immediately stored at − 80 °C until use. The sections were fixed in cold methanol for 30 min at − 20 °C. The fixed samples were stained with hematoxylin and eosin (H&E) and imaged using the Nanozoomer scanner (Hamamatsu) at 40x magnification. Samples were permeabilized for 30 min at 37 °C with the permeabilization enzyme provided with the Visium Spatial Gene Expression Reagent Kits (10× Genomics). Following permeabilization, reverse transcription; cDNA second strand synthesis, denaturation, and amplification; and library construction were performed according to the manufacturer’s instructions. All libraries were sequenced with a depth of at least 50,000 reads per spot (minimum of ~ 250 millions per sample) at the NovaSeq (Illumina).

Spatial transcriptomics data analysis

Sequencing data was processed using the Space Ranger software (10× Genomics) for demultiplexing and FASTQ conversion of barcodes and reads data, alignment of barcodes to the stained tissue image, and generation of read count matrices. The processed sequencing data were inputs for the analyses using the Seurat (version 4.1.1) [15]. Data preprocessing with Seurat involved the initial visualization of the counts onto the tissue image to discriminate technical variance from histological variance (e.g., collagen-enriched regions present lower cellularity that reflects in low counts), removal of necrotic regions and tissue folding spots, and normalization with SCTransform [16]. Following filtering, all samples were merged into a single dataset, and the embedding was corrected for batch effects among samples using Harmony [17]. The data dimensionality was reduced with PCA and then clustered with Leiden [18]. Cell types to cluster assignments were performed based on the most variable features (genes) and revised by a pathologist (R.A.) to confirm with the liver histological regions.

Differential expression analysis of tumor spots from responders versus non-responders was performed using pseudo-bulking on DESeq2 [19]. Genes were considered differentially expressed when log two fold change was at least ± 2 and p-adjusted value ≤ 0.01. The gene set enrichment analysis was performed using the genes ranked according to the log-fold change (LFC) determined by the differential expression analysis. We evaluated in each of the groups, non-responders and responders, the enriched expression of genes belonging to the MSigDB Hallmark Pathways [20].

Transcriptional factor regulatory module activity analysis

Inference of regulatory modules between transcription factors (TF) and downstream regulated genes was performed using SCENIC (version 0.11.2) [21]. The annotation of cis-regulatory motif and genome ranking database are acquired from the cisTarget database available at GRNBoost2 was used to fit the putative regulatory module based on the co-expression between every expressed gene and predefined TF list. Each TF gene with an importance measure higher than 95 percentiles is included in the module. Modules with < 20 genes are excluded from the downstream analysis. RcisTarget identifies TF-binding motif enrichments and prune non-direct binding modules. Finally, regulatory module activities in every spot of Visium data are quantified by AUCell. The dependencies cisTarget, GRNBoost2, and AUCell are wrapped into SCENIC.

Ligand-receptor cell–cell signaling network reconstruction

We used Domino (version 0.1.1) [22] to analyze the signaling networks based on the gene regulatory module activities in the Visium spots. Prior to calculating signaling relationship, we subset clusters to define spatial regions for analysis. Genes that expressed less than 2.5% of spots are excluded from our analysis. Similar to the implementation of Domino for single-cell RNA-seq, we adapted this approach to spatial transcriptomics computing the spot-based Pearson correlation between regulatory modules and normalized, z-score scaled expression of receptors obtained from CellphoneDB (version 2.0) [23]. Correlations are set to zero if the receptor is targeted by the transcription factor. TF regulatory modules are signaled by receptors with a Pearson correlation coefficient larger than 0.3. Cognate ligands identified by CellphoneDB are required for all receptors to be included in the signaling network.

TCGA and CIBERSORT analysis

Gene-level RNA sequencing (RNAseq) data was downloaded from Genomic Data Commons harmonized database for The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) using the TCGAbiolinks package (v2.26.0) [24]. As for gene expression, we used gene counts from STAR alignment that were log2-tranformed for further analysis. We also used CIBERSORT-predicted cell proportions for TCGA-LIHC (n = 371) from Thorsson et al. [25]. We used Spearman correlation to find the association between gene expression and T cell proportions estimated by CIBERSORT. The analysis was done using R/Bioconductor computational environment (v4.2.2). The heatmap was plotted using ComplexHeatmap (v2.14.0) [26]


Spatial transcriptomics identifies HCC cell composition differences between responders and non-responders to immunotherapy

We profiled HCC samples obtained from patients on a recently reported clinical trial of neoadjuvant CABO/NIVO using ST to examine tumor intrinsic molecular and cellular features of response and resistance to the treatment. Among the 12 resection specimens analyzed, 5 had a major or complete pathologic response. ST was performed for all 12 frozen surgical HCC specimens (Fig. 1A), of which 7 (4 responders and 3 non-responders) passed pre-determined quality control parameters (Fig. 1B, C). For 5 samples, due to extensive necrosis of the tumor, the sequencing data did not pass the quality control (low number of counts and genes detected per ST spot) for the analysis. The remaining 7 samples with adequate quality were carried forward for downstream analysis. Unsupervised clustering of the gene expression data obtained from the merged data after Harmony batch correction [17] recapitulates the sample architecture. Gene expression clusters map to their respective major cell types and cell subtypes commonly present in HCC samples: cancer cells, cancer-associated fibroblasts (CAFs), and immune cells (Fig. 1B, C). The marker genes identified during the clustering analysis are established markers for these cell types (Additional file 1: Fig. S1). The cell types assigned to the spatial gene expression clusters were confirmed by a pathologist (R.A.).

Fig. 1
figure 1

Spatial transcriptomics analysis of HCC samples treated with cabozantinib and nivolumab. A Experimental workflow. B Hematoxylin and eosin (H&E)-stained images of the samples profiled and the spatial clusters for responders. C H&E and spatial clusters for non-responders. D Tumor-, immune-, and cancer-associated fibroblast (CAF) composition in each responder sample as determined by spatial transcriptomics. E Tumor-, immune-, and CAF proportions in non-responders samples. To quantify the cell proportions of tumor, immune, and CAF in each patient, the multiple clusters that were classified as the same type were merged into one major category

The assignment of gene expression profiles to cell type composition enables the comparison of cellular proportions between responders and non-responders. To compare cell proportions, we grouped the identified types into three major cell types: tumor, immune, and CAF. These broad cell subtypes were obtained for each patient by collapsing subclusters of these cell types on a per-patient basis. This analysis demonstrated that responders have an increased presence of immune cells while non-responders present with a higher abundance of cancer cells (Fig. 1D, E). This observation corroborates the previous spatial and single-cell proteomics profiling of these samples that showed enrichment for cancer cells in non-responders and higher infiltration by immune cells in responders [10, 27]. In addition, we were able to map the CAF clusters in all 4 responders and one non-responder. The prevalence of HCC cells among non-responders suggests that the lack of response to neoadjuvant CABO/NIVO therapy is tightly correlated with the absence of immune cells infiltrating the tumor.

Gene expression analysis of tumor compartments shows activation of immune-related pathways in responders and of cell growth pathways in non-responders.

To investigate the tumor intrinsic molecular changes associated with response and resistance to the neoadjuvant treatment, we performed differential expression analysis on the ST clusters mapping to HCC cells only. As the ST approach provides genome-wide information, it is a suitable technology to discover the distinct transcriptional signatures between cancer cells from responders and non-responders. Additionally, the ability to map these signatures to the tissue architecture allows the selection of the ST spots that map to HCC cells for a controlled analysis of the cancer components in each sample.

The clusters mapping to tumor areas (Fig. 2A, B) were extracted from the ST data and pseudo-bulked for the differential expression analysis between responders and non-responders. A total of 508 genes are upregulated in responders (Fig. 2C, red dots), and 47 genes are upregulated in non-responders (Fig. 2C, blue dots). The top differentially expressed genes in responders are immune-related genes (e.g., CCL19, CXCL14, IGHM, CXCL6), while in the non-responders, there is increased expression of tumor markers (e.g., AFP, IGF2, WNK4) (Additional file 2: Table S1), suggesting that in patients responding to CABO/NIVO, there is an active immune or inflammatory response not observed in samples from non-responders.

Fig. 2
figure 2

Differential expression analysis of A tumor clusters (subset of spots classified as a tumor in the clustering) from responders versus B tumor clusters from non-responders across all patients. C Volcano plot of the differential expression analysis showing the most differentially expressed genes in responders (red dots) and non-responders (blue dots), showing the upregulation of immune genes in responders relative to the upregulation of hepatocellular markers among non-responders. D Pathway enrichment analysis between responders (red dots) and non-responders (blue dots) reveals activation of immune-related pathways in responders’ tumors, while non-responders' tumors have activation of proliferation and metabolic pathways

Subsequently, gene set enrichment analysis was performed to identify the pathways enriched in the responders’ HCC samples versus non-responders. Samples from patients that responded to CABO/NIVO are enriched for the expression of genes that belong to the pathways associated with active immune response (HALLMARK_INFLAMMATORY_RESPONSE, HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_ALLOGRAFT_REJECTION, HALLMARK_COMPLEMENT) and elevated antigen process and presentation machinery genes (Fig. 2D), suggesting that among responders, tumor intrinsic transcriptional features are leading to the increased immune cell infiltration observed in these samples relative to non-responders samples. This hypothesis is supported by the upregulation of genes involved in antigen processing and presentation (KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION) among responders relative to non-responders (Additional file 3: Fig. S2). On the other hand, non-responders lack the expression of the immune-related pathways and are transcriptionally enriched by genes from signaling the pathways that are involved in maintaining cell proliferation (HALLMARK_E2F_TARGETS, HALLMARK_MYC_TARGETS_V1) and metabolism (HALLMARK_OXIDATIVE_PHOSPHORYLATION, HALLMARK_CHOLESTEROL_HOMEOSTASIS), suggesting that non-responders’ cancer cells have their growth features maintained and activated (Fig. 2D).

Overall, the gene expression analysis suggests that in HCC samples that did not respond to CABO/NIVO, the cancer cells are not affected by the treatment and are able to maintain their proliferative and metabolic functions. On the other hand, the samples from patients that responded to the therapy show that in their cancer cells, these functions are overcome by the activation of the immune-related pathways that could explain the tumor immune infiltration and tumor shrinkage in response to the ICI component of the treatment.

Intercellular interaction analyses identify transcription factor regulatory networks associated with response and resistance to CABO/NIVO

Intercellular interaction analyses are facilitated for ST datasets since tissue architecture is maintained simplifying the selection of neighboring cell types to examine the potential cell-to-cell communication. The active interaction between neighboring cell types is critical to understand cancer biology and response to therapies since the TME has a great influence on tumor biology, progression, and response to therapies. Examining these cellular relationships and the molecular outcomes is critical to determine the potential pathways of intervention for alternative therapeutic options with increased efficacy. The interaction analyses were performed to understand the cell-to-cell crosstalk of HCC-immune and HCC-CAF cells in response to the CABO/NIVO neoadjuvant treatment. For each of the samples, we used the identified spatial clusters and used the Domino software to determine the signaling pathways that are activated because of the intercellular interactions. Due to HCC intrinsic high heterogeneity, interaction analysis was performed for each patient individually.

Among the responders, we observed the activation of the PAX5 network in the immune regions adjacent (spatial gene expression clusters in direct contact) (Additional file 4: Fig. S3) to tumor clusters (Fig. 3A, Additional file 5: Fig. S4), while FOS and JUN modules are highly active in CAFs surrounding the tumor spots (Fig. 3B, C, Additional file 5: Fig. S4). PAX5 is a transcription factor that is central to B cell differentiation [28, 29]. The PAX5 activity, which was investigated in immune-enriched or CAF-enriched regions only, co-localizes with spots showing high expression of the genes CD19, CD22, CD79A, and CIITA (B cell markers) (Fig. 3D), relative to spots corresponding to HCC or CAF clusters confirming that this is an essential factor for B cell lineage activation and maturation. Moreover, it suggests that B cells are a critical component of the tumor immune response in HCC activated by the CABO/NIVO neoadjuvant therapy. The activation of FOS and JUN from the HCC-CAF interactions is related with high expression of ECM remodeling markers (COL1A1, COL3A1, VIM) (Fig. 3E).

Fig. 3
figure 3

Intercellular interaction analysis. A HCC-immune interaction analysis identified the activation of the PAX5 network. B, C HCC-CAF interaction analysis pointed to the activation of the FOS and JUN networks. D PAX5 is a transcription factor that regulates B cell activity, and the PAX5 network identified co-localizes with the distribution of B cells as determined by the spatial distribution of B cell markers. E FOS and JUN are transcription factors that can regulate genes involved in extracellular matrix remodeling, and the networks regulated by these genes colocalize with CAF marker genes

This spatial distribution of tumor immune response and ECM remodeling networks suggests that in the presence of immune cells, there is an active response mostly driven by B cells that could be the initial trigger for tumor cell killing by cytotoxic immune cells and recruitment of other effector immune cells. The presence of the stroma and active remodeling with increased collagen production (COL1A1, COL3A1) could be a result of fibrosis as a response to cell death that recruit CAFs and initiate immune exclusion in collagen/CAF-rich regions and so creates a niche characterized by a lack of immune cells. These findings suggest that drugs that initiate or maintain B cells activity combined with CAF inhibitors could be alternatives to increase the efficacy of immunotherapies to treat HCC.

Spatial transcriptomics analysis reveals specific HCC heterogeneity related with recurrence after neoadjuvant therapy

Among the 5 patients that responded to CABO/NIVO neoadjuvant therapy, four patients remain without disease recurrence at least 3 years after surgery, whereas a single patient developed recurrent disease after 1 year of therapy. We utilized ST to investigate the distinct features of this single clinical outlier (hereby named HCC1-R) that might explain the unexpected early recurrence observed in this patient. The initial pathology examination of this sample identified striking histological intratumoral heterogeneity, unlike any other sample in the cohort. Two distinct tumor regions were apparent in this sample: one that is immune rich (Fig. 4A, cyan) and another that is immune poor (Fig. 4A, dark blue). The histological distinction is recapitulated by the spatially resolved gene expression profile that identifies these two tumor regions by distinct ST clusters (Fig. 4B). As highlighted previously, the major advantage of ST is that it provides genome-wide gene expression while maintaining tissue architecture, thus an ideal approach to analyze such a unique sample.

Fig. 4
figure 4

Intra-sample heterogeneity analysis with spatial transcriptomics. A Responder sample with remarkable heterogeneity with two tumor regions, one immune-rich (cyan) and one immune-poor (dark blue). B Each of these regions presents distinct transcriptional profiles as determined by clustering analysis. CE The immune-rich tumor region highly expresses immune-related pathways that were initially observed to be enriched across responders’ samples. FH Proliferation and metabolic pathways, which are enriched across non-responders’ tumors, are expressed in high levels at the immune-poor region. I-J The immune-rich tumor region expresses high levels of the antigen processing and presentation machinery genes, which is associated with the more efficient attraction of immune cells

We performed differential expression analysis to compare the immune-rich and immune-poor regions to identify the differences that could explain the distinct immune infiltration between these two tumor regions. The analysis shows that among the top upregulated genes from the immune-rich region, the majority are immune markers while the immune-poor region expresses increased levels of HCC-specific tumor markers (Additional file 6: Fig. S5 and Additional file 7: Table S2). The gene set enrichment analysis, similar to the analyses of all responders versus non-responders tumors, reveals the enrichment for immune-related pathways in the tumor region that is highly immune infiltrated and for growth-related pathways in the tumor area that is not infiltrated by immune cells (Additional file 6: Fig. S5). The expression and pathway analysis suggest that the immune-rich region from the sample collected from patient HCC1-R recapitulates the transcriptional profile observed across all responders’ tumors while the immune-poor area is similar to non-responders cancer cells.

To test our hypothesis that the two distinct regions resemble responder (immune rich) and non-responder (immune poor) samples, we examined if the expression of pathways enriched across responders and non-responders is reproducible on the immune-rich and immune-poor regions, respectively, of patient HCC1-R. Using a module score analysis of the signaling pathways enriched in both patient groups (responders and non-responders), we mapped the averaged expression of the set of genes belonging to each of the enriched pathways. Thus, we determined if a pathway is enriched in a spot if the module score is high. The significantly enriched MSigDB pathways in responders are highly expressed in the immune-rich tumor region of patient HCC1-R (Fig. 4C–E and Additional file 8: Fig. S6A, B, C, D, and E), while those enriched in patients from non-responders are upregulated in the immune-poor tumor region (Fig. 4F–H and Additional file 8: Fig. S6F, G, H, I, and J), thus suggesting that the HCC heterogeneity in this sample recapitulates the features of tumors from responders (immune rich) and non-responders (immune poor). To verify if immune infiltration was associated with more efficient antigen processing and presentation, we also examined the module score for the expression of the antigen processing and presentation machinery (KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION). At the immune-rich tumor region, the module score for this pathway is significantly higher when compared to the module score at the immune-poor region (Fig. 4I), suggesting that one of the mechanisms driving the immune cell infiltration is effective antigen processing and presentation. Finally, the intercellular interaction analysis recapitulates what was observed in other responders. PAX5 module is more active in the immune-rich region, whereas FOS and JUN modules have significantly higher activation in CAFs (immune-poor region) (Fig. 5A–C; Additional file 9: Fig. S7). The activation of PAX5 is again associated with the areas enriched for B cells, but the ECM activity is not frequent in those regions (Fig. 5D, E).

Fig. 5
figure 5

Intercellular interaction analysis in the context of intra-sample heterogeneity. A Intercellular interaction analysis in the immune-rich tumor region (cancer-immune) reveals activation of PAX5. B, C FOS and JUN are the active networks from the interaction analysis at the immune-poor tumor region (HCC-CAF). D PAX5 activation co-localizes with the expression of B cell markers concentrated at the borders of the immune-rich tumor region. E FOS and JUN networks are co-expressed with CAF markers adjacent to the immune-poor tumor region

To understand the lack of immune infiltration into one of the tumor regions from patient HCC1-R, we looked at signatures associated with immune evasion. Interestingly, a signature associated with immune evasion and resistance to different types of cancer therapies, including immunotherapies, and that is highly expressed by the immune poor tumor region is from cancer stem cells (CSC) (Fig. 6A). This CSC signature comprises previously published genes (SP, PROM1, ALDH1A1, THY1, EPCAM, ANPEP, CD44, CD24, CXCL12, ICAM1, CACNA2D1, CD47, LGR5, KRT19, ABCG2, AFP, SOX9, POU5F1, SOX2, NANOG, NOTCH, ZIC2, PBX3, ZNF148, TCIM, SALL4, ZFP42, YY1A1, FOXM1, MYCN, SCD, MRPS5, PMPCB, ANGPTL4, PDK4, XDH, IRAK1, PTPN11, ANXA3, IL6, IL8, OSM, IGF1, FGF2, ANGPTL1, CTSS, UBE2T, EPHB2) and was compiled in a recent review [30]. CSCs are associated with therapeutic resistance, immune escape, recurrence, and metastasis. The presence of HCC cells with stemness features has been previously observed in HCC and also associated with therapeutic resistance [30].

Fig. 6
figure 6

Cancer stem cell detection in hepatocellular carcinoma heterogeneous sample. A Cancer stem cell markers are highly expressed at the immune-poor tumor region of the responder sample with intrasample heterogeneity as noted by the spatial expression and the levels of expression are significantly different between the distinct tumor regions. B Heatmap with the expression of the CSC genes in the TCGA HCC samples. C Correlation heatmap with the Spearman correlation coefficient between expression of CSC markers and T cell proportions (CIBERSORT) in HCC samples from TCGA showing a negative correlation between T CD8 cell proportions and high expression of some CSC genes

To validate that the CSC molecular signature is associated with HCC that fails to elicit an effective immune response and so is more prone to resistance to ICIs, we verified the expression of the CSC genes in HCC samples from The Cancer Genome Atlas (TCGA) (Fig. 6B). Using CIBERSORT, we obtained the T cell proportions from the same samples. Finally, we examined the Spearman correlation between the cell proportions from CIBERSORT with the expression of the genes in the CSC signature (Fig. 6C). Overall, the high expression of a significant number of CSC markers in HCC tumors (Fig. 6B) was correlated with low proportions of detected T CD8+ cells (negative correlation) (Fig. 6C), a cell type that is associated with better outcomes and prolonged response to ICIs (Fig. 6B). This finding corroborates the hypothesis that CSCs evade the immune system resulting in lower infiltration of effector T cells. The presence of CSC in the only patient that recurred from the CABO/NIVO neoadjuvant therapy after pathological response suggests that the presence of these cells is associated with therapeutic resistance and that this population will then persist and lead to recurrence.


Modern systemic treatment for HCC combining anti-VEGF therapy with ICI therapy prolongs survival and is widely used in the first-line treatment of patients with HCC [5]; however, little is known about the mechanisms of response and resistance to this new treatment paradigm. Prior analysis from our group using samples from a clinical trial of neoadjuvant CABO/NIVO demonstrates that clinical responses are associated with higher infiltration of immune cells while resistance was associated with diminished presence of immune cells that were in its majority immunosuppressive [10, 27]. However, due to the limited marker panels in the experimental approaches in that previous study, the molecular pathways that underlie immune response and changes in the cancer cells themselves were not investigated to determine their role in the response to neoadjuvant therapy. ST analysis is a novel approach that enables whole transcriptomics analysis within the tissue architecture context. We applied ST to unveil the cancer cells’ intrinsic features that are associated with response and resistance to CABO/NIVO. Since the ST approach used to profile the specimens collected is genome-wide, the analysis is not restricted to markers or pathways previously selected.

Using ST, we identified gene expression clusters that map to specific cell types in the HCC samples. With that, the gene expression analysis of the subset of cancer cells from each patient demonstrates that in responders, there is an active anti-tumor immune response characterized by the high expression of immune-related genes and recruitment of immune-related pathways. Moreover, in non-responders’ cancer cells, we observe the activation of proliferation and metabolic pathways. The intercellular interaction analysis of cancer-immune cells demonstrates that the immune response observed is driven by a strong B cell response. The cancer-immune cell interaction is associated with the activation of a PAX5 module that is enriched in regions with an increased presence of B cells. PAX5 is a transcription factor essential for B cell maturation, commitment, immunoglobulin rearrangements, and activity [28]. In tumors, B cells are frequently found in tertiary lymphoid structures where they act by presenting antigens to T cells that could lead to their activation into effector cytotoxic cells [31]. Thus, the presence of B cell activation suggests that there is also induction of a T cell cytotoxic response and effective response to treatment with an ICI, in this case, nivolumab. Another observation from the intercellular interaction analysis is that in the regions with HCC-CAF communication, there is an activation of ECM remodeling. This was also an observation across the responders’ samples. The active ECM remodeling in the presence of therapeutic response indicates that in the areas of tumor cell death, there is a fibrosis process that replaces the space once occupied by the cancer cells. Further analysis with a temporal model or sample collection would reveal if that process will lead to future resistance to the therapy, as the presence of desmoplasia with a dense collagen ECM may create a barrier to immune cell infiltration into the remaining tumor. Overall, these findings suggest that therapeutic resistance is associated with the inability of cancer cells to trigger an anti-tumor immune response that would increase immune cell infiltration and that this response relies on the presence of B cells that will prime or activate cytotoxic T cells. Thus, B cell presence is a potential cellular biomarker of response to therapies with ICIs.

With the spatially resolved analysis of this cohort of neoadjuvant-treated HCC samples, it was possible to examine a sample from a patient that initially presented with a therapeutic response but that later recurred. The surgical specimen collected from this patient shows a notable tumor heterogeneity pattern with an immune-poor tumor region separated from an immune-rich tumor region. This tissue distribution would have been lost if the samples were profiled using single-cell approaches as those depend on sample dissociation. The immune-poor tumor region molecularly resembles the tumors from non-responders, while the immune-rich counterpart is similar to the responders’ samples that have high immune activity. The immune infiltration is a result of more effective antigen presentation by the cancer cells from the immune-rich region. Successful antigen presentation is a well-known component of prosperous tumor immune response and could be used as a molecular biomarker to predict responses to ICIs [32]. On the other hand, the immune poor tumor area does not express high levels of the antigen-presenting machinery and thus evades the immune system. One of the reasons that could explain the downregulation of antigen presentation is that the cancer cells in that cluster express CSC markers. Due to their stemness features, these cells are expected to be more successful in immune escape and are also more prone to therapeutic resistance [30, 33]. Furthermore, the presence of CSC could explain the tumor recurrence in this patient since these unique cells are frequently present in disease relapse. Hence, the presence of CSC in the tumor is another biomarker that could potentially predict response to therapies and the risk of recurrence after treatment. Ultimately, with the analysis of this clinical outlier, we were also able to demonstrate a distinction between intrinsic resistance and recurrence. The CSC signature is only upregulated in patient HCC1-R. We did not observe the CSC signature in the samples from patients that did not respond to CABO/NIVO. Thus, the mechanisms of intrinsic resistance and acquired resistance, leading to recurrence, are different.


Although our study is limited to a small sample size, this is an exceptional cohort that used a therapeutic combination to treat HCC with promising clinical benefits to a fraction of patients. Regarding the experimental approach, ST also has limitations as there is no single-cell resolution due to the fact that each spot can capture the signal of more than one cell limiting the ability of specifically identifying individual immune cell types and subtypes. However, it demonstrates its usefulness since it allows heterogeneity identification in the spatial context. As previously mentioned, we would have not been able to identify the unique tumor architecture of the single responder that recurred after therapy. The analysis of that sample allowed the identification of an important tumor feature that is potentially associated with the recurrence in that specific patient. In summary, the ST analysis of HCC clinical trial samples identified molecular and cellular mechanisms of therapeutic response and resistance. The mechanisms unveiled by the spatially resolved gene expression analyses would not be found by bulk or single-cell analysis as the lack of the spatial component would not have permitted the interpretations related to the cellular distributions and interactions we were able to extract. Further analysis of other clinical cohorts is essential to corroborate our findings and could potentially identify alternative therapeutic interventions for cancer types that still do not benefit from treatment with ICIs.

Availability of data and materials

All spatial transcriptomics preprocessed data from this study are available from NCBI GEO GSE238264 ( [34]. The raw data is available under controlled access according to the institution’s policy for personal patient data protection. Raw data access can be obtained by contacting Dr. Luciane T. Kagohara ( All codes for the analyses presented are available at GitHub ( [35].



Alpha fetoprotein




Cancer-associated fibroblast


Chemokine ligand 19


Class II major histocompatibility complex transactivator


Collagen type I alpha 1 chain


Collagen type III alpha 1 chain


Cancer stem cell


Chemokine ligand 6


Chemokine ligand 14


Complementary deoxyribonucleic acid


Cancer stem cells


Extracellular matrix


Fos proto-oncogene


Hepatocellular carcinoma


Hematoxylin and eosin


Immune checkpoint inhibitor


Insulin-like growth factor 2


Immunoglobulin heavy constant mu


Jun proto-oncogene


Log-fold change


Liver hepatocellular carcinoma




Optimal cutting compound


Paired box 5


Principal component analysis


Programmed cell death 1


Spatial transcriptomics


The Cancer Genome Atlas


Transcription factor


Tumor microenvironment


T regulatory cell


Vascular endothelial growth factor




WNK lysine-deficient protein kinase 4


  1. Llovet JM, Zucman-Rossi J, Pikarsky E, Sangro B, Schwartz M, Sherman M, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2016;2:16018.

    Article  PubMed  Google Scholar 

  2. Villanueva A. Hepatocellular carcinoma. N Engl J Med. 2019;380:1450–62.

    Article  CAS  PubMed  Google Scholar 

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

    Article  PubMed  Google Scholar 

  4. Gordan JD, Kennedy EB, Abou-Alfa GK, Beg MS, Brower ST, Gade TP, et al. Systemic therapy for advanced hepatocellular carcinoma: ASCO guideline. J Clin Oncol. 2020;38:4317–45.

    Article  CAS  PubMed  Google Scholar 

  5. Finn RS, Qin S, Ikeda M, Galle PR, Ducreux M, Kim T-Y, et al. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N Engl J Med. 2020;382:1894–905.

    Article  CAS  PubMed  Google Scholar 

  6. Finn RS, Ikeda M, Zhu AX, Sung MW, Baron AD, Kudo M, et al. Phase ib study of lenvatinib plus pembrolizumab in patients with unresectable hepatocellular carcinoma. J Clin Oncol. 2020;38:2960–70.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Kelley RK, Yau T, Cheng AL, Kaseb A, Qin S, Zhu AX, et al. VP10-2021: Cabozantinib (C) plus atezolizumab (A) versus sorafenib (S) as first-line systemic treatment for advanced hepatocellular carcinoma (aHCC): results from the randomized phase III COSMIC-312 trial. Ann Oncol. 2022;33:114–6.

    Article  Google Scholar 

  8. Abou-Alfa GK, Chan SL, Kudo M, Lau G, Kelley RK, Furuse J, et al. Phase 3 randomized, open-label, multicenter study of tremelimumab (T) and durvalumab (D) as first-line therapy in patients (pts) with unresectable hepatocellular carcinoma (uHCC): HIMALAYA. JCO. 2022;40 4_suppl:379–379.

    Article  Google Scholar 

  9. Xu J, Zhang Y, Jia R, Yue C, Chang L, Liu R, et al. Anti-PD-1 antibody SHR-1210 combined with apatinib for advanced hepatocellular carcinoma, gastric, or esophagogastric junction cancer: an open-label, dose escalation and expansion study. Clin Cancer Res. 2019;25:515–23.

    Article  CAS  PubMed  Google Scholar 

  10. Ho WJ, Zhu Q, Durham J, Popovic A, Xavier S, Leatherman J, et al. Neoadjuvant cabozantinib and nivolumab converts locally advanced HCC into resectable disease with enhanced antitumor immunity. Nat Cancer. 2021;2:891–903.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhu AX, Abbas AR, de Galarreta MR, Guan Y, Lu S, Koeppen H, et al. Molecular correlates of clinical response and resistance to atezolizumab in combination with bevacizumab in advanced hepatocellular carcinoma. Nat Med. 2022;28:1599–611.

    Article  CAS  PubMed  Google Scholar 

  12. Aldea M, Andre F, Marabelle A, Dogan S, Barlesi F, Soria J-C. Overcoming resistance to tumor-targeted and immune-targeted therapies. Cancer Discov. 2021;11:874–99.

    Article  CAS  PubMed  Google Scholar 

  13. Rao A, Barkley D, França GS, Yanai I. Exploring tissue architecture using spatial transcriptomics. Nature. 2021;596:211–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Davis-Marcisak EF, Deshpande A, Stein-O’Brien GL, Ho WJ, Laheru D, Jaffee EM, et al. From bench to bedside: single-cell analysis for cancer immunotherapy. Cancer Cell. 2021;39:1062–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184:3573–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 2019;20:296.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16:1289–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep. 2019;9:5233.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14:1083–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Cherry C, Maestas DR, Han J, Andorko JI, Cahan P, Fertig EJ, et al. Computational reconstruction of the signalling networks surrounding implanted biomaterials from single-cell transcriptomics. Nat Biomed Eng. 2021;5:1228–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. Cell PhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. 2020;15:1484–506.

    Article  CAS  PubMed  Google Scholar 

  24. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44:e71.

    Article  CAS  PubMed  Google Scholar 

  25. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang T-H, et al. The immune landscape of cancer. Immunity. 2018;48:812-830.e14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.

    Article  CAS  PubMed  Google Scholar 

  27. Mi H, Ho WJ, Yarchoan M, Popel AS. Multi-scale spatial analysis of the tumor microenvironment reveals features of cabozantinib and nivolumab efficacy in hepatocellular carcinoma. Front Immunol. 2022;13:892250.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Holmes ML, Pridans C, Nutt SL. The regulation of the B-cell gene expression programme by Pax5. Immunol Cell Biol. 2008;86:47–53.

    Article  CAS  PubMed  Google Scholar 

  29. Cobaleda C, Schebesta A, Delogu A, Busslinger M. Pax5: the guardian of B cell identity and function. Nat Immunol. 2007;8:463–70.

    Article  CAS  PubMed  Google Scholar 

  30. Lee TKW, Guan XY, Ma S. Cancer stem cells in hepatocellular carcinoma - from origin to clinical implications. Nat Rev Gastroenterol Hepatol. 2022;19:26–44.

    Article  PubMed  Google Scholar 

  31. Kinker GS, Vitiello GAF, Ferreira WAS, Chaves AS, Cordeiro de Lima VC, Medina T da S. B cell orchestration of anti-tumor immune responses: a matter of cell localization and communication. Front Cell Dev Biol. 2021;9:678127.

  32. Sharma P, Allison JP. Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential. Cell. 2015;161:205–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Yoon SK. The biology of cancer stem cells and its clinical implication in hepatocellular carcinoma. Gut Liver. 2012;6:29–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhang S, Yuan L, Danilova L, Mo G, Zhu Q, Deshpande A, et al. Spatial transcriptomics analysis of neoadjuvant cabozantinib and nivolumab in advanced hepatocellular GEO Accession viewer. Accessed 1 Aug 2023.

  35. Zhang S, Yuan L, Danilova L, Mo G, Zhu Q, Deshpande A, et al. Spatial transcriptomics analysis of neoadjuvant cabozantinib and nivolumab in advanced hepatocellular codes. 2023. Zenodo.

Download references


We thank the Johns Hopkins University School of Medicine Experimental and Computational Genomics Core for performing the sequencing of the spatial transcriptomics libraries and for the pre-processing of the raw sequencing data.


SU2C/AACR DT-14–14 (E.M.J.), the Emerson Cancer Research Fund (E.M.J., E.J.F.), an Allegheny Health Network (AHN) grant (E.J.F.), U01CA212007 (E.J.F. and A.SP.), U01CA253403 (E.J.F.), SPORE GI P50CA062924-24A1 (E.M.J., M.Y., E.J.F., and L.T.K.), R01CA138264 (ASP), Exelixis (M.Y.), Bristol Myers Squibb (M.Y.), P30CA006973, and Johns Hopkins Bloomberg-Kimmel Institute for Cancer Immunotherapy. This work was made possible in part through the support of the Maryland Cancer Moonshot Research Grant to the Johns Hopkins Medical Institutions (FY24).

Author information

Authors and Affiliations



EJF and LTK instigated and supervised the study. SZ, MY, EJF, and LTK planned, designed, and wrote the manuscript with input from all authors. LTK and GM performed the experimental methods. SZ, LY, and LD performed the computational analysis and data interpretation. GM, QZ, AD, ATFB, JE, ASP, and RA provided technical and material support. MY and EMJ provided clinical expertise. All authors approved the final manuscript.

Corresponding authors

Correspondence to Elana J. Fertig or Luciane T. Kagohara.

Ethics declarations

Ethics approval and consent to participate

All patients provided written informed consent prior to enrollment, and the trial was registered under as NCT03299946. The protocol was approved by the Institutional Review Board (IRB) at Johns Hopkins University under the number IRB00149350. All patients provided written informed consent prior to enrollment for the use of clinical and genomic data for research. The study conformed to the principles of the Helsinki Declaration.

Consent for publication

Not applicable.

Competing interests

R.A.A. reports receiving a commercial research support from Bristol-Myers Squibb and is a consultant/advisory board member for Bristol-Myers Squibb, Merck, AstraZeneca, Incyte, and RAPT Therapeutics. E.M.J. reports other support from Abmeta, personal fees from Genocea, personal fees from Achilles, personal fees from DragonFly, personal fees from Candel Therapeutics, other support from the Parker Institute, grants and other support from Lustgarten, personal fees from Carta, grants and other support from Genentech, grants and other support from AstraZeneca, personal fees from NextCure, and grants and other support from Break Through Cancer outside of the submitted work. M.Y. reports receiving research grants from Incyte, Bristol-Myers Squibb, and Exelixis and is a consultant for AstraZeneca, Eisai, Exelixis, and Genentech. E.J.F. is on the Scientific Advisory Board of Viosera Therapeutics/Resistance Bio and is a consultant to Mestag Therapeutics. The remaining authors declare that they have no competing interests. Cabozantinib was supplied by Exelixis, and nivolumab was supplied by Bristol-Myers Squibb.

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: Fig. S1.

Heatmap of cluster markers across samples.

Additional file 2: Table S1.

Results from the differential expression analysis between responders and non-responders.

Additional file 3: Fig. S2.

Overall expression of the antigen processing and presentation signaling pathway in the subset of spots in the spatial transcriptomics data from non-responders and responders.

Additional file 4: Fig. S3.

Representation of adjacent regions selection for interaction analysis in one HCC sample.

Additional file 5: Fig. S4.

Heatmap of transcription factor network scores from patient HCC3-R as a representation of the DOMINO networks obtained from the interaction analysis.

Additional file 6: Fig. S5.

Differential expression analysis between the immune poor and immune rich regions on sample from patient HCC1-R.

Additional file 7: Table S2.

Results from the differential expression analysis between immune poor versus immune rich tumor regions from patient HCC1-R.

Additional file 8: Fig. S6.

Expression of MSigDB Hallmark pathways across all spots on the immune poor and immune rich tumor regions from patient HCC1-R.

Additional file 9: Fig. S7.

Heatmap of transcription factor network scores from patient HCC1-R from the DOMINO interaction analysis.

Additional file 10.

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 The Creative Commons Public Domain Dedication waiver ( 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

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, S., Yuan, L., Danilova, L. et al. Spatial transcriptomics analysis of neoadjuvant cabozantinib and nivolumab in advanced hepatocellular carcinoma identifies independent mechanisms of resistance and recurrence. Genome Med 15, 72 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Hepatocellular carcinoma
  • Spatial transcriptomics
  • Therapeutic resistance
  • Neoadjuvant therapy
  • Immunotherapy
  • Tumor recurrence