Trial of TACE + ICB for Biomarker Discovery
To evaluate the biological determinants of HCC responses following TACE with PD-1 antibody treatment, we evaluated specimens from patients enrolled in a clinical trial of DEB-TACE and sintilimab in locally advanced HCC (15). Patients had HCC of BCLC stage A who exceeded the Milan criteria or had BCLC stage B (Fig. 1A, Table S1-S2, NCT04174781). All patients had Child-Pugh A liver function and did not have vascular invasion, extrahepatic metastasis, or clinically significant portal hypertension. Sixty-one patients enrolled, all of whom underwent at least one treatment cycle combining DEB-TACE and sintilimab (Fig. 1A). The clinical characteristics of trial patients are summarized in Table S1. The majority of the patients were male (52 out of 61, representing 85%), with a median age of 58 years (ranging from 26 to 75 years old). Of these patients, 85% (52 individuals) had two or fewer lesions, while 56% (34 individuals) were diagnosed with BCLC stage B HCC. The median size of the target tumor was 7 cm (ranging from 2.1 to 14.1 cm).
In total, 82% (50 out of 61) of the patients underwent only one treatment cycle with sintilimab combined with DEB-TACE, while 18% (11 patients) received two or three treatment cycles. At the data cutoff date, out of the 60 patients assessed for effectiveness, 37 (62%) exhibited an objective response to DEB-TACE combined with sintilimab; 3 patients had a complete response (5%), 34 a partial response (57%) and 20 experienced stable disease. Overall, the disease control rate was 95% (57/60) with a median follow up of 26 months. Patients showing complete and partial response were categorized as having an objective remission (OR), while patients showing steady (SD) were classified as such (Fig. 1B and Fig. S1A). Of the treated patients, adequate biopsy specimens were collected from 38, enabling analysis of biological correlates (Fig. 1C). Representative computer tomograph of OR (objective response) and SD (stable disease) patients were shown in Fig. 1D.
Immune Suppressive Epithelial Cells Elevate in Responders
We hypothesized that biological markers known to be active in modulating the HCC TME might be associated with treatment response. We generated a comprehensive spatial multi-omics dataset from patient samples. We excluded patients who underwent multiple rounds of TACE or did not have both pre- and post-treatment samples (Fig. 1C). As an initial evaluation, we devised a 40-marker IMC panel (39 protein markers and DNA) (Fig. 2A) based on a list compiled for the study of cellular neighborhoods in HCC (16). To ensure comprehensive coverage of each slide during IMC scanning, we performed scans across multiple regions (Fig. S1B). All antibodies were validated through immunohistochemistry (IHC) before heavy metal conjugation, and each metal-bound antibody was further verified by IMC (Fig. S2-4). The panel included markers for epithelial, endothelial, and stromal cells, various immune cells, as well as the cytokines IL-1β, TNF-α, and IL-6, the proliferation marker Ki-67, and the apoptotic indicator cleaved caspase-3 (Table S2). In total, mass cytometry analysis was done on 459 tumor ROIs from pre and post therapy samples (SD, n=15; OR, n=23).
Marker expression levels in each cell were quantified and converted into a matrix. Following batch correction, approximately 1,895,221 cells were grouped into 19 cell meta-clusters identified by PhenoGraph (Fig. 2A to C). These meta-clusters were then annotated based on marker expression (Fig. 2C, Fig. S5A and B). Macrophage subsets, including resident macrophages (CD11blow CD68+) and infiltrating macrophages (CD11b+ CD68+ CD14+) were separately identified (Fig. 2A). T cell subsets (CD4+ or CD8+), were further refined based on functional markers (Fig. 2A). Additionally, we observed epithelial cells (Pan-Cytokeratin+) and immune-suppressive epithelial cells (PD-L1+/B7H4+ Pan-Cytokeratin+). Fibroblasts were further divided into Collagen Ihi fibroblasts and α-SMAhi myofibroblasts (17) (Fig. 2A). Some clusters contained multiple cell types, such as epithelial cells and T cells, immune-suppressive epithelial cells and resident macrophages, and resident macrophages and fibroblasts, indicating close interactions of neighboring cells. One cluster from lineage-negative cells was also detected (Fig. 2A and Fig. S5). The frequencies of these cells in evaluated patients are shown in Fig. 2D.
We observed a counterintuitive pattern in immune-suppressive epithelial cells, which increased after the treatment only within the OR group (Fig. 2E). In contrast, we identified opposing trends in α-SMAhi myofibroblasts and Collagen Ihi fibroblasts following the combination treatment. The number of α-SMAhi myofibroblasts decreased, while Collagen Ihi fibroblasts increased post-treatment. The increase in Collagen Ihi fibroblasts was less pronounced in the OR group (Fig. 2E). The combination treatment also revealed opposing patterns in Tregs and CD4/8+ T cells, as Tregs experienced a decline while CD4/8+ T cells exhibited a rise. In addition, both resident and CD11c+ macrophages demonstrated a frequency increase after the combined therapy (Fig. 2E). These observations suggest that therapy initiates an unexpected adaptive resistance mechanism. A decrease in α-SMAhi myofibroblasts along with an increase in Collagen Ihi fibroblasts could reflect remodeling of the extracellular matrix that may either support or hinder an antitumor immune response, depending on the context. The observed decrease in Tregs alongside an increase in CD4/8+ T cells aligns with the intended immunostimulatory effects of PD-1 blockade, which aims to bolster antitumor immunity. Additionally, the rise in both resident and CD11c+ macrophages after the combination therapy underscores the dynamic nature of innate immune responses.
Treatment Response is Marked by Organizational Changes in the TME
Beyond tissue composition, we hypothesized that cellular organization may also contribute to or predict for HCC responses to combined therapy. We utilized high-resolution spatial data to identify in situ alterations to cellular neighborhood organization after the combined therapy (18, 19), which we defined as the ten closest neighbor cells surrounding a central cell (Fig. 3A). We analyzed the cellular composition of the neighbors and annotated these neighborhoods based on the primary cellular clusters, resulting in various functional cellular neighborhood units (Fig. 3B). Additionally, we employed a Voronoi diagram, which was applied to each IMC image using the CN composition, to visualize the correlation between the defined CN functional units and the original IMC images (Fig. 3C) (20). The IMC topology map, depicted by the Voronoi diagram with different CNs represented by distinct colors (as shown in Fig. 3C), aligned well with the IMC image, confirming that cellular neighborhood assignments were correct.
Furthermore, we adapted previously described methods to identify important regions within the HCC TME (21, 22). An expansive homotypic cellular cluster, as identified by uniform cellular phenotypes, was demarcated as a discrete spatial domain within the TME. Notable regions included the epithelial region, immune-suppressive epithelial region, fibroblast region, and the contiguous interfaces where immune-suppressive epithelial cells/immune cells, as well as regions of epithelial cells and fibroblast co-localization (Fig. 3C). Similar to the bulk frequency analysis, we found that the immune-suppressive epithelial region increased in the OR group after treatment, while the epithelial region decreased in the OR group. The fibroblast region showed an increase in both the SD and OR groups, although the magnitude of increase was less pronounced in the OR group. Additionally, the immune region experienced an increase in the OR group post-treatment. The epithelial/fibroblast intersection region and the immune suppressive epithelial/immune region displayed contrasting patterns in the OR group after treatment, with the former decreasing and the latter increasing (Fig. 3C and D).
We performed interaction analysis to uncover potential mechanisms of resistance in the SD group. In that group, we observed a higher number of cellular interactions involving fibroblasts and epithelial/fibroblast interactions. Conversely, cellular interactions involving infiltrating macrophages, NK cells, and various T cells were notably absent in the epi-fibroblast intersection regions (Fig. 3E). On the other hand, in the immune suppressive-immune intersection regions, we observed higher cellular interactions involving immune suppressive epithelial cells and various T cells, while fibroblast interactions were comparatively reduced (Fig. 3F). This analysis suggests that distinct cellular subsets orchestrate environmental factors that give rise to contrasting regional features. Specifically, immune suppressive epithelial cells are in position to mediate cellular interactions with immune cells, particularly various T cells, in the OR group. In contrast, the SD group lacked these cellular interactions. These patterns suggest that immune-suppressive epithelial cells recruit T cells, whereas epithelial cells promote fibroblast enrichment, which may reduce the efficacy of therapeutic responses.
Spatial and Single-Cell Transcriptomics Identify cGAS as a Regulator of the TME
Our IMC analysis allowed us to visualize and quantify how cellular neighborhoods respond to TACE + ICB. To discover the biological processes that regulate these changes, we used single-cell RNA sequencing and spatial transcriptomics to evaluate samples (Fig. S6A). For this post-treatment samples from 6 patients in the OR and 4 patients in the SD cohorts were evaluable. We captured a total of 92,785 single cells from 6 OR and 4 SD patients post-surgery, with high sequencing quality. We identified 10 primary cell populations (Fig. 4A) based on specific markers (Fig. 4B), such as EpCAM+ epithelial cells, VWF+ endothelial cells (23) , and ACTA2+ fibroblasts (24) (Fig. 4B). We also identified various immune cell subsets, including CD3+ T cells, CD79B+ B cells, IGHG1+ plasma B cells, CD68+ macrophages, GNLY+ natural killer (NK) cells (25), and CPA3+ mast cells (26) (Fig. 4B). Additional clustering of T cells revealed 5 subsets (Fig. 4C and D): Granzyme A/K+ CD8+ effector T cells, CD127+ CCR7+ central memory CD8+ T cells (CD8+ Tcm), Granzyme A/K+ IL7R+ effector memory CD8+ T cells (CD8+ Tem) (27) , CTLA4+ exhausted CD4+ T cells, and FOXP3+ regulatory T cells (Treg) (Fig. 4C and D). We identified two major subsets of macrophages: resident macrophages and infiltrating macrophages (Fig. 4E and F). Resident macrophages exhibited distinct markers, such as C1QA/B, which distinguished them from infiltrating macrophages (28). In contrast, infiltrating macrophages displayed higher levels of CCR2 and S100A8/9, than resident macrophages (Fig. 4E and F)
Similar to the IMC analysis, we discerned two distinct subsets of epithelial cells within the TME: immune suppressive epithelial cells expressing PD-L1, and a separate population of epithelial cells devoid of PD-L1 expression (Fig. 4G and H). The PD-L1 population plays a regulatory role in immune recognition, whereas the PD-L1 null population, represents a phenotypically different group of epithelial cells within the same tissue context. We classified fibroblasts into two subtypes based on the IMC results: Collagen Ihi fibroblasts and α-SMAhi myofibroblast cells (Fig. 4I and J). Our scRNA-seq landscape encompassed the key cellular components of the HCC TME and are consistent with the IMC results. These findings provide a detailed accounting of the diverse cellular populations in the TME and frame their potential roles in tumor immune responses.
To unravel the mechanisms underlying the spatial and cellular interactions that dictate the regional response to TACE and ICB treatment, we performed an integrative analysis of single-cell RNA sequencing (scRNA-seq) and corresponding spatial transcriptomics (29). We employed a robust cell type decomposition (RCTD) approach, a computational method that harnesses the cell type profiles acquired from scRNA-seq data to decompose cell type mixtures while accounting for discrepancies across sequencing technologies. We used RCTD to decompose the spots within the spatial transcriptomics of all tissue sections, based on the refined immune and non-immune cellular clusters. After single-cell spatial deconvolution, we evaluated each spot for critical cell types identified by the IMC analysis, including epithelial cells, immune suppressive epithelial cells, fibroblasts, and T cells (Fig. S6B-E). This approach let us to accurately assess the presence and distribution of these key cell types across the tissue sections, providing insights into the spatial organization of the TME. As Fig. S6B-E shows, epithelial cells and fibroblast spots were significantly enriched in the SD patients. By contrast, immune suppressive epithelial cells and T cell spots were heavily enriched in the OR patients. Notably, we observed double positive spots (green spots) indicating the co-occurrence of immune suppressive epithelial cells and T cells was prominently observed in the OR patients (Fig. 4K and Fig. S6F). Conversely, double positive spots for epithelial cells and fibroblasts were predominantly enriched in the SD patients (Fig. 4L and Fig. S6G). These spatial transcriptomic insights suggest that immune-suppressive epithelial cells recruit T cells, enhancing therapeutic efficacy, while epithelial cells promote fibrosis, impeding treatment response.
To understand the biological determinants this organization, we conducted KEGG pathway analyses on the spots with dual positivity for immune suppressive epithelial cells and T cells (Fig S7A). In the OR group, immune activation pathways such as chemokine signaling, T cell differentiation/trans-endothelial migration, and NF kappa B were evident. We also observed immune checkpoint signaling following immune activation (Fig. S7A). We also noted activation of cytosolic DNA sensing (cGAS pathway) pathways in the top group (Fig. S7A). In contrast, in spots with dual positivity for epithelial cells and fibroblasts (Fig. S7B), we observed metabolism-related pathways and oxidative phosphorylation pathways. Particularly noteworthy was the identification of the PPAR signaling pathway, which plays a key role in fibrosis (30), in the SD group (Fig. S7B).
To validate integrated single-cell and spatial transcriptomic findings, we performed bulk RNA transcriptomic analysis on post-treatment samples obtained from a separate cohort of 19 patients treated with the same regimen (Table S1). Among the top differential genes identified, one was cGAS, a major regulator of double stranded DNA recognition (31) (Fig. S7C). To demonstrate the role of cGAS-STING in regulating immune suppression in our patient samples, we conducted several analyses. First, we performed correlation analysis for all the pathways identified in Fig. S7A, achieved by employing gene signatures specific to each pathway for Gene Set Variation Analysis (GSVA). The gene signatures served as input for calculating the GSVA scores. Analysis focused on immune suppressive epithelial and T cell spots across all pathways outlined in Fig. S7A. Subsequently, these GSVA scores were utilized to compute the correlations among the pathways. We noticed the cGAS-STING pathway correlated positively with the most enriched molecular pathways, including those for T-cell recruitment and migration (Fig. S7D), with the exception of the pyroptosis pathway (Fig. S7D). Second, we used SCENIC (single-cell regulatory network inference and clustering) analysis to identify key transcription factors that form the gene regulatory network within immune suppressive epithelial cells (32). We detected cGAS-STING related transcription factors such as IRF3 and NF-κB (Fig. S7E) within immune suppressive epithelial cells. We validated the SCENIC results through Western blotting of individual patient tumor samples which showed elevated levels of phosphorylated versions of STING pathway proteins, TBK1, IRF3, STING, and RelA in the OR group but not the SD group (Fig. 4M). These findings suggest the STING pathway is activated in responsive samples, reinforcing the idea that STING participates in the treatment response. Consistently, mIHC staining also showed more cGAS signal and downstream RelA (33) signal in the OR group, while these signals were absent in the SD group (Fig. 4N). Epithelial cells in the OR group displayed high levels of PD-L1 expression (Pan-CK+) compared to the absence of PD-L1 expression in epithelial cells in the SD group (Fig. 4O). Using mIHC we observed colocalization of PD-L1 and the cGAS expression (Fig. 4N). We note that this connection between cGAS-STING pathway and immune checkpoint expression has been observed previously (34).
To further confirm the association between cGAS signaling and PD-L1 expression, we exposed a panel of three human HCC cell lines to the topoisomerase inhibitor epirubicin, which causes double-stranded DNA (dsDNA) breaks. This caused cGAS signaling and upregulation in PD-L1 expression (Fig. S8), consistent with a previous report (35). This may partially explain why single administration of STING agonists have not yielded satisfactory results in clinical trials since upregulation of PD-L1 could theoretically attenuate the anti-tumor immune response (36). Nevertheless, PD-L1 expression alone does not faithfully predict for immunotherapeutic response in hepatocellular carcinoma (37).
We further analyzed the samples to find actionable hypotheses related to therapeutic vulnerabilities. Utilizing spatial transcriptomics to analyze samples from OR and SD cohorts, we delineated distinct spatial features unique to each group. In the OR cohort, regions of overlap between immune suppressive epithelial cells and T cells were pronounced, indicating a potential niche for immune modulation. Conversely, the SD cohort was characterized by a notable co-localization of epithelial cells and fibroblasts, suggesting a tissue composition of fibrosis. To further elucidate the cellular interactions underlying the recruitment of T cells by immune-suppressive epithelial cells and the promotion of fibrosis by epithelial cells, we performed a ligand-receptor analysis on both OR and SD patients, integrating data from both Spatial Transcriptomics (ST) and single-cell RNA sequencing (scRNA-seq). We used NicheNET, a tool that combines expression data from intercommunicating cells with pre-existing knowledge of signaling and gene regulatory networks, to predict potential ligand-target associations (38). First, we selected our genes of interest from the differentially expressed genes, which we calculated based on ST datasets for specific areas identified via IMC analysis. Next, we applied scRNA-seq for cellular interaction analysis via NicheNET. Our observations indicated that PD-L1+ epithelial cells could secrete multiple cytokines and chemokines, such as CCL5 (39), MIF (40), and LTB (41), thereby recruiting immune cells such as T cells (Fig. S9A). In contrast, in SD patients, epithelial cells transmitted numerous inhibitory signals for PPAR signaling pathways to fibroblasts (Fig. S9B), such as ATF3 (42), BCL6 (43), LCOR (44), EZH2 (45), NR0B2 (46), SIRT1 (47), GPS2 (48), and NRIP1 (49). Accordingly, based on ST analysis, we detected reduced PPAR signaling within the epithelial-fibroblast regions associated with SD patients (Fig. S9C), especially in areas showing dual-positive signals from both epithelial cells and fibroblasts. Our examination revealed downregulation of genes, such as PLIN2 (50), FABP6 (51), PLTP (52), UBC (53), and HMGCS1 (54) (Fig. S9C). These findings indicate a suppression of PPAR signaling in the non-responder (SD) patients.
To validate NicheNET findings, we evaluated levels of double-stranded DNA (dsDNA) and subsequent activation of the cGAS signaling pathway using mIHC and found these were elevated in the OR group, along with expression of PD-L1. This activation was associated with increased recruitment and activation of T cells (CD45RO) (Fig. 4O and Fig. S9D). T cells also showed resident marker CD103 (Fig. S10). Notably fibrosis, as indicated by Collagen I staining, was relatively low in the OR group. Conversely, the SD group exhibited high levels of fibrosis, accompanied by low dsDNA and cGAS signaling. As a result, activated T cells were scarce in the SD group (Fig. 4O and Fig. S9D). These observations were consistent with the pathway analysis obtained through the integration of scRNA-seq and ST, which showed an interaction between dsDNA, cGAS signaling, fibrosis, and T cell activation.
In conclusion, the above findings provide a model for how cGAS signaling can modulate responses to TACE + ICB in HCC. Activation of cGAS signaling by double-stranded DNA (dsDNA) results in increased recruitment and activation of T cells, as well as upregulation of PD-L1 expression in the responsive (OR) group. Conversely, non-responders (SD) group exhibit fibrosis, low dsDNA, and limited cGAS signaling, leading to reduced T cell recruitment and PD-L1 expression. These results highlight the significance of cGAS signaling and fibrosis in shaping the TME and suggest its potential as a therapeutic target in HCC (Fig. 4P).
Addition of Anti-Fibrotic Agents Enhance Therapy for HCC
The observation that fibrosis limited the efficacy of TACE+anti-PD-1 treatment suggests that anti-fibrotic agents might enhance therapy. We evaluated this concept in a mouse HCC model induced by cMyc and Nras expression (16) since in murine models of HCC, activation of both Myc and Ras oncogenes has been implicated in the promotion of hepatic fibrosis (55, 56). Treatments were conditioned starting on day 14 with Elafibranor, a PPAR agonist known to have anti-inflammatory and anti-fibrotic effects. Elafibranor has shown promising results in reducing liver inflammation, oxidative stress, and fibrosis in preclinical liver disease models (57). On day 21, we mimicked TACE by treating with the STING agonist cGAMP; we also administered anti-PD-1 therapy (Sintilimab). After 30 days, we harvested the liver with tumors for further analysis (Fig. 5A).
Both Elafibranor and the dual combination therapy of TACE+anti-PD-1 demonstrated reductions in tumor burden, with the dual combination therapy showing a greater reduction. However, the triple combination therapy involving Elafibranor, TACE, and anti-PD-1 exhibited the least tumor burden (Fig. 5B and C). Treatment with Elafibranor alone resulted in reduced fibrosis. The dual combination therapy of TACE+anti-PD-1 showed an increase in CD8+ T cells. Notably, the triple combination therapy led to reduced fibrosis and increased CD8+ T cells (Fig. 5D and Fig. S11).
Graphical Classification Based on Spatial Similarity (GCSS) Predicts Treatment Response
An enticing clinical strategy would be to predict which patients are likely to have a poor response to therapy due to fibrosis, then add an anti-fibrotic agent for that population. We hypothesized that it might be possible to generate a predictive model using imaging mass cytometry (IMC) data obtained from needle biopsies. We used pre-treatment IMC images and associated clinical outcomes from the TACE+ICB trial as the training dataset for our model (Fig. 6A).
The model was constructed by dividing segmented cells within each Region of Interest (ROI) into separate images (N=29), each reflecting a mask of specific cellular identity. Each image represents the spatial localization of a specific cell type. Subsequently, we calculate the pairwise spatial similarity between different cell types for each ROI, using a structural similarity matrix (SSIM) measurement (Fig. 6B). We then employed a network classification method (58), Graphical Classification based on Spatial Similarity (GCSS), to predict for association between spatial patterns of cells, including interactions, and the efficacy of therapy (Fig. 6A). The GCSS model considers both the spatial similarity information of different cell types and the network structure of the data in a computationally efficient way. To confirm the GCSS method was optimal, we conducted a comparison between the GCSS method and the classical support vector machine (SVM) and Least Absolute Shrinkage and Selection Operator (Lasso) methods. The results showed that GCSS exhibited higher sensitivity and accuracy compared to SVM and Lasso (Fig. 6C-D).
Incidentally, GCSS also provided confirmatory mechanistic data regarding cellular interactions within each group. In the SD group, we observed interactions between epithelial/fibroblast and CD45RA+ CD8+ T cells, effector CD8+ T cells, and epithelial/fibroblast cells, as well as neutrophil and epithelial cell interactions. In the OR group, we found interactions between NK cells and infiltrating macrophages, effector CD8+ T cells and immune suppressive epithelial cells, and GATA3+ CD4+ T cells and resident macrophages. These interactions are consistent with our integrated spatial multi-omics analysis and prior mechanistic exploration (Fig. 6E).
We used dominant patterns in the model to create a minimized mIHC panel for predicting treatment response consisting of six markers, including CD4, CD8, CD68, Collagen I, Pan-CK, and PD-L1 (Fig. 6F). These markers represented the crucial cell interactions responsible for different treatment responses within the OR and SD groups (Fig. 6E). This approach is intended to enable clinical feasibility since it reduces assay costs and can be done more rapidly.
To validate the predictive capability of the minimal panel, we collected 61 pre-treatment specimens from HCC patients who underwent dual combination treatment. Using mIHC input into the GCSS model, we achieved over 80% accuracy in predicting treatment response, similar to results on IMC data (compare Fig. 6C to 6G). These results suggest the mIHC panel will be effective for predicting treatment responses for HCC to TACE and ICB at minimal cost. This result was achieved despite a relatively small sample size, suggesting robustness in the method.