Abstract

Background & Aims

   Despite atezolizumab plus bevacizumab being a standard treatment for
   advanced hepatocellular carcinoma (HCC), a significant proportion of
   patients do not achieve durable benefit. This study aimed to identify
   predictive biomarkers for this therapy by investigating the role of
   immune activation within the tumor microenvironment (TME).

Methods

   We characterized the intratumoral TME of patients with advanced HCC
   treated with atezolizumab plus bevacizumab using single cell
   transcriptomics on pretreatment tumor biopsies from 12 patients. To
   complement and support these findings, we integrated our single cell
   data with publicly available bulk RNA-sequencing data from independent
   clinical trial cohorts.

Results

   Patients who responded to combination therapy with atezolizumab plus
   bevacizumab demonstrated an immune-activated TME, marked by enhanced
   cytotoxicity and a tumor-specific T cell response. These patients also
   exhibited an increased proportion of inflammatory cytokine-enriched
   tumor-associated macrophage clusters with stronger interactions with T
   cells, an increased population of conventional dendritic cells, and
   activated antigen-presenting function in tumor endothelial cells. When
   publicly available bulk RNA-sequencing data from independent clinical
   trial cohorts were analyzed, these immune activation features were
   associated with improved progression-free survival (median 10.8 months,
   95% CI: 7.3–not reached versus 5.5 months, 95% CI: 4.0–6.7; p <0.001).

Conclusions

   These findings suggest that the existence of an activated immune TME
   before treatment is crucial for a favorable clinical response in
   patients with HCC treated with atezolizumab plus bevacizumab.

Impact and implications:

   Only a subset of patients with HCC benefit from combination therapy
   with atezolizumab plus bevacizumab, limiting its clinical utility. In
   this study, we used single cell RNA analysis to identify TME features
   associated with a clinical response to this therapy. Our findings
   suggest that a pre-existing immune-activated TME is crucial for
   predicting the response to atezolizumab plus bevacizumab. Specifically,
   features such as enhanced T cell cytotoxicity, inflammatory
   cytokine-enriched macrophage clusters, active antigen presentation in
   endothelial cells, and an increased presence of dendritic cells may aid
   patient selection and inform therapeutic strategies.

   Keywords: RNA, Single cell, Immunotherapy, Immune activation,
   Biomarkers, Macrophages, Tumor associated

Graphical abstract

   [49]Image 1
   [50]Open in a new tab

Highlights:

     * •
       An activated immune microenvironment in advanced HCC is associated
       with treatment response to atezolizumab plus bevacizumab
       combination therapy.
     * •
       Responders exhibit increased cytotoxic T cells, inflammatory
       cytokine-rich TAMs, cDC2 cells, and antigen presentation by TECs.
     * •
       The 10-gene sc_ABRS signature, identified through single-cell
       analysis, is associated with PFS in the GO30140/IMBRAVE150
       database.

Introduction

   The combination of atezolizumab and bevacizumab is now a standard
   first-line therapy for advanced hepatocellular carcinoma (HCC).[51]^1
   This immunotherapy and anti-angiogenic combination has shown
   significant survival benefits compared with previous standards of
   care,[52]^2^,[53]^3 because atezolizumab enhances antitumor immunity
   and bevacizumab targets angiogenesis, a key driver of tumor growth and
   immune evasion in HCC.[54]^4^,[55]^5 However, the objective response
   rate to this combination therapy is only ∼30%, highlighting the crucial
   need to identify biomarkers that can predict response and optimize
   treatment strategies in advanced HCC.

   Although several biomarkers, such as programmed death-ligand 1 (PD-L1)
   expression, tumor mutational burden, and microsatellite instability,
   have been proposed as predictors of immunotherapy response in various
   cancers,[56][6], [57][7], [58][8] their predictive value in HCC has yet
   to be fully established.[59]^9 This uncertainty may stem from the
   inherent heterogeneity of HCC, characterized by diverse systemic
   immunity, pathogenesis, histopathology, and molecular
   features.[60]^10^,[61]^11 Moreover, the complex interplay within the
   tumor microenvironment (TME) can modulate treatment response and
   potentially contribute to drug resistance.[62][12], [63][13], [64][14]
   Therefore, a thorough understanding of these intricate interactions is
   essential for accurately predicting immunotherapy response in HCC.

   Single cell-level analysis is a highly effective method for evaluating
   the heterogeneity of cellular components in the TME.[65]^15 However, a
   comprehensive investigation of the correlation between treatment
   response to atezolizumab plus bevacizumab and TME at the single cell
   level in advanced HCC has been hindered by the challenging nature of
   procuring fresh tumor tissue. In this study, we conducted single cell
   RNA sequencing (scRNA-seq) of freshly acquired biopsy samples from
   patients with advanced HCC to analyze cellular and molecular mechanisms
   and identify predictive markers for clinical response to atezolizumab
   plus bevacizumab.

Materials and methods

Sample collection

   This study was approved by the Institutional Review Board of the
   Samsung Medical Center. Written informed consent was obtained from each
   patient (IRB No. 2021-01-021-009). Between June 2020 and February 2023,
   tumor biopsy samples were collected from 18 patients with advanced HCC,
   including patients with Barcelona Clinic Liver Cancer (BCLC) stage B
   HCC with diffuse, bilobar involvement for whom locoregional treatments
   were not feasible, and patients with BCLC stage C HCC. Pretreatment
   core-needle biopsies were performed on 17 patients. For one patient who
   had undergone surgical resection but subsequently developed pulmonary
   metastases, tumor-infiltrating lymphocytes (TILs) were collected for
   analysis.

ScRNA and T cell receptor V(D)J sequencing

   The Chromium Single Cell V(D)J Solution from 10x Genomics was used for
   single-cell 5′ gene expression and V(D)J sequencing. Cell suspensions
   were loaded onto the Chromium Single Cell Controller to create single
   cell gel beads in emulsion (GEMs) using Chromium Single Cell V(D)J
   Reagent Kits. The resulting libraries for single-cell 5′ gene
   expression and V(D)J were sequenced on the Illumina Novaseq 6000
   platform with 150 paired-end reads. Alignment of the scRNA and V(D)J
   sequencing data to the hg38 genome and the generation of gene-barcode
   metrics were performed using Cell Ranger software (v.5.0.0).

ScRNA-seq analysis

   To ensure high-quality cells, those expressing fewer than 300 genes and
   those with >20% unique molecular identifier (UMI) tags from
   mitochondria gene-expressed cells were filtered out. After this
   quality-filtering step, the remaining counts underwent normalization,
   dimension reduction using principial components analysis (PCA),
   integration, shared nearest neighbor (SNN) clustering, and uniform
   manifold approximation and projection (UMAP) visualization of
   scRNA-sequencing data using the Seurat R package (v4.3.0.1) with
   default options. Cluster annotations were based on marker gene lists,
   including PTPRC, CD3D, LYZ, VWF, COL1A2, and other relevant genes. For
   the analysis of differentially expressed genes (DEGs), the
   FindAllMarkers or FindMarkers function from Seurat was used with the
   MAST test method.

TCR V(D)J analysis and clonotype identification

   Using filtered annotation T cell receptor (TCR) data from Cell Ranger
   (v5.0.0), only productive CDR3 sequences matched with UMIs from
   scRNA-seq data were selected for subsequent analysis. New clonotype IDs
   were generated by concatenating the CDR3 protein sequences of alpha and
   beta chains, taking into account all possible chain combinations.
   Instances of multiple chains, such as two beta chains and one alpha
   chain, or two alpha chains and one beta chain, were all considered in
   the subsequent analysis. To quantify the abundance of clonality, the
   number of cells with the same clonotype ID was counted. The Gini index,
   used as a measure of clonality, was calculated using the Vegan package
   in R.

Gene set enrichment pathway analysis

   For functional pathway analysis, we performed gene set enrichment
   pathway analysis (GSEA) using the DEG list. This analysis was conducted
   using the clusterProfiler and fgsea package in R, with the curated gene
   sets (C2) or ontology gene sets (C5) from the Molecular Signature
   Database v7.2 ([66]www.gsea-msigdb.org/gsea/msigdb). For specific gene
   sets, such as tissue-resident T cell markers[67]^16 or tumor-specific
   markers,[68]^17 we used gene sets listed in related publications.

Cell–cell interaction analysis

   To investigate cell–cell interactions between cell types, we used
   CellphoneDB[69]^18 (v4) to identify significant ligand–receptor
   combinations using CellphoneDB Data (v4.1.0). The interaction score was
   calculated based on the mean expression values of individual
   ligand–receptor pairs across corresponding cell types within clinical
   groups, where at least 10% of cells expressed both the ligand and
   receptor genes. Statistical significance was assessed through a
   permutation test with 1,000 iterations, with other parameters set to
   their default options.

Gene expression signature analysis

   The phagocyte, M1, M2, and angiogenesis signature scores in the
   scRNA-seq analysis were calculated using the gene sets described in
   [70]Table S1 and the AddModuleScore function in the Seurat package. To
   examine the correlation between specific cell clusters and patient
   survival in bulk RNA-seq data, a mean of marker gene expression was
   used, and the groups were categorized based on the 50th percentile.

Multiplex immunofluorescence

   Formalin-fixed paraffin-embedded (FFPE) tissue sections (n = 9) from
   biopsy samples collected for single cell analysis were subjected to
   multiplex immunofluorescence (mIF) staining. Sections were cut at 4 μm,
   deparaffinized, rehydrated, and stained following the manufacturer’s
   protocol. mIF staining targeted DAPI, CD68, CD163, and pan-cytokeratin
   (CK) using the Opal 7-color fluorescence IHC kit (PerkinElmer, Waltham,
   MA, USA) and the BOND RX system (Leica Microsystems, Vista, CA,
   USA).[71]^19 Antibodies and fluorophores included CK (AE1/AE3, Dako;
   Opal 620), CD68 (KP1, Dako; Opal 690), CD163 ([72]EPR19518, Abcam; Opal
   570), and CCL3 (1H20L19, Invitrogen; Opal 520). Spectral images of
   multiplex-stained slides were captured with the Vectra Polaris
   Automated Quantitative Pathology Imaging System (PerkinElmer).

   Of the nine slides evaluated, one was excluded because of suboptimal
   staining, leaving eight slides for analysis. For these, 15 regions of
   interest (ROIs) were selected (1–3 ROIs per slide). Cell segmentation
   and annotations for CD68, CD163, CCL3, and CK were performed in the
   ROIs, with 8–19 cells annotated per antibody per ROI to train the
   software in cell classification. The inForm software (version 3.0.0,
   Akoya Biosciences, Marlborough, MA, USA) determined the positivity or
   negativity of each antibody. Cells positive for CD163 and negative for
   CK were classified as tumor-associated macrophages (TAMs), and the
   proportion of CCL3-positive TAMs was calculated. All analyses were
   performed by experienced pathologists.

Integration of the publicly available bulk RNA-seq dataset from the GO30140
phase Ib and IMbrave150 phase III clinical trials

   To further explore the clinical relevance of our single cell-derived
   findings on progression-free survival (PFS) and overall survival (OS),
   we analyzed a publicly available bulk RNA-seq dataset. This included
   311 tumor samples[73]^20 from patients with advanced HCC treated with
   either atezolizumab plus bevacizumab (n = 253) or sorafenib (n = 58) in
   the phase Ib (GO30140)[74]^5 and phase III (IMbrave150)[75]^2 clinical
   trials. Raw fastq data (EGAD00001008128) and corresponding clinical
   data (EGAD00001008130) were obtained from the European Genome-Phenome
   Archive. The raw reads were aligned to the human reference genome
   (GRCh37) using STAR v2.5.1b with default parameters, and gene
   expression quantification was performed using RSEM software. Samples
   with low library sizes (<1 million reads) were excluded from the
   analysis, resulting in 217 patients treated with atezolizumab plus
   bevacizumab being included in the final dataset. Clinical
   characteristics are summarized in [76]Table S2.

Survival and statistical analysis

   Survival analyses were conducted using the Cox proportional hazards
   model with a log-likelihood ratio test to evaluate overall model
   significance. Kaplan-Meier survival analysis stratified samples into
   high- and low-signature groups. Statistical analyses were performed in
   R (version 4.1.2). The Wilcoxon test compared continuous variables
   between two groups. Statistical significance was indicated as follows:
   ∗p <0.05, ∗∗p <0.01, and ∗∗∗p <0.001.

Results

Profiling of 10 major cell types in advanced hepatocellular carcinoma

   To delineate the single cell landscape of advanced HCC, we analyzed 18
   tissue samples, including 17 tumor biopsies and one TIL sample, with
   scRNA-seq using the 10x Genomics platform ([77]Fig. 1A). Clinical
   characteristics are detailed in [78]Table S3; 12 patients were treated
   with first-line atezolizumab plus bevacizumab. All but two of these
   patients were treatment naive. In total, 63,936 high-quality cells
   (average 3,552 per sample) were analyzed, revealing 10 distinct
   cell-type clusters based on prior HCC studies[79]^21^,[80]^22
   ([81]Fig. 1B,C): C1 myeloid cells (PTPRC and CD68); C2 plasmacytoid
   dendritic cells (DCs; CLEC4C and CD74); C3 B cells (PTPRC, CD79A, and
   MS4A1); C4 plasma cells (CD79A and MZB1); C5 T/natural killer (NK)
   cells (PTPRC, CD3D, CD4, CD8A, and KLRF1); C6 fibroblasts (ACTA2 and
   COL1A1); C7 endothelial cells (FCN3, PECAM1, and VWF); C8 malignant
   hepatocytes (HP and KRT7); C9 mast cells (KIT); and C10 cycling cells
   (STMN1 and MKI67).

Fig. 1.

   [82]Fig. 1
   [83]Open in a new tab

   Single-cell transcriptomic profiling and tumor response to atezolizumab
   plus bevacizumab in HCC.

   (A) Flowchart of the study. (B) UMAP visualization displaying 10
   distinct cell clusters, each represented by a unique color. (C) Dot
   plot of canonical markers for each cluster. Circle size indicates the
   percentage of cells expressing each marker, and color represents the
   scaled average expression level. (D) Waterfall plot of the maximum
   change in tumor size. Statistical analysis was performed using the
   Wilcoxon test. (E) Proportion of major clusters by immunotherapy
   response group, excluding malignant clusters (hepatocytes). The x-axis
   represents the response group, with statistical analysis performed
   using the Wilcoxon test. UMAP, uniform manifold approximation and
   projection.

   To identify biomarkers, we performed transcriptomic analysis on 12
   patients treated with atezolizumab plus bevacizumab. Patients were
   classified as responders (partial response [PR], n = 4) or
   nonresponders (stable disease [SD]/progressive disease [PD], n = 8)
   based on modified Response Evaluation Criteria in Solid Tumors
   (mRECIST) criteria[84]^23 ([85]Fig. 1D). The overall response rate was
   33.3% (complete response [CR]: 0; PR: 4), with 25.0% SD and 41.7% PD.
   The proportions of cell types varied substantially across samples,
   possibly reflecting underlying differences in tumor phenotypes or
   biopsy locations ([86]Figs. S1A and B). No significant differences were
   observed in major cell type proportions between responders and
   nonresponders ([87]Fig. 1E), although immune cell density tended to be
   higher in responders ([88]Fig. S1C).

Enriched cytotoxic and tissue-resident T cell signatures in responders to
atezolizumab plus bevacizumab

   To further investigate the mechanisms underlying the response to
   atezolizumab plus bevacizumab therapy in HCC, we focused on dissecting
   T cell clusters. A total of 27,717 T or NK cells were identified and
   categorized into 14 clusters, including five CD4 T cell subsets, four
   CD8 T cell subsets, two gamma delta T cell subsets, one
   mucosal-associated invariant T, and two NK cell subsets ([89]Fig. 2A).
   We also performed single cell TCR sequencing (scTCR-seq) to assess T
   cell diversity and clonality, yielding 41,884 TCR sequences with a 90%
   match to T cells ([90]Fig. S2A). The analysis of number of T cells
   (abundance) sharing unique CDR3 (Clone ID) indicated an increase in
   clonality score (Gini index) within CD8 T cells (C6–C9 clusters) in the
   responder group, although this did not reach statistical significance
   (p = 0.11) because of the small sample size ([91]Fig. 2B,C).

Fig. 2.

   [92]Fig. 2
   [93]Open in a new tab

   Enhanced cytotoxicity and tissue-resident T cell presence in
   immunotherapy responders.

   (A) UMAP visualization of T cells from 18 patients with HCC, showing
   14 T and NK cell clusters, each represented by a distinct color. (B)
   UMAP visualization of clonal abundance, indicating the number of cells
   with the same TCR clone ID. (C) Boxplots showing the degree of
   clonality (Gini index) across different response groups. Statistical
   comparisons were performed using the wilcoxon tests. (D) Heatmap of
   relative expression levels for canonical T and NK cell markers. (E)
   Volcano plot displaying differential gene expressions in CD8 T cells
   between responder and nonresponders, specifically for clusters 6–9. (F)
   Enrichment plot of GO terms related to cell killing, leukocyte-mediated
   cytotoxicity, tissue-resident memory gene signatures, and
   tumor-specific T cell signatures (Melanoma antigen vs. EBV or CMV
   virus) in the responder group compared with the nonresponder group in
   CD8 T cells. CMV, cytomegalovirus; EBV, Epstein–Barr virus; HCC,
   hepatocellular carcinoma; NK, natural killer; TCR, T cell receptor;
   TRM, tissue-resident memory T cell; UMAP, uniform manifold
   approximation and projection.

   Within the CD8+ T cell clusters, distinct subpopulations were
   delineated, including central memory T cells (Tcm) expressing IL7R (C6,
   CD8_Tcm_IL7R^+), effector memory T cells (Tem) expressing CCL4 (C7,
   CD8_Tem_CCL4^+), Tem cells expressing GZMK and exhaustion markers, such
   as PDCD1, ICOS, and HAVCR2 (C8, CD8_Tem_GZMK^+), and cytotoxic T cells
   (Tcyt) expressing GZMH (C9, CD8_Tcyt_GZMH^+) ([94]Fig. 2D). When
   comparing T cell subclusters between responders and nonresponders, the
   CD4_Tcm_IL7R+ (C2) and CD8_Tcm_IL7R+ (C6) clusters were significantly
   lower in responders than in nonresponders (p <0.05) ([95]Fig. S2B).

   To further explore the differences between responders and
   nonresponders, we performed DEG analysis within CD8+ T cell clusters.
   Responders exhibited higher expression of exhaustion molecules (such as
   PDCD1, ICOS, and HAVCR2) and cytotoxic molecules (such as FGFBP2, GNLY,
   and GZMB) compared with nonresponders ([96]Fig. 2E). GSEA revealed
   significant enrichment of tumor cell killing-related pathways, such as
   leukocyte-mediated cytotoxicity, in responders ([97]Fig. 2F). Moreover,
   gene sets associated with tissue-resident T cells were significantly
   enriched in the responder group, suggesting a higher presence of such
   cells in responder patients (Normalized Enrichment Score [NES] = 1.74,
   p = 0.034). In addition, considering a predominant infection by HBV in
   our cohort, we validated that these differential gene sets were tumor
   specific rather than related to viral infection (NES = 1.53, p =
   0.004). These findings collectively emphasize T cell activation of
   cytotoxicity-related pathways, abundance of tissue-resident T cells,
   and enrichment of tumor-specific signatures in the responder group
   ([98]Fig. 2F). When we computed signature scores for cytotoxicity and
   exhaustion, we found that responders had significantly higher
   cytotoxicity and exhaustion scores in CD8+ T cells (both p <0.001)
   ([99]Fig. S2C).

Proinflammatory TAM_activating_CCL^+ cluster enhances T cell recruitment and
predicts response to atezolizumab plus bevacizumab

   We next performed myeloid cell clustering and generated merged UMAP
   plots ([100]Fig. 3A), identifying two distinct C1Q+TAM clusters. The
   resting TAM cluster (C1, TAM_resting_C1Q^+) expressed C1Q genes (C1QA,
   C1QB, and C1QC) and GPNMB, whereas the activating TAM cluster (C2,
   TAM_activating_CCL^+) showed proinflammatory chemokine expression
   (CCL3, CCL4, and CXCL2) along with C1Q genes ([101]Fig. S3A).
   Additional clusters included MDSC-like_EREG^+ cells (C3) expressing
   EREG, FCN1, and S100A4, and Neutrophil_CSF3R+ cells (C4) with high
   neutrophil marker expression (CSF3R, G0S2, and NAMPT) ([102]Fig. 3B).
   Both TAM_resting_C1Q^+ and TAM_activating_CCL^+ clusters expressed
   phagocytosis-associated genes (MRC1, CD163, MERTK, and C1QB), although
   TAM_activating_CCL^+ showed distinctly higher M1-associated features
   ([103]Fig. 3C).

Fig. 3.

   [104]Fig. 3
   [105]Open in a new tab

   Enhanced inflammatory features in myeloid cells of responders.

   (A) UMAP of the distribution of myeloid cells, including macrophages,
   dendritic cells, and other myeloid clusters. Each cluster is indicated
   by a unique color. (B) Heatmap illustrating the expression levels of
   marker genes across different myeloid cell clusters. Colors represent
   the relative expression levels within each cell. (C) Gene-set
   signatures related to phagocytosis, M1 (proinflammatory), M2
   (anti-inflammatory), and angiogenesis scores within the four
   macrophage/monocyte-related clusters (clusters 1–4). (D) Boxplot of the
   proportion of each cell cluster within macrophage-related clusters
   (clusters 1–4). Statistical significance was evaluated using the
   Wilcoxon test: ∗p <0.05, ∗∗p <0.01. (E) Violin plot depicting the
   expression of inflammation-related genes (CCL3, CCL4, IL1B, and KLF2)
   according to the response group in cluster 1 (TAM_resting_C1Q^+) and
   cluster 2 (TAM_activating_CCL^+). (F) Dot plot of the KEGG pathway
   enrichment analysis results for the top 100 upregulated genes
   associated with each macrophage-related cluster (clusters 1–4). Dot
   color represents FDR values, whereas circle size indicates the gene
   ratio for each KEGG pathway. (G) Bubble plot of ligand–receptor pairs
   of CCL3 chemokines between myeloid cells and other cell types. CCL, C–C
   motif chemokine ligand; FDR, false discovery rate; KEGG, Kyoto
   Encyclopedia of Genes and Genomes; KLF2, Kruppel-like factor 2; TAM,
   tumor-associated macrophage; UMAP, Uniform Manifold Approximation and
   Projection.

   Analysis of the myeloid cell profiles between responder and
   nonresponder groups revealed a distinct distribution of myeloid cell
   clusters. The responder group had a higher proportion of the
   TAM_activating_CCL^+ cluster (C2), whereas the nonresponder group had
   higher proportions of the resting TAM cluster (C1) and the
   Neutrophil_CSF3R^+ cluster (C4) ([106]Figs. S3B and D). Responders also
   showed significantly elevated expression of proinflammatory mediators,
   including chemokines CCL3, CCL4, the cytokine IL1B, and the
   transcriptional regulator KLF2 ([107]Fig. 3E). Further analysis via
   Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment
   demonstrated activation of immune-recruiting pathways, including
   chemokine signaling and NF-κB signaling, in the TAM_activating_CCL+
   cluster ([108]Fig. 3F). Cell–cell interaction analysis identified
   CCL3–CCR5 as a key ligand–receptor pair between myeloid cells and T
   cells, with stronger interactions observed in responders
   ([109]Fig. 3G). This was further supported by CCL3 expression analysis
   in an independent cohort (GO30140 + IMbrave150, n = 217), where strong
   positive correlations were observed between CCL3 and markers of CD8^+ T
   cells (CD8A; correlation coefficient [R] = 0.71, p = 2.2e-12) and
   cytotoxic T cells (GZMH; Pearson’s r [R] = 0.67, p = 8.4e-11) in
   responders. In nonresponders, weaker correlations were observed between
   CCL3 and these markers (CD8A: R = 0.46, p = 5.2e-05; GZMH: R = 0.38,
   p = 0.0011) ([110]Fig. S3C). mIF analysis demonstrated that responders
   exhibited a high degree of colocalization between CD163^+ and CCL3^+
   cells, whereas this coexistence was significantly reduced in
   nonresponders ([111]Figs. S4A and B). The number of CD163+CCL3^+ TAMs
   was significantly higher in responders than in nonresponders
   ([112]Fig. S4C), and a strong correlation with PFS was observed (R =
   0.998, p <0.0001; [113]Fig. S4D).

cDC2_CD1C^+ dendritic cells are enriched in responders

   Within the DC population, we identified four distinct clusters: (1)
   cDC1_CLEC9A^+ expressing markers such as IDO1, XCR1, and CLEC9A; (2)
   cDC2_CD1C^+ expressing CD1C, CLEC10A, and FCER1A; (3) mDC_LAMP3^+
   representing migrating DCs and expressing markers such as CCR7, BIRC3,
   and LAMP3; and (4) plasmacytoid DCs (pDC) expressing markers IRF7 and
   LILRA4 ([114]Fig. 3B). Upon comparing DC populations, we observed a
   higher proportion of cDC2_CD1C+ in responders than in nonresponders
   ([115]Fig. S3D).

Tumor-associated endothelial cells exhibit enhanced antigen presentation in
responders

   To delineate the distinct features of endothelial cells, we performed
   subcluster analysis on the 1,360 endothelial cells identified. This
   analysis revealed two normal endothelial clusters: Endothelial_PROX1^+
   and liver sinusoidal endothelial cells (LSECs), along with a third
   cluster representing tumor-associated endothelial cells (TECs)
   ([116]Fig. 4A). The Endothelial_PROX1^+ cluster expressed TFF3 and
   PROX1, whereas the LSEC cluster expressed CLEC4G. The TEC cluster
   exhibited high expression of VWF, PLVAP, and PECAM1 ([117]Fig. 4B).

Fig. 4.

   [118]Fig. 4
   [119]Open in a new tab

   Antigen-presenting characteristics in tumor-associated endothelial
   cells.

   (A) UMAP visualization of endothelial cells, color coded by endothelial
   cell subclusters. (B) Expression levels of endothelial marker genes.
   (C) Bar plot of significant GSEA results in GO pathways in TEC. Bar
   color represents the significance of enrichment, measured by the log p
   value. (D) Expression of antigen-presenting related genes, cell
   projection, and angiogenesis based on response in TEC. Statistical
   significance was evaluated using the Student’s t test: ∗∗p <0.01, ∗∗∗p
   <0.001, ∗∗∗∗ p <0.0001. GSEA, gene set enrichment analysis; GO, Gene
   Ontology; LSEC, liver sinusoidal endothelial cells; TEC, tumor
   endothelial cells; UMAP, uniform manifold approximation and projection.

   To further understand the function of TECs in atezolizumab and
   bevacizumab combination therapy, GSEA was performed using ranked DEG
   lists. Among 2,421 significantly enriched pathways (p <0.05), key
   pathways included antigen processing, cell projection, and angiogenesis
   ([120]Fig. 4C; [121]Table S4). Notably, antigen-presenting pathways
   were enriched in the responder group, with high expression of genes
   crucial for antigen presentation, such as CD74, HLA-DRA, B2M, and HLAC
   ([122]Fig. 4D). By contrast, the nonresponders showed enrichment in
   pathways related to cell projection, including processes such as
   assembly, organization, and regulation. Angiogenesis pathways were also
   enriched in nonresponders, particularly those involved in the positive
   regulation of blood vessel endothelial cell proliferation, sprouting
   angiogenesis, and blood vessel maturation. Genes such as EMP1, ICAM1,
   SEPTIN2, and TLN1, known to have roles in cell projection, were highly
   expressed in the nonresponder group. In addition, genes contributing to
   angiogenesis, including CD40, HSPG2, RAMP2, and TGFBR3, exhibited
   higher expression in nonresponders compared with responders.

Integration of bulk RNA-seq data reveals an activating TME-associated gene
signature for predicting response to atezolizumab plus bevacizumab in HCC

   To evaluate the predictive value of gene signatures from single cell
   analysis, we calculated cytotoxicity, exhaustion, activated TAM, and
   cDC2 scores using transcriptomic data from 217 pretreatment tumor
   biopsies (GO30140 + IMbrave150) and assessed associations with PFS and
   OS. Patients with high cytotoxic or exhaustion scores (top 50th
   percentile) had significantly longer PFS under atezolizumab plus
   bevacizumab (p = 0.014 and p = 0.027, respectively; [123]Fig. 5A,B).
   High activated TAM and DC scores were similarly linked to longer PFS
   (p = 0.0039 and p = 0.0022, respectively; [124]Fig. 5C,D). These
   signatures were not predictive for PFS in patients treated with
   sorafenib ([125]Fig. 5E–H). For OS, only the cytotoxic T cell signature
   was significantly associated with this measure (p = 0.024;
   [126]Figs. S5A–D). Based on these findings, we developed a 10-gene
   atezolizumab plus bevacizumab response signature (sc_ABRS) encompassing
   cytotoxic/exhausted T cells, activated TAMs, and cDC2 cells
   ([127]Table S1). High sc_ABRS expression was associated with an
   extended PFS of a median of 10.8 months (95% CI: 7.3 months–not
   reached) compared with 5.5 months (95% CI: 4.0–6.7 months; p <0.001)
   ([128]Fig.5I), but did not correlate with OS (p = 0.29; S5E). With an
   AUC of 0.75 for PFS, sc_ABRS performed comparably to established
   biomarkers, such as CD274 (AUC = 0.72), effective T cell signature
   (AUC = 0.76) and atezolizumab + bevacizumab response signature[129]^20
   (AUC = 0.76) ([130]Fig. 5J). Although sc_ABRS did not outperform these
   markers, its single cell origin provides a more detailed view of the
   tumor immune microenvironment.

Fig. 5.

   [131]Fig. 5
   [132]Open in a new tab

   Predictive value of single-cell derived gene signatures for response to
   atezolizumab plus bevacizumab in HCC.

   (A–D) Kaplan-Meier plots of PFS in 217 patients from the IMBRAVE150 +
   GO30140 trial treated with atezolizumab plus bevacizumab, stratified by
   high versus low (median split) expression of the following gene
   signatures: (A) cytotoxic, (B) exhaustion, (C) activated TAM, and (D)
   cDC2. Patients with higher expression of these signatures showed
   significantly longer PFS. (E–H) Kaplan-Meier plots of PFS in 58
   patients treated with sorafenib in the IMBRAVE150/GO30140 trial,
   showing no significant differences between high and low expression
   groups for the same gene signatures: (E) cytotoxic, (F) exhaustion, (G)
   activated TAM, and (H) cDC2. (I) Violin plots comparing sc_ABRS (a
   10-gene signature derived from single cell analysis) scores across
   different radiological response categories (CR, PR, SD, and PD) based
   on RECIST criteria, with p values from two-sided Student’s t tests. (J)
   Kaplan-Meier plot showing PFS stratified by high versus low sc_ABRS
   scores in patients treated with atezolizumab plus bevacizumab. (K) ROC
   curve analysis demonstrating the predictive performance of the sc_ABRS
   score compared with previously reported biomarkers, including effective
   T cell signatures, CD274, and ABRS. ABRS, atezolizumab + bevacizumab
   response signature; cDC2, type 2 conventional dendritic cells; CR,
   complete response; PFS, progression-free survival; PD, progressive
   disease; PR, partial response; RECIST, Response Evaluation Criteria in
   Solid Tumors; ROC, receiver operating characteristic; sc_ABRS, single
   cell atezolizumab plus bevacizumab response signature; SD, stable
   disease; TAM, tumor-associated macrophage.

Discussion

   In this comprehensive single cell analysis of the immune landscape in
   advanced HCC, we identified immune cell populations and molecular
   signatures associated with favorable outcomes, particularly an
   activated TME characterized by tumor-specific cytotoxic T cells,
   cytokine-enriched TAMs, increased cDC2 cells, and enhanced antigen
   presentation by TECs. This was further supported by findings from an
   independent cohort.

   Consistent with previous findings from the GO30140 and IMbrave150
   trials,[133]^20 we observed that responders to atezolizumab–bevacizumab
   exhibited higher levels of exhausted and cytotoxic gene signatures
   compared with nonresponders. Upregulation of these genes is crucial for
   identifying patients likely to benefit from immune checkpoint inhibitor
   therapy in HCC. By inhibiting the PD-1/PD-L1 axis, these therapies can
   restore exhausted CD8^+ T cells and enhance their effector
   function.[134]^24 Although baseline levels of cytotoxic CD8^+ T cells
   were not significantly different, responders demonstrated enhanced
   clonal expansion, particularly within the C9 cluster of terminally
   differentiated CD8^+ T cells expressing GZMH and CX3CR1. This cluster
   aligns with CD8 TEMRA cells, which have been reported to
   exhibit patrolling activity in the antitumor response under
   atezolizumab–bevacizumab therapy.[135]^5

   Beyond CD8^+ T cells, TAMs have a crucial role in modulating the
   response to atezolizumab and bevacizumab combination therapy. Spatial
   transcriptomic analysis of HCC surgical specimens revealed an
   immunosuppressive TME characterized by SPP1^+ TAMs interacting with
   cancer-associated fibroblasts, which limit immune cell infiltration and
   reduce anti-PD-1 efficacy.[136]^22 However, other studies have shown
   that PD-L1^+ CXCL10^+ TAMs are enriched in responders, suggesting a
   role in recruiting T cells into the TME.[137]^25 These conflicting
   results illustrate the multifaceted role of TAMs within the TME and
   emphasize the need for further investigation into the diverse TAM
   populations and their impact on treatment responses. We found a unique
   TAM subtype characterized by chemokines, such as CXCL2, CCL3, CCL4, and
   IL1B, which was enriched in responders. This TAM_activating_CCL^+
   subtype closely interacts with T cells, suggesting a role in CD8^+ T
   cell recruitment. These features align with the inflammatory
   cytokine-enriched TAM subtype described by Ma et al.,[138]^26
   characterized by expression of inflammatory cytokines, such as IL1B,
   CXCL1/2/3/8, and CCL3, which are involved in active immune cell
   recruitment and regulation during the tumor-associated inflammatory
   response. In addition, this TAM subtype exhibited a high myeloid
   inflammation signature, previously associated with improved PFS in
   patients treated with atezolizumab plus bevacizumab compared with
   atezolizumab monotherapy.[139]^20 However, additional research is
   required to elucidate the function of bevacizumab in modulating these
   TAMs to improve the antitumor response.

   In contrast to the findings of Cappuyns et al.,[140]^25 we observed
   minimal PD-L1 and CXCL10 expression in TAMs ([141]Figs. S6A–C),
   possibly because of limited myeloid cell numbers in our analysis or
   differences in patient characteristics. Our cohort primarily comprised
   advanced viral HCC cases (83.3%) with high-risk features (50%), whereas
   Cappuyns et al.[142]^25 investigated nonviral HCC (93%) with less
   advanced disease (32% with macrovascular invasion). These results
   suggest that TAM interactions in the TME vary by HCC etiology and
   stage, potentially influencing response to atezolizumab plus
   bevacizumab.

   Growing evidence suggests that TECs have a multifaceted role in the
   TME, not only promoting tumor angiogenesis, but also mediating immune
   regulation.[143]^27 In this study, we observed that TECs in responders
   were enriched for antigen-presenting pathways, indicating their active
   involvement in T cell priming and activation. However, we were unable
   to identify specific TEC subpopulations associated with treatment
   response because of the limited number of TECs in the biopsy samples.
   These findings underscore the importance of an immunologically
   activated TME that can effectively recruit and prime T cells,
   contributing to favorable responses to atezolizumab plus bevacizumab in
   advanced HCC.

   This study has several limitations. First, although RNA scope provided
   valuable single cell gene expression data, it lacked spatial resolution
   and protein quantification. Second, the small sample size restricted
   analysis of non-immune cells and differentiation between viral and
   nonviral HCC. In addition, we could not validate findings with other
   scRNA-seq datasets because of time and resource constraints, requiring
   cautious interpretation of our results. Third, although radiological
   reviews confirmed that most biopsies were representative, they
   inherently capture only a small portion of the tumor, potentially
   underestimating intratumor heterogeneity. Lastly, longitudinal biopsies
   were not feasible because of bleeding risks in patients with HCC.
   Despite these limitations, our findings offer a detailed
   transcriptional view of advanced HCC treated with atezolizumab plus
   bevacizumab, supported by bulk transcriptome data from an independent
   cohort.

   In conclusion, single cell analysis of the TME in advanced HCC reveals
   that an activated immune microenvironment before treatment is crucial
   for a favorable response to combination therapy with atezolizumab plus
   bevacizumab. These insights highlight the potential for identifying
   patients with advanced HCC most likely to benefit from such treatment;
   however, further research with larger, longitudinal cohorts is needed
   to explore dynamic changes in the TME.

Abbreviations

   Atezolizumab plus bevacizumab, atezolizumab and bevacizumab combination
   therapy; BCLC, Barcelona Clinic Liver Cancer; CCL, C–C motif chemokine
   ligand; CK, pan-cytokeratin; CMV, cytomegalovirus.CR, complete
   response; DC, dendritic cell; DEG, differentially expressed gene; EBV,
   Epstein–Barr virus; FDR, false discovery rate; FFPE, Formalin-fixed
   paraffin-embedded; GEM, gel bead in emulsion; GO, gene ontology; GSEA,
   gene set enrichment analysis; HCC, hepatocellular carcinoma; KEGG,
   Kyoto Encyclopedia of Genes and Genomes; KLF2, Kruppel-like factor 2;
   LSEC, liver sinusoidal endothelial cell; mIF, Multiplex
   immunofluorescence; mRECIST, modified Response Evaluation Criteria in
   Solid Tumors; NES, Normalized Enrichment Score; NK, natural killer; OS,
   overall survival; PCA, principial components analysis; PD-L1,
   programmed death-ligand 1; PD, progressive disease; pDC, plasmacytoid
   dendritic cell; PFS, progression-free survival; PR, partial response;
   ROI, regions of interest; ScRNA-seq, single-cell RNA sequencing;
   scTCR-seq, single-cell T cell receptor sequencing; SD, stable disease;
   SNN, shared nearest neighbor; TAM, tumor-associated macrophage; Tcm,
   central memory T cell; TCR, T cell receptor; Tcyt, cytotoxic T cell;
   TEC, tumor-associated endothelial cell; Tem, effector memory T cell;
   TIL, tumor-infiltrating lymphocyte; TME, tumor microenvironment; TRM,
   tissue-resident memory T cell; UMAP, uniform manifold approximation and
   projection; UMI, unique molecular identifier.

Financial support statement

   This research was supported by a grant of the Korea Health Technology
   R&D Project through the Korea Health Industry Development Institute
   (KHIDI), funded by the Ministry of Health & Welfare, Republic of Korea
   (grant number:RS-2020-[144]KH088687).

Authors’ contributions

   Conceptualization: Y-HP and WP. Methodology: JL, W-YP, Y-HP, MJG, BGS,
   SYH, IH. Formal analyses: JL, MJG, WK, DIC, KG, SYH, IH. Resources: JL,
   DHS, G-YG, DIC, KG. Writing - original draft: JL, MJG. Writing - review
   & editing: JL, MJG, Y-HP, W-YP, BGS, DHS, WK, G-YG, DIC, KG, SYH, IH.
   Supervision: Y-HP, W-YP, JHL, MSC. Project administration: Y-HP, W-YP,
   JHL, MSC. Funding acquisition: Y-HP, WK. All authors read and approved
   the final manuscript.

Data availability statement

   The datasets used in this study are available from the corresponding
   authors upon reasonable request. All datasets and software used for
   analysis are listed in the CTAT table. Further inquiries can be
   directed to the corresponding authors.

Conflicts of interest

   The authors declare no competing interests.

   Please refer to the accompanying ICMJE disclosure forms for further
   details.

Footnotes

   Author names in bold designate shared co-first authorship

   Supplementary data to this article can be found online at
   [145]https://doi.org/10.1016/j.jhepr.2024.101304.

Contributor Information

   Woong-Yang Park, Email: woongyang.park@samsung.com.

   Yong-Han Paik, Email: yh.paik@samsung.com.

Supplementary data

   The following are the supplementary data to this article:
   Multimedia component 1
   [146]mmc1.pdf^ (1.5MB, pdf)
   Multimedia component 2
   [147]mmc2.docx^ (39KB, docx)
   Multimedia component 3
   [148]mmc3.pdf^ (619.5KB, pdf)
   Multimedia component 4
   [149]mmc4.pdf^ (5.9MB, pdf)

References