Introduction

Hepatocellular carcinoma (HCC) is a prevalent and deadly malignancy that is diagnosed at an advanced stage and has a poor 5-year survival rate.1 2 Recurrence or metastasis still occurs within 5 years in the majority of patients even after surgical resection.3 The application of immune checkpoint blockade (ICB), such as PD-1/L1 inhibitors and CTLA4 inhibitors, has led to breakthroughs in the treatment of liver cancer patients.4–6 However, only a small proportion of patients can benefit from these treatments, indicating significant challenges in the field of immunotherapy for HCC.7 8 Cancer-associated fibroblasts (CAFs), which are critical components of the tumor microenvironment (TME), promote immune evasion, tumor metastasis, and therapy resistance by remodeling the extracellular matrix (ECM), secreting growth factor and cytokines.9 10 Genetic depletion or pharmacological targeting of distinct CAF populations leads to different outcomes,11 12 highlighting the CAFs are highly heterogeneous in types and functions affecting cancer progression and immunomodulation.

CAF populations are diverse and dynamic and originate from various cell sources, which contributes to their functional diversity.13 Understanding the heterogeneity and subpopulation characteristics of CAFs is essential for targeting protumorigenic CAFs or interference with their activity as potential strategies in anticancer therapy. Recent advancements in techniques such as single-cell sequencing and flow cytometry have enabled precise analysis of CAF subtypes and their biological functions in different cancer types.13–16 These studies provide evidence supporting the presence of two distinct functional populations of CAFs: procancer CAFs and anticancer CAFs. Procancer CAFs orchestrate ECM remodeling and tumor-promoting inflammation, and modulate the immune microenvironment toward immunosuppression,17 such as myofibroblastic CAFs and inflammatory CAFs (iCAFs). Conversely, anticancer CAFs are able to stimulate T-cell activation associated with MHC class II expressing to improve immunotherapeutic strategies, such as antigen-presenting CAFs (apCAFs).18 CAFs are highly interrelated with sensitivity to anticancer therapies. Neoadjuvant chemotherapy can remodel CAFs to enhance immune activation and inhibit tumor progression.11 However, CAFs can also remodel the ECM, making it denser and potentially forming a physical barrier that impedes the delivery of chemotherapy and immunotherapy drugs.19 Overall, various CAF subtypes may activate distinct molecular pathways to influence the TME and resistance to immunotherapy and chemotherapy.

However, the intricate crosstalk between CAFs, HCC cells, and immune cells, particularly within different CAF subtypes, remains unclear. A comprehensive understanding of the multidimensional interactions between CAFs and infiltrating immune cells in the TME would help us elucidate the mechanisms through which CAFs induce immunosuppression. Further exploration of these interactions may reveal potential molecular targets for CAF-targeted therapy, providing new strategies for subsequent immunotherapy.

We integrated single-cell RNA sequencing (scRNA-seq) datasets for HCC tumors. Our analysis led to the identification of six common subtypes of CAFs in HCC tumors: STMN1+ CAFs, CXCL12+ CAFs, MYH11+ CAFs, SEPT7+ CAFs, POSTN+ CAFs, and CD36+ CAFs. Among these subtypes, we specifically focused on the newly discovered POSTN+ CAF subpopulation. POSTN+ CAFs play a crucial role in promoting the progression of HCC through the activation of the ECM, hypoxia and transforming growth factor beta (TGF-β) signaling pathways. Furthermore, they contribute to the formation of an immune barrier that suppresses the infiltration of CD8+ T cells into the TME. Additionally, POSTN+ CAFs interact with SPP1+ macrophages via the IL-6/STAT3 axis. Our study revealed that the abundance of POSTN+ CAFs can serve as a reliable predictor of the response to immunotherapy in HCC patients. This indicates that targeting POSTN+ CAFs could enhance the effectiveness of immunotherapeutic approaches for treating HCC. In summary, our study elucidated the novel intercellular communication network between POSTN+ CAFs and SPP1+ macrophages and between POSTN+ CAFs and CD8+ T cells. These findings suggest that POSTN+ CAFs play a significant role in promoting the development of an immunosuppressive TME.

Material and methodsHuman liver samples

Liver samples were obtained from patients who underwent hepatectomy or percutaneous liver biopsy. Clinical information is summarized in online supplemental table 1.

Cell culture

The human myeloid leukemia mononuclear cell line (THP-1) was obtained from the Cell Resource Center at the Chinese Academy of Medical Sciences and authenticated by the same. THP-1 cells were cultured in RPMI1640 supplemented with 10% FBS (Gibco) and 1% penicillin-streptomycin (Invitrogen, USA). To conduct differentiation assays, THP-1 cells were differentiated into macrophages through 24-hour culturing with 100 ng/mL phorbol myristate acetate (Sigma). CAFs transfected with si-CREB3L1 or si-NC for 48 hours were seeded in six-well plates and cultured under hypoxia or normoxia for 24 hours. siRNAs targeting CREB3L1 are listed in online supplemental table 2.

Plasmid transfections

The pcDNA3-based plasmid encoding human POSTN was acquired from Gentlegen (Suzhou, China). CAFs were seeded and allowed to grow in 24-well culture plates for 48 hours before transfection with the plasmids using Lipofectamine 3000 (Invitrogen).

Chemotaxis assay

Macrophages at a concentration of 4×104 cells per well were plated in the upper chamber with 0.1% gelatin-coated membranes, while CAFs were plated in the lower chamber at a concentration of 1×105 cells per well. Tocilizumab (200 µg/mL) was added to the upper chamber. After incubating for 24 hours at 37°C, the cells were washed with 1x PBS, fixed in 10% paraformaldehyde for 30 min, and stained with 0.05% crystal violet for another 30 min. Finally, the number of cells that had migrated to the bottom of the membrane was assessed by counting them in four different regions.

Protein extraction and western blot

The cells were lysed using RIPA solution (Beyotime Biotechnology, China) supplemented with protease and phosphatase inhibitors (Beyotime Biotechnology). Western blot analysis was conducted following previously described protocols.20 The protein concentration in cell lysates was determined using a BCA assay (Beyotime Biotechnology). Subsequently, the cell lysates were subjected to 10% SDS-polyacrylamide gel electrophoresis and transferred onto a polyvinylidene fluoride membrane (Bio-Rad, Hercules, CA) for antibody incubation.

Antibodies

Primary antibodies used in this study were as follows: rabbit anti-POSTN mAb (ab215199, Abcam); rabbit anti-CD8 mAb (245118, Abcam); rabbit anti-MYH11 mAb (ab133567, Abcam); rabbit anti-STMN1 mAb (ab52630, Abcam); rabbit anti-CD36 mAb (ab252923, Abcam); rabbit anti-α-SMA mAb (19 245S, Cell Signaling Technology); rabbit anti-STAT3 mAb (A1192, abclonal Technology); rabbit anti-IL-6 mAb (A1570, abclonal Technology); rabbit anti-SPP1 mAb (8242S, Cell Signaling Technology); rabbit anti-p-STAT3 mAb (9145S, Cell Signaling Technology); rabbit anti-Osteopontin mAb (A1499, abclonal Technology); rat anti-CD68 mAb (ab53444, Abcam); rabbit anti-Osteopontin mAb (ab214050, Abcam); rabbit anti-α-SMA mAb (19 245S, Cell Signaling Technology); rabbit anti-CREB3L1 mAb (11 235–2-AP, Proteintech); mouse anti-β-actin mAb (3700S, Cell Signaling Technology).

Secondary antibodies used in this study were as follows: HRP-conjugated goat anti-rabbit IgG (ab205718, Abcam); HRP-conjugated goat anti-rat IgG (ab205720, Abcam); Alexa Fluor 488-conjugated goat anti-rabbit IgG (A-11008, Invitrogen); Alexa Fluor 594-conjugated goat anti-mouse IgG (A-11005, Invitrogen).

Quantitative real-time PCR

Total RNA was extracted from the samples using TRIzol reagent (Invitrogen) and then precipitated with isopropyl alcohol. The concentration and quality of the RNA were assessed using a NanoDrop system. Complementary DNA was synthesized using a commercially available kit (Vazyme, Nanjing, China) following the manufacturer’s instructions. Real-time quantitative PCR was performed using SYBR-green-based kits (Vazyme) to quantify mRNA levels. The relative expression of the target gene was normalized to β-actin in each sample. The following primer sequences were used: β-actin (forward, 5′-CTCGCCTTTGCCGATCC-3′; reverse, 5′-ATCCTTCTGACCCATGCCC-3′), POSTN (forward, 5′-CTCATAGTCGTATCAGGGGTCG-3′; reverse, 5′-ACACAGTCGTTTTCTGTCCAC-3′), IL-6 (forward, 5′-TGAGGAGACTTGCCTGGTGAA-3′; reverse, 5′-CAGCTCTGGCTTGTTCCTCAC-3′), CREB3L1 (forward, 5′-GGAGAATGCCAACAGGACC-3′; reverse, 5′-GCACCAGAACAAAGCACAAG-3′), COL1A1 (forward, 5′-CCCCTGGAAAGAATGGAGATG-3′; reverse, 5′-TCCAAACCACTGAAACCTCTG-3′), COL3A1 (forward, 5′-AGGAAATGATGGTGCTCCTG-3′; reverse, 5′-GTTCCCCAGGTTTTCCATTT-3′), COL5A1 (forward, 5′-TCGCTTACAGAGTCACCAAAG-3′; reverse, 5′-GTTGTAGATGGAGACCAGGAAG-3′).

Animal studies

All mice were housed in an SPF level animal facility with controlled temperature and humidity on a 12–12 hours light-dark cycle with food and water supplied ad libitum.

Adeno-associated virus (AAV)-mediated gene transfer for knockdown POSTN in vivo, the following AAV was used in this study: AAV8-control and AAV-shPOSTN (Shanghai Genechem). AAV8-shPOSTN was injected into 6 weeks old male C57BL/6 mice via the tail vein at a dose of 1×1099 infectious units (IFU) per 200 µL per mouse 2 weeks before orthotopic tumor transplantation. For orthotopic implantation, 3×106 Hepa1-6 cells were injected into the left lobes of the livers of C57BL/6 mice. Mice were treated with either 10 mg/kg anti‐mouse PD‐1 (BE0146, BioXCell) or isotype IgG control (BE0089, BioXCell) through intraperitoneal injection twice weekly starting on day 7. Mice were sacrificed on day 21. shRNAs targeting POSTN are listed in online supplemental table 2.

Tissue processing and fibroblast isolation

Primary HCC tissues were minced using a scalpel in a tissue culture dish. Subsequently, the minced tissues were enzymatically dissociated in 10 mL of PBS containing 0.1% collagenase I at 37°C for 1 hour with gentle agitation. The resulting suspension was neutralized with a complete medium and then centrifuged at 300×g for 5 min. The collected cells were plated in culture dishes and allowed to grow. After 48 hours, the non-adherent cells and tissue debris were removed by washing the dishes twice with PBS. The adherent fibroblasts were further incubated for 6–10 days until they reached 80%–90% confluence.

Flow cytometry

Briefly, the cells were resuspended in the staining buffer and incubated with the following antibodies on ice for 30 min: anti-CD45 (368521, Biolegend, USA), anti-EMCAM (324221, Biolegend), Live/Dead blue fluorescent dye (L34963, Thermo Scientific) and anti-CD29 (303023, Biolegend). Among live cells, EpCAM and CD45 were used among live cells to eliminate epithelial and immune cells, respectively, while CD29 was indicative of stromal cells. For intracellular staining, the cells were fixed with fixation buffer (Biolegend) on ice for 15 min and then washed twice with an Intracellular Staining Permeabilization Wash Buffer. Antibodies against POSTN (sc-398631, Santa Cruz) were added and incubated for 1 hour on ice.

Multiplex immunohistochemistry staining

Multiplex immunohistochemistry (mIHC) staining was employed to visualize the expression of α-SMA, STMN1, MYH11, CD36, POSTN, CD8, and CD68 in tumor tissues. Paraffin blocks were used to obtain three consecutive sections, each 3 µm thick, with one section designated for H&E staining. The remaining two FFPE tumor slides, also 3 µm thick, underwent a 3-hour dehydration process at 70°C. Subsequently, the paraffin sections were then de-paraffinized using xylene and rehydrated with alcohol. Heat-induced antigen retrieval was performed in EDTA buffer at pH 9.0 using a microwave oven. The sections were then blocked with a commercially available blocking buffer for 10 min. Following this, the slides were sequentially incubated with primary antibodies and horseradish peroxidase-conjugated secondary antibodies, followed by tyramide signal amplification (TSA). After each TSA round, the slides underwent antigen retrieval and antibody stripping. Finally, nuclei were stained with DAPI. This comprehensive procedure facilitated the visualization of multiple antigens in the tumor tissues, thereby enabling the assessment of expression patterns for the targeted proteins.

Data collection

The public scRNA-seq datasets used in this research were obtained from various sources. Specifically, the datasets GSE149614,21 GSE151530,22 GSE156625,23 GSE189903,24 and GSE20264225 were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Additionally, dataset HRA00174826 was obtained from the Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA001748). For the analysis of normalized gene expression data and microarray datasets of HCC samples, publicly available data from The Cancer Genome Atlas (TCGA) data portal (http://gdac.broadinstitute.org/) and GEO were employed. The specific datasets accessed from GEO include GSE10143,27 GSE192912,28 GSE109211,29 GSE10186,30 and GSE54236.31 To acquire the expression matrix and clinical information of the TCGA-LIHC dataset and GSE14520,32 the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse14520) was accessed. Furthermore, expression profile and clinical information from the Fudan HBV-HCC cohort33 were collected and collated from the biosino NODE database (OEP000321,33 https://www.biosino.org/node/project/detail/OEP000321).

Spatial transcriptome sequencing data of 13 samples including HCC tumors, normal tissues, and leading tissues were acquired from the study conducted by Wu et al34 (http://lifeome.net/supp/livercancer-st/data.htm). Furthermore, spatial transcriptome data from tumor sections of eight HCC patients who underwent anti-PD-1 treatment (five non-responders and three responders) were downloaded from the study conducted by Liu et al35 (https://data.mendeley.com/datasets/skrx2fz79n/1).

Single-cell transcriptome analysis and differential expression analysis

CellRanger (V.6.0.2) was used for read mapping and gene expression quantification. Seurat package (V.4.1.1) was employed for downstream analysis. Cells with fewer than 1000 UMIs or greater than 15% mitochondria genes were excluded from the analysis. Doublets were assessed using the DoubletFinder (V.2.0.3) algorithm each sample. To address batch effects and prevent the analysis from being dominated by individual patient characteristics, the Harmony algorithm was applied to remove batch effects between samples in the single-cell data. To perform dimensionality reduction and identify cell subtypes, the top 4000 most variable genes were selected using the FindVariableFeatures function and the data were scaled accordingly. Principal component analysis (PCA) was conducted using these variable genes. Nearest neighbors for graph clustering based on the top 50 principal components were determined using the FindNeighbors function. The resulting clusters were obtained using the FindCluster function, and cell visualization was achieved using the uniform manifold approximation and projection (UMAP) algorithm. Furthermore, gene signatures specific to various cell types were scored in the identified clusters. These gene signatures included markers for B cells (CD79A and CD79B), endothelial cells (PECAM1 and VWF), hepatocytes (ALB and APOA2), fibroblasts (COL1A1 and COL1A2), myeloid cells (CSF3R, C1QB, CD1C, and CLEC9A), plasma cells (MZB1 and IGHG1), and T/NK cells (CD3D, CD3E, KLRD1, and KLRC1). Differential gene expression analysis between clusters was performed using the “FindAllMarkers” function. The parameters used for this analysis were min.pct=0.15, logfc.threshold=0.15, and only.pos=TRUE. The Wilcoxon’s rank-sum test with the Benjamini-Hochberg method was employed to obtain p values and adjusted p values for the comparisons.

Gene sets level analysis for scRNA-seq data

Single-cell signature scoring was conducted employing the RunAUCell function in pochi R package (V.0.1.0). Differential signature score enrichment between groups was determined using a two-sided Wilcoxon rank-sum test with Benjamini-Hochberg FDR correction. The gene sets used for this analysis were provided in online supplemental table 3.

Pathway and cytokine signaling signatures analysis

We performed pathway and cytokine signaling analyses on CAF subgroups using PROGENy and CytoSig, respectively. Scores were computed using the run_aucell function from decoupleR package (V.2.2.2) based on PROGENy network model and Cytosig python packages with default parameters. The top 1000 target genes of the progeny model were used, as recommended for single-cell data.

Function enrichment and gene set enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment and gene set enrichment analysis (GSEA) of differentially expressed gene sets were implemented in the clusterProfiler R package (V.4.7.1.2). Enriched pathways with adjusted p value below 0.05 were considered as significantly enriched by DEGs.

Trajectory analysis

To better understand the cellular differentiation among various subtypes of CAFs, we used the Monocle 2 algorithm V.2.24.1. Initially, we honed in on the top 100 differentially expressed genes distinguishing between the CAF subtypes. Employing the reduceDimension function with the DDRtree method and limiting it to a maximum of two components facilitated dimensionality reduction. Subsequently, we used the “differentialGeneTest” function to identify genes exhibiting significant differential expression correlated with pseudotime values. Additionally, we used branch expression analysis modeling (BEAM) to unravel gene expression patterns influenced by branch fate. For visualizing the expression dynamics of branch-specific genes along the differentiation trajectories, we employed the plot_genes_branched_pseudotime function, generating curves that depicted the Loess-smoothed expression profiles of genes along each trajectory.

Transcription factor regulon analysis

We used pySCENIC for the analysis of the regulatory network and regulon activity. Through pySCENIC analysis, we identified subgroup specific transcription factors (TFs) using the Wilcoxon rank-sum test. Additionally, we generated regulons-associated specific score (RSS) within each specific cell type by computing Jensen-Shannon divergence, facilitated by philentropy (V.0.6.0) package.

Cell–cell communications

To explore potential interactions among different cell types within the HCC TME, we conducted cell–cell communication analysis using the CellPhoneDB Python package (V.2.1.1). This approach relies on a publicly available repository of curated receptors and ligands along with their interactions. Enriched receptor-ligand pairs between two cell types were predicted based on the expression of a receptor by one cell type and the corresponding ligand by another. By default, we considered only ligands and receptors expressed in at least 10% of cells within a given cell type. Additionally, we selected significant pairs with a p value <0.05 and a mean value ≥1 for subsequent analyses.

To further identify the key mediators of two cell subgroups, we applied Nichenet package to infer the interaction between POSTN+ CAFs and SPP1+ macrophages (Mphs). For ligands and receptor interactions, genes which are expressed in larger than 10% cells of clusters were considered. Top 20 ligands and top 100 targets of differential expressed genes of “sender cells” and “receive cells,” were extracted for paired ligand-receptor activity analysis. We used the active_ligand_target_links function to compute the potential intensity of regulation between the ligand and target. Activity scores ranged from 0 to 1. The expression of differential expressed ligands and receptors were also shown in heatmap by calculating the average genes expression in indicated cell types and scaled across indicated subtypes.

Cell types deconvolution for RNA-seq datasets

To assess the function of cell types in larger compendiums of samples, we applied the online tool CIBERSORTx to generate a reference signature matrix from our single-cell RNA-seq dataset and estimate cell-type proportions from the public bulk RNA-seq and microarray datasets based constructed cell-type reference. For creating signature matrices, CIBERSORTx was run with quartile normalization was disabled for RNA-seq datasets and was enabled for microarray datasets, and all other parameters with their default settings. For the imputation of cell fractions, the quartile normalization was disabled for RNA-seq datasets and was enabled for microarray datasets, the permutation parameter was set to 500 times, and all other parameters were kept at their default settings. Spearman’s correlation analysis was performed to assess the relationship among the proportions of cell-type infiltration and was plotted by corrplot package (V.0.92).

Survival analysis

Survival analysis was performed by R package survival. HR was calculated by Cox proportional hazards model and 95% CI was reported, and Kaplan-Meier survival curve was modeled by the survfit function. The two-sided long-rank test was used to compare Kaplan-Meier survival curves. The comparison of the percentage of patients who respond to ICB treatment between different groups was determined by the χ2 test.

Spatial transcriptomics analysis

The gene-spot matrices generated after ST data processing from ST and Visium samples were analyzed with the Seurat package (V.4.1.1) in R. Spots were filtered for minimum detected gene count of 200 genes while genes with fewer than 10 read counts or expressed in fewer than three spots were removed. Normalization across spots was performed with the LogVMR function. Dimensionality reduction and clustering were performed with independent component analysis (PCA) with the first 30 PCs. To better exhibit spatial expression of features, the spots were enhanced using the “spatialEnhance” function of the BayesSpace package (V.1.6.0), the expression features were enhanced with the “enhanceFeatures” function. The signature score derived from the scRNA-seq dataset was added to the “metadata” of the ST dataset with the “AddModulScore” function with default parameters in Seurat. Spatial feature expression plots were generated with the “SpatialFeaturePlot” function in the Seurat package. We applied the SpaGene method to identify ten spatial gene expression patterns for each sample and assessed the similarity between patterns based on the Jaccard index. The result was shown using the ComplexHeatmap package (V.2.12.1).

Statistical analysis

Statistical analyses were performed as described in Seurat (V.4.1.1) in R (V.4.2.0) software. Spearman correlation was used to estimate correlations. Survival analysis was measured using the Kaplan-Meier method. Statistical significance was determined by the Kruskal-Wallis test and Wilcoxon’s rank-sum test with the Benjamini-Hochberg method for multiple comparison correction.

ResultsGeneration of an HCC single-cell transcriptomic atlas

To systematically investigate the cellular composition of HCC, we compiled an HCC atlas by compiling scRNA-seq datasets from six datasets for a total of 220 samples (figure 1A and online supplemental figure S1B–C). This comprehensive scRNA-seq atlas of HCC incorporates quality-controlled and preanalyzed transcriptomic data from publicly available datasets on tumor samples and adjacent normal tissues. We employed the Harmony algorithm to correct for batch effects in scRNA-seq datasets from tumor and adjacent normal tissues. Overall, we retained nearly 1.1 million cell transcriptomes for subsequent analysis, and the cells were clustered into 32 distinct clusters representing seven major cell subtypes. These subtypes were identified based on classical single-cell marker expression. Specifically, we identified 25,946 B cells marked by CD79A and CD79B, 90,809 endothelial cells marked by PECAM1 and VWF, 137,809 hepatocytes marked by ALB and APOA2, 39,230 fibroblasts marked by COL1A1 and COL1A2, 212,308 myeloid cells marked by LYZ and C1QB, 15,986 plasma cells marked by MZB1 and IGHG1, and 579,944 T/NK cells marked by CD3D and CD3E (figure 1B–E and online supplemental figure S1A). Additionally, a subset of myeloid cells expressed CSF3R (neutrophil marker), CD1C (cDC2 marker), and CLEC9A (cDC1 marker), while a subset of T/NK cells expressed KLRD1 and KLRC1 (NK cell markers) (figure 1D). Although all seven major cell types were present in both tumor and adjacent normal tissues, their infiltration levels varied, possibly reflecting different stages of HCC progression. Remarkably, in addition to hepatocytes, a greater proportion of fibroblasts originated from tumor tissue to other liver cells (figure 1E), suggesting a potential role for CAFs in the progression of HCC.

Figure 1

Panoramic scRNA-seq analysis of liver samples of hepatocellular carcinoma patients. (A) Schematic representation of the data collection process. (B) Uniform manifold approximation and projection (UMAP) plot showing the seven major cell types. Dots represent individual cells, and colors represent different cell populations. (C) Dot plot showing the top two differentially expressed genes of major cell types. (D) UMAP plots showing the expression of classical molecular markers of corresponding cell types; the color represents the marker expression value. (E) Proportion of the sample types in seven major cell types shown in bar plots (left panel) and the total cell number of each cell type (right panel) are shown.

Identification and distinct features of six CAF subtypes in human HCC

The heterogeneity and plasticity of CAFs play pivotal roles in tumor progression. In our study, we identified six subpopulations of CAFs based on their gene expression profiles: CD36+ CAFs, CXC12+ CAFs, MYH11+ CAFs, POSTN+ CAFs, SEPT7+ CAFs, and STMN1+ CAFs (figure 2A,D and online supplemental figure S2A). CD36+ CAFs, POSTN+ CAFs, and STMN1+ CAFs were significantly more abundant in tumors than in non-tumor tissues, while the opposite trend was observed for CXC12+ CAFs, MYH11+ CAFs, and SEPT7+ CAFs (figure 2B,C and online supplemental figure S2B). In addition to changes in cell abundance, scRNA-seq revealed distinct cellular characteristics of different CAF subtypes. According to the literature, we reclassified CAFs from HCC patients into iCAFs, myofibroblast-like CAFs (myCAFs), apCAFs, vascular CAFs (vCAFs), lipid process-associated CAFs (lpCAFs), and proliferative CAFs (pCAFs) based on their molecular features. POSTN+ CAFs exhibited characteristics of iCAFs and myCAFs, as they expressed inflammatory signature genes (CXCL9, CXC10, and CXCL11) and ECM-associated signature genes (POSTN, MMP14, and MMP11). Notably, POSTN+ CAFs in tumor tissues expressed higher levels of inflammatory and ECM signatures than those in normal tissues. MYH11+ CAFs were identified as vCAFs due to their high expression of vascular-related signatures (MYH11, MUSTN1, and DSTN). CXCL12+ CAFs were designated apCAFs and exhibited high levels of antigen presentation signature genes (CD74, CXCL12, and HLA-DRB1). Interestingly, CXCL12+ CAFs in tumor tissues showed lower levels of antigen presentation signature genes than those in normal tissues, suggesting a potential tumor-suppressive phenotype. CD36+ CAFs and SEPT7+ CAFs were categorized as lpCAFs due to their high expression of lipid metabolism-related genes (APOC3, APOC1, and FABP1). In tumor tissues, these CAF subtypes showed even higher levels of lipid process signature genes than did normal tissues. Finally, the STMN1+ CAFs were classified as pCAFs because they exhibited high expression of genes associated with cell proliferation (STMN1, TOP2A, and MKI67). This finding was consistent with the results of the cell cycle analysis (figure 2E,F and online supplemental figure S2B–D).

Figure 2

Heterogeneity and plasticity analysis of cancer-associated fibroblasts in hepatocellular carcinoma (HCC). (A) Uniform manifold approximation and projection plot showing the classification of cancer-associated fibroblasts (CAFs), with different colors representing different subtypes. (B) Proportion of the sample types in six CAF subtypes shown in bar plots, with different colors representing different sample types. (C) Sankey plot showing variation in the proportions of different CAF subtypes between the sample types, with different colors representing different CAF subtypes. (D) Dot plot showing the expression of the top three different expressed genes in the six CAF subtypes. The color depth represents the average expression value and the size of the point represents the percentage of the gene expression. The Wilcoxon rank-sum test was used to assess differences. (E) Radar plot showing the mean expression of six CAF subgroup signatures calculated by the AUCell algorithm in the reclassified CAF subtypes. (F) Heatmap showing the mean expression of associated gene signatures in each CAF subtype. (G) Kyoto Encyclopedia of Genes and Genomes terms of differentially expressed genes significantly enriched in each CAF subtype. The colors represent the scaled value of the −log10 p value. (H) Differential activation of (up) PROGENy pathways and (bottom) CytoSig cytokine signaling signatures in each CAF subtype. The heatmap colors indicate the deviation from the overall mean. The black dots indicate significance at different p value thresholds. Only cytokine signatures with a p value ≤0.05 in at least one group are shown. (I) Representative IF staining of human HCC tissue. The levels of α-SMA (red), STMN1 (silver), CD36 (purple), POSTN (green), and MYH11 (yellow) in individual and merged channels are shown. Bar, 100 µm. (J) Box plots showing the absolute infiltration of each CAF subtype in tumor and normal samples in GSE10143. (K) Box plots showing the ssGSEA score of each CAF subtype in paired CAFs (n=9) and para-cancer fibroblasts (PAFs, n=9) in GSE192912. apCAFs, antigen-presenting CAFs; ECM, extracellular matrix; iCAFs, inflammatory CAFs; lpCAFs, lipid process-associated CAFs; myCAFs, myofibroblast-like CAFs; pCAFs, proliferative CAFs; vCAFs, vascular CAFs.

To investigate the functionality of different CAF subpopulations, we conducted KEGG functional enrichment analysis on each subpopulation. Notably, KEGG analyses revealed significant enrichment of the PPAR signaling pathway and cholesterol metabolism pathway in CD36+ CAFs and SEPT7+ CAFs, consistent with the findings of the lpCAF signature. Moreover, pathways related to ECM-receptor interactions, focal adhesion, and PI3K-Akt signaling were found to be significantly enriched in POSTN+ CAFs. In line with the vCAF signature, MYH11+ CAFs were significantly enriched in pathways associated with vascular smooth muscle contraction, calcium signaling, and cardiac muscle contraction. CXCL12+ CAFs exhibited significant enrichment in pathways related to cytokine–cytokine receptor interactions, complement and coagulation cascades, and viral protein interaction with cytokines and cytokine receptors. Furthermore, STMN1+ CAFs, referred to as pCAFs, were found to be significantly enriched in the cell cycle and DNA replication pathways based on KEGG pathway enrichment analysis (figure 2G and online supplemental figure S2E).

Interestingly, the TGF-β pathway, which has been implicated in CAF-mediated cancer progression involving cell proliferation, invasion, and metastasis, was significantly upregulated in POSTN+ CAFs (figure 2H). Consistently, the CytoSig tool revealed significant upregulation of the TGF-β1 cytokine signature in POSTN+ CAFs, and genes in this signature are known to stimulate growth and ECM production in fibroblasts (figure 2H). This finding suggested that POSTN+ CAFs may play a crucial role in tumor progression and fibrosis through TGF-β signaling. To validate the presence of major CAF subsets in human HCC samples, we performed mIHC staining (figure 2I). Remarkably, POSTN+ CAFs were more highly enriched in the peripheral region of HCC tumors than in the core region (figure 2I and online supplemental figure S2H). This finding was further supported by the analysis of a liver tissue array dataset and fibroblast transcriptome sequencing data, which both confirmed the specific enrichment of POSTN+ CAFs in tumor tissues (figure 2J,K).

Cell-state transition trajectory of different CAF subpopulations

To gain insights into the dynamic processes underlying CAF subtypes at the single-cell level, we derived the pseudotime cell trajectory of the various CAF subtypes based on the Monocle 2 algorithm (figure 3A). Remarkably, statistical analyses of the combination of pseudotime and inferred state data revealed two distinct trajectories, with CD36+ CAFs (lpCAFs) representing the progenitor state and subsequently evolving into POSTN+ CAFs or MYH11+ CAFs separately (figure 3B,C and online supplemental figure S3A). Additionally, the majority of proliferative STMN1+ CAFs were present in the initiation state, while CXCL12+ CAFs were present in the terminal state within the differentiation branch, leading to differentiation of POSTN+ CAFs. Moreover, we explored gene expression patterns associated with the two differentiation branches involved in CAF state transitions. As depicted in the smoothed heatmap, the functional enrichment of highly expressed genes in the early stage of the trajectory (cluster 4) in both differentiation branches was enriched in protein secretion and G2M checkpoint pathways. In contrast to the MYH11+ CAF differentiation branch, the terminal state of the POSTN+ CAF differentiation branch exhibited higher expression levels of cluster 3 genes, which are implicated in epithelial-mesenchymal transition (EMT), TGF-β signaling, and hypoxia pathways known to be strongly associated with tumor progression. Nevertheless, the terminal state of the MYH11+ CAF differentiation branch displayed increased expression of cluster 1 genes, which were enriched in TNF-α signaling via the NF-KB, G2M checkpoint, and adipogenesis pathways (figure 3D).

Figure 3

Trajectory analysis of cancer-associated fibroblast (CAF) subtypes in hepatocellular carcinoma. (A) The developmental trajectory of CAFs was inferred by monocle 2 and colored according to the subtype. (B) The pseudotime of the developmental trajectory. The colors from purple to yellow indicate values from low to high, respectively. (C) The inferred states of the trajectory. The proportions of CAF subtypes are shown in the circles with different colors representing different CAF subtypes, and each circle color represents a different state. (D) Expression heatmap of significant (q<1e-3) genes based on branch expression analysis comparing the two CAF cell fates (left). Box plots showing the top five significantly enriched hallmark pathways in each gene cluster (right). (E) Dot plot showing Spearman’s correlation of the pseudotime and hypoxia scores, with different colors representing different CAF subtypes. (F) Pseudotime projections of transcriptional changes in extracellular matrix-associated genes and hypoxia-associated genes between the two trajectory branches. (G) Heatmap showing the mean activity of the top differentially activated regulons in each CAF subtype inferred by pyscenic. (H) Dot plots showing the top 15 specific activated transcription factors ranked by the regulon-specific score in each CAF subtype. (I) Uniform manifold approximation and projection (UMAP) plot (up) and violin plot (bottom) showing the expression of CREB3L1. (J) UMAP plot (up) and violin plot (bottom) showing the activity of the CREB3L1 regulon.

We further evaluated the correlation between the pseudotime axis and a set of progression-related pathways, including the EMT, TGF-β signaling, and hypoxia pathways. Interestingly, the hypoxia pathway showed a significant positive correlation with pseudotime (Spearman’s correlation test, rho=0.324, p<0.001) and was specifically activated in POSTN+ CAFs (figure 3E and online supplemental figure S3B,C). Hypoxia is a critical characteristic of the TME and plays a role in regulating CAF heterogeneity and promoting tumor progression. To identify potential differentially expressed genes between the two branches, we performed BEAM. We observed that antigen presentation-related genes (CD74 and CXCL12) and cytokines (CCL19, CCL21, and IL32) were significantly upregulated in the POSTN+ CAF differentiation branch (online supplemental figure S3D) and increased expression of ECM-related genes (COL11A1, COL6A1, COL6A3, COL8A1, and VCAN) and hypoxia-associated features (GLI2, NFKB1, STAT1, STAT2, TWIST1, and CREB3L1) in the POSTN+ CAF differentiation branch (figure 3F). Moreover, to assess the spatial organization of POSTN+ CAFs and hypoxic regions, spatial transcriptomics analysis of tumor-derived sections from two HCC patients was performed. The results revealed colocalization of POSTN+ CAFs and malignant cells in hypoxic zones within the same spot (online supplemental figure S3E). This finding suggested that hypoxia might stimulate CAF-mediated collagen secretion to modulate ECM properties.

To further determine the master regulator of POSTN+ CAFs, we conducted pySCENIC analysis to evaluate the most highly expressed TFs and their activities in the TF regulatory network. Specifically, we observed that CREB3L1 exhibited high expression and activity levels in the regulatory network of POSTN+ CAFs, and we constructed an RSS model for the POSTN+ CAF cell type using Jensen-Shannon divergence and identified CREB3L1 as one of the 15 most prevalent TFs based on the RSS, suggesting that CREB3L1 may serve as the key TF driving this differentiation pathway (figure 3G–J). Similarly, we found that hypoxia stimulated CAFs expressing POSTN and CREB3L1 (online supplemental figure S3F). Additionally, we next observed that the protein level of CREB3L1 was significantly greater in OE-POSTN-CAFs than in OE-NC-CAFs and that POSTN expression was strongly correlated with CREB3L1 expression in the TCGA-LIHC cohort (online supplemental figure S3G,H).

A previous study demonstrated that CREB3L1 is highly correlated with ECM production,36 which facilitates remodeling of the tumor stromal microenvironment. We next detected the mRNA levels of ECM-related genes in CAFs treated with si-NC or si-CREB3L1 under hypoxia, and the results showed that CREB3L1 promoted collagen secretion to modulate ECM properties (online supplemental figure S3I). These findings suggested that CREB3L1 may be a key molecule driving this differentiation pathway through its regulation of downstream signaling pathways.

Tumor-specific POSTN+ CAFs are associated with HCC progression

The differences in the proportions of infiltrated cell types between tumor-adjacent tissues and tumor tissues indicate that the TME undergoes dynamic remodeling, which plays a crucial role in HCC progression. CAFs have long been suggested to be a critical stromal cell type involved in the regulation of tumorigenesis and cancer progression. However, the heterogeneous nature of CAF populations has not been fully elucidated. To further investigate this possibility, we used publicly available data on HCC to validate the presence of tumor-specific POSTN+ CAFs. We compared the differential infiltration of POSTN-expressing fibroblasts between tumor and adjacent normal tissues from each donor. Tissue immunofluorescence staining revealed significantly greater POSTN expression in tumor tissues than in non-tumor tissues (online supplemental figure 4A). This phenomenon was also verified at the cellular level through immunofluorescence analysis and flow cytometry (online supplemental figure 4B,C). Notably, HCC patients with more POSTN+ CAFs had shorter overall survival and progression-free survival (online supplemental figure 4D), indicating that POSTN+ CAFs may promote tumor progression and adversely impact patient prognosis.

High infiltration of POSTN+ CAFs is associated with T-cell exclusion

To gain further insight into the potential functional signals that induce POSTN+ CAFs, we conducted an analysis to identify genes that differentially expressed between POSTN+ CAFs and other CAF subgroups. Subsequently, we performed GSEA to examine the enriched pathways. online supplemental figure S5A illustrates the results of this analysis. The upregulated genes in POSTN+ CAFs were significantly enriched in several pathways associated with ECM remodeling in the TME, including the ECM-receptor interaction, the PI3K-AKT signaling pathway, focal adhesion, and the TGF-β signaling pathway (figure 4A). Numerous studies have highlighted the crucial role of ECM-related CAFs in promoting resistance to immunotherapy by facilitating T-cell exclusion.37–39 To validate this finding, we assessed the correlation between the infiltration of POSTN+ CAFs and T-cell exclusion in three independent cohorts downloaded from the Tumor Immune Dysfunction and Exclusion Database (http://tide.dfci.harvard.edu/) (figure 4B–D). Our findings revealed a significant positive correlation between the infiltration of POSTN+ CAFs and T cell exclusion. Furthermore, we observed a significant negative correlation between the infiltration of POSTN+ CAFs and CD8+ T cell infiltration (online supplemental figure S5B–C). To investigate the spatial organization of POSTN+ CAFs and CD8+ T cells, we performed spatial transcriptomic analysis. Interestingly, our observations revealed that POSTN+ CAFs were responsible for excluding T cells from the malignant area within HCC tumor tissue sections(figure 4E). Similar to the findings reported by Liu et al, POSTN+ CAFs formed a tumor immune barrier structure that excluded T cells in non-responders to ICB therapy but not in responders (figure 4F). Bayesian-enhanced spatial data highlighted that regions with high POSTN expression exhibited elevated levels of CD68 and COL1A1 in non-responders, while regions with high CD3D expression demonstrated lower levels of POSTN in responders (figure 4G). In addition, we explored the differences in the expression of immune checkpoint-associated genes between POSTN+ CAFs with high infiltration and POSTN+ CAFs with low infiltration in the TCGA-LIHC cohort. The results suggested that the presence of POSTN+ CAFs is associated with immunotherapy efficacy (figure 4H). Importantly, analysis of patients with progressive disease (PD) compared with those with a partial response (PR), based on the dataset from Chuan-Yuan Wei et al40, revealed significantly greater infiltration of POSTN+ CAFs in PD patients. This suggests that patients with a high proportion