Integration of single-cell and spatial transcriptomics reveals fibroblast subtypes in hepatocellular carcinoma: spatial distribution, differentiation trajectories, and therapeutic potential

Abstract

Background: Cancer-associated fibroblasts (CAFs) are key components of the hepatocellular carcinoma (HCC) tumor microenvironment (TME). regulating tumor proliferation, metastasis, therapy resistance, immune evasion via diverse mechanisms. A deeper understanding of the l diversity of CAFs is essential for predicting patient prognosis and guiding treatment strategies.

Methods: We examined the diversity of CAFs in HCC by integrating single-cell, bulk, and spatial transcriptome analyses.

Results: Using a training cohort of 88 HCC single-cell RNA sequencing (scRNA-seq) samples and a validation cohort of 94 samples, encompassing over 1.2 million cells, we classified three fibroblast subpopulations in HCC: HLA-DRB1 + CAF, MMP11 + CAF, and VEGFA + CAF based on highly expressed genes of which, which are primarily located in normal tissue, tumor boundaries, and tumor interiors, respectively. Cell trajectory analysis revealed that VEGFA + CAFs are at the terminal stage of differentiation, which, notably, is tumor-specific. VEGFA + CAFs were significantly associated with patient survival, and the hypoxic microenvironment was found to be a major factor inducing VEGFA + CAFs. Through cellular communication with capillary endothelial cells (CapECs), VEGFA + CAFs promoted intra-tumoral angiogenesis, facilitating tumor progression and metastasis. Additionally, a machine learning model developed using high-expression genes from VEGFA + CAFs demonstrated high accuracy in predicting prognosis and sorafenib response in HCC patients.

Conclusions: We characterized three fibroblast subpopulations in HCC and revealed their distinct spatial distributions within the tumor. VEGFA + CAFs, which was induced by hypoxic TME, were associated with poorer prognosis, as they promote tumor angiogenesis through cellular communication with CapECs. Our findings provide novel insights and pave the way for individualized therapy in HCC patients.

Keywords: Angiogen esis; Hepatocellular carcinoma; Hypoxic microenvironment; Tumor-associated fibroblasts; VEGFA + CAF.

Fig. 1

Isolation of fibroblasts in HCC. A UMAP shows cell clustering of cell types in the discovery cohort. B UMAP plot showing highly expressed genes in cell types. C Scatterplot showing the classification of fibroblasts and mural cells based on canonical gene markers. D Bar graph showing the percentage of fibroblasts/mural cells as well as undetermined cells in each mesenchymal subpopulation based on the classification results. E UMAP plot showing the results of fibroblasts and mural cells demarcation in the discovery cohort using the consensus gene signature

Fig. 2

Subpopulation and functional enrichment analysis of fibroblasts. A UMAP shows that fibroblasts can be clustered into three subpopulations. B UMAP plot showing the expression of three typical genes (MMP11/HLA-DRB1/VEGFA) in three fibroblast subpopulations. C Heatmap showing top 5 highly expressed genes in each subpopulation of fibroblasts. D–F Bar plot showing top 5 terms or pathways significantly enriched for HLA-DRB1 + CAF, MMP11 + CAF, VEGFA + CAF. G Heatmap showing the different subtypes of fibroblasts' progeny pathway scores. H Box line plots showing the scores of basement membrane, interstitial collagen, collagen, ECM glycoproteins, and proteoglycans in comparison for three subpopulations of fibroblasts. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001

Fig. 3

Spatial distribution of the three fibroblast subtypes. A Bar graphs show the percentage of the three fibroblasts among various sample types. B RO/E showing the abundance of the three fibroblasts to be enriched in AL and tumour tissues C Accuracy validation using a pseudobulk dataset generated from a single-cell transcriptome. Bar graphs show the correlation between CIBERSORTx estimates and true abundance for each cell type (Pearson r). D Box plot showing the comparison of the abundance of the three fibroblasts in normal and tumour tissues based on the deconvolution results of CIBERSORTx. E Spatial distribution of the three fibroblasts in the HCC1_L slide. F Spatial distribution of the three fibroblasts in the HCC2_L slide. G CellTrek deconvolution based on the HCC sections in the distribution of the three fibroblast subtypes in,HCC3_L,HCC4_L. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001

Fig. 4

Cell-state transition trajectory inference of CAFs. A UMAP plot demonstrating possible differentiation trajectories between the three fibroblasts. B UMAP plot demonstrating the distribution of fibroblasts in pseudotime. C Violin plot demonstrating the comparison of the three fibroblasts in pseudotime. D Distribution of the three fibroblasts as well as fibroblasts in normal and tumour tissues in pseudotime. E Heatmap demonstrating genes in fibroblasts during their differentiation over time and clustering analysis dividing them into three modules: progenitor, activation and differentiation. F Trends in gene expression of the three modules over time. G Trends in gene expression of the three modules over time in normal and tumour tissues. H Hypoxia, TGFb, VEGF, WNT pathway scores over time in normal and tumour tissues

Fig. 5

VEGFA + CAF was associated with poor prognosis. A UMAP plot demonstrating the distribution of Scisso + as well as Scissor − cells among the three fibroblasts. B Bar graph demonstrating the percentage of Scisso + cells in the fibroblast subpopulation. C Kaplan–Meier survival curve demonstrating the differences in OS as well as RFS in patients with different Scisso + as well as Scissor − cell abundance. D Demonstration of survival-independent cells, Scisso + as well as Scissor − cell subpopulations of highly and lowly expressed genes. E Forest plot demonstrating the relationship between the three fibroblast abundances and patient OS (in TCGA-LIHC). F Kaplan–Meier survival curve demonstrating the differences in OS in patients with different VEGFA + CAF abundances. G Kaplan–Meier survival curve demonstrating the differences in OS in patients with different VEGFA + CAF abundances in the two external validation sets, VEGFA + CAF was not only associated with poorer OS but also poorer RFS

Fig. 6

Cell communication between VEGFA + CAF and CapECs contribute to angiogenesis.in HCC. A UMAP plot demonstrating subpopulations of endothelial cells. B Heatmaps illustrating the top five most highly expressed genes in the various subpopulations of endothelial cells. C OR plots demonstrating the tissue distribution preferences of different fibroblasts as well as endothelial cell subpopulations. D VEGFA + CAF and CapEC scores, colocalized spots, and abundance distribution in the HCCT_5_3 sample. E VEGFA + CAF and CapEC scores, colocalized spots, and abundance distribution in the HCCT_5_2 sample. F Scatter plots demonstrating the activity of all pathways, as well as the VEGF pathway alone. G Heatmap demonstrating the strength of outgoing and incoming interactions of VEGF pathways between endothelial cells and fibroblast subpopulations

Fig. 7

Predicting prognosis and making therapeutic decision choices for HCC patients by machine learning using highly expressed genes in VEGFA + CAF. A Kaplan–Meier survival curve between VEGFA + CAF scores and patients' OS and RFS in four large cohort bulk seq datasets. B Using VEGFA + CAF highly expressed genes, 101 combinations of 10 machine learning predictions were used to predict the OS of patients with HCC, and the top 20 patients with the highest average C-index in the four datasets were selected for presentation. The top20 with the highest average C-index in the four datasets are shown. C C-index of the model with the best predictive performance in the four datasets. D Risk scores generated by machine learning are associated with poorer prognosis in HCC patients. E Demonstration of the prediction of response or non-response to sorafenib in patients with HCC using 7 machine learning using VEGFA + CAF highly expressed genes (in auc form)