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

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

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

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

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. C–E The immune-rich tumor region highly expresses immune-related pathways that were initially observed to be enriched across responders’ samples. F–H 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

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

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.