The HCC scRNA-seq data were obtained from Gene Expression Omnibus (GEO) database GSE156625 (Platform: GPL16791) (Sharma et al., 2020). A total of 43 tumor tissue samples and 14 adjacent control samples from 14 patients with HCC were included in this study. Healthy normal liver samples, fetal liver samples, and mouse liver samples were excluded. In addition, Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) data were obtained from the TCGA database (https://portal.gdc.cancer.gov/projects/TCGA-LIHC), including 374 HCC tissues with bulk RNA-seq data and the corresponding clinical information.

Single-cell data were processed in quality control, filtering the 1% of cells with the highest and lowest feature numbers of expression, as well as cells with more than 10% mitochondrial gene expression. After quality control, the data were integrated and analyzed based on SCTransform standardization (Hafemeister and Satija, 2019). The merged data were subjected to cell clustering analysis using the default parameters of the R package Seurat (Butler et al., 2018). A single-cell atlas was constructed as described previously (Liang et al., 2022b). Cell clusters were first identified based on the FindNeighbors and FindClusters functions of the Seurat package, and then downscaled and visualized as a single-cell atlas based on the Unified Modal Approximation and Projection (UMAP) algorithm (Becht et al., 2018). Cluster-specific marker genes were identified using the FindAllMarkers function of the Seurat package, and cell types were identified based on the expression of highly specific genes in the cell clusters as well as classical cell markers.

The copy number variations (CNVs) in tumor cells derived from the HCC tissues were calculated using the R package inferCNV (v1.6.0; inferCNV of the Trinity CTAT Project, https://github.com/broadinstitute/inferCNV) to identify subclones of diseased cells and to infer tumor evolution.

To explore biological status and functional differences between different cell subpopulations, we performed functional enrichment analysis based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using the R package clusterProfiler (Yu et al., 2012). The hypergeometric test was used to compare differences, and p < 0.05 was considered statistically significant.

Gene regulatory networks (GRNs) determine and maintain cell-type-specific transcriptional states; active transcription factors (TFs) and their target genes are commonly denoted as GRNs. To explore the key regulators that drive and maintain cellular state behavior, we constructed GRNs with TFs as their cores, divided them into several groups of co-expression modules, and inferred GRNs and cell states based on single-cell expression profiles with single-cell regulatory network inference and clustering (SCENIC) (Van de Sande et al., 2020). The binding patterns of the TFs were obtained from the JASPAR database (https://jaspar.genereg.net).

Development or external stimuli can change the functions of different cells in different ways, resulting in differential gene expression within the cells. Sorting these cells according to gene expression enabled further inference of cell development trajectory. In this study, the cell lineage trajectory of each cell type was inferred using Monocle3 (Trapnell et al., 2014), and the cells were projected into a low-dimensional space using UMAP.

Tumor cells can interact with surrounding cells through the circulatory and lymphatic systems, thereby influencing the development and progression of cancer. Thus, to further investigate cellular interactions in the HCC microenvironment, we used the R package iTALK for cellular communication analysis. iTALK identifies high-confidence ligand–receptor interactions between cells by identifying genes that are highly or differentially expressed in cell clusters and matching these genes to a built-in ligand–receptor database (Wang et al., 2019).

The R package limma (Ritchie et al., 2015) was used to normalize the TCGA-LIHC data. Receiver operating characteristic (ROC) curve analysis was performed using the R package pROC (Robin et al., 2011). In addition, to further explore the prognostic potential of highly expressed genes specific to each cell subpopulation, overall survival (OS) and recurrence-free survival (RFS) analyses of HCC were performed using the R package Survminer (https://rdocumentation.org/packages/survminer), and differences were calculated using the log-rank test with p < 0.05 considered significant.

The workflow of the study is shown in Figure 1A. To investigate the tumor ecological landscape in patients with primary HCC, we analyzed scRNA-seq data from GSE156625, which consisted of 43 tumor tissue samples and 14 adjacent control samples from 14 patients with HCC. After cell clustering and downscaling, 112,317 cells were grouped into 25 clusters manually identified according to 10 cell types with expressed markers consistent with previous laboratory and scRNA-seq studies (Sharma et al., 2020; Sun et al., 2021; Liang et al., 2022a) (Figure 1B), including HCC cells (AMBP, TTR, ARG1, GPC3, and ALB), hepatocyte (Hep; ORM1), epithelial cells (Ep; KRT19 and MUC1), endothelial cells (En; VWF, PECAM1, ENG, CD34, CDH5, KDR, and ICAM1), fibroblasts (Fib; ACTA2, COL1A1, and PDGFRB), naïve T cells (CD3D, CD3E, and CD3G), CD8⁺ T cells (IFNG, CD8A, and CD8B), natural-killer cells (NK; GNLY and KLRF1), B cells (MZB1, CD79A and MS4A1), and macrophages (Mac; HLA-DRA, HLA-DRB1, HLA-DRB5, CD68, and CD14) (Figure 1C). The annotation of HCC cells was validated by inferring CNVs based on the scRNA-seq data (Figure 1D). Subsequently, by comparing the differences in cell type composition between HCC and control samples, we found naïve T cells to be more abundant in HCC than in adjacent control tissues (Figure 1E). This result suggested significant lymphatic immune infiltration of tumor cells. Additionally, the abundance of fibroblasts was higher and that of endothelial cells was lower in HCC than in adjacent control tissues. In summary, we constructed a single-cell atlas of patients with HCC, revealing differences in cell-type composition in the HCC tumor microenvironment.

We performed a subpopulation analysis on the HCC cells and obtained seven subpopulations, including APOA1 ⁺ HCC, SPINK1 ⁺ HCC, APOA1 ⁺ SPINK1 ⁺ PLA2G2A ⁺ HCC, DDX5 ⁺ HCC, APOC2 ⁺ HCC, AIF1 ⁺ HCC, and CLDN10 ⁺ HCC (Figures 2A, B). Among these markers, APOA1 was not only highly expressed in the APOA1 ⁺ HCC subpopulation, but was also commonly expressed in other HCC subpopulations, and SPINK1 was also expressed in most HCC cells (Figure 2C). We conducted ROC curve analysis of APOA1 and SPINK1 based on the TCGA-LIHC data and found that high APOA1 and SPINK1 expression levels could serve as potential biomarkers for HCC, with AUC values of 0.782 and 0.690, respectively (Supplementary Figure S1). Enrichment analysis showed that SPINK1 ⁺ HCC and APOA1 ⁺ SPINK1 ⁺ PLA2G2A ⁺ HCC cells were significantly involved in bacterial infection-related pathways (Figure 2D). Notably, all HCC subpopulations except APOC2 ⁺ HCC were significantly involved in drug metabolism-related pathways, and most of these subpopulations were at the end of the cell developmental trajectory (Figure 2E). This result suggested that as HCC progresses, tumor cells may develop drug resistance by activating drug metabolism pathways, leading to poor outcomes in patients. Therefore, targeting these subpopulations could prevent the development of drug resistance in patients. We constructed GRNs to further explore the TFs regulating the subpopulations (Figures 2F, G). The results showed that subpopulations with drug-resistant potential were mainly regulated by TFs such as CRZB5, AR, HLF, and RARA, implying that tracking these genes is helpful for understanding HCC progression. In addition, we explored the prognostic potential of specific markers highly expressed in each subpopulation of HCC cells. The results suggested that high SPINK1 and DDX5 expression levels were significantly associated with lower HCC survival Figure 2H. In summary, drug metabolism pathways may be activated by tumor cells in the HCC microenvironment, resulting in increased difficulty in treatment.

CAFs are the most abundant cell type in the TME and are the center of cross-communication between various cells in the tumor mesenchyme. Thus, we performed a subpopulation analysis and identified seven major fibroblast subpopulations expressing different markers (Figures 3A–C). We found that the abundances of NDUFA4L2 ⁺ Fib and NDUFA4L2 ⁺ CRYAB ⁺ CAF cells were significantly increased in HCC, suggesting that these subpopulations may be involved in promoting tumor development (Figure 3D). Enrichment analysis revealed that all seven subpopulations were significantly associated with the focal adhesion pathway (Figure 3E), which is the point of contact between cells and the ECM and usually closely associated with cell migration. GRN analysis of the fibroblast subpopulations revealed that the target genes were divided into seven modules regulated by TFs including NR2F2, ETS1, TBX19, and HNF4A (Figures 3F, G). Ranking of the regulator specificity scores of the fibroblast subpopulations revealed that NR2F2 was the most specific regulator of both the CXCR4 ⁺ PMAIP1 ⁺ Fib and MT1M ⁺ Fib subpopulations and has the ability to promote tumor cell proliferation, epithelial-mesenchymal transition, and invasive features (Mauri et al., 2021) (Figure 3G). As a gene therapy tool for HCC, forced re-expression of HNF4A inhibits proliferation and eliminates cancer-specific features of target cells (Takashima et al., 2018). Therefore, the combined transduction of NR2F2 and HNF4A may be more stable for cellular reprogramming in HCC.

Pseudotime analysis showed that NDUFA4L2 ⁺ Fib cells were positioned earlier on the differentiation trajectory than were NDUFA4L2 ⁺ CRYAB ⁺CAF cells, presumably because NDUFA4L2 ⁺ Fib cells differentiated into NDUFA4L2 ⁺ CRYAB ⁺CAFcells over time, further promoting tumor development (Figure 3H). Subsequently, we evaluated the prognostic potential of specific markers highly expressed in each subpopulation of HCC and found that patients with higher CYRAB, FAP, and NDUFA4L2 expression had poorer OS (Figure 3I), whereas patients with lower APOC3 expression had poorer RFS (Figure 3J). Therefore, high NDUFA4L2 expression in fibroblasts may promote HCC development. To address the location of fibroblasts in the TME, we need to further investigate their connections with other cells.

The transformation of endothelial cells into mesenchymal fibroblasts after stimulation, known as the endothelial-mesenchymal transition (EMT), is an important source of fibroblasts. Our analysis confirmed the increased abundance of fibroblasts and decreased abundance of endothelial cells in the HCC group (Figure 1E). Thus, we performed a subpopulation analysis of the endothelial cells and obtained eight subpopulations (Figures 4A, B), named according to the markers they specifically expressed (Figure 4C). The abundances of STC1 ⁺ En and VWF ⁺ Encells were significantly increased in the HCC group (Figure 4D), and STC1 ⁺ En was the most abundant subtype. Most endothelial cell subpopulations were significantly associated with focal adhesion and cell adhesion molecule pathways (Figure 4E). GRN analysis of the endothelial cells revealed that the target genes were divided into three modules and regulated by TFs including KLF10, MEIS1, HOXB5, NR1I3, and TFE3 (Figure 4F). Among them, KLF10 was the most specific regulator of two subpopulations, FABP4 ⁺ En and SLC9A3R2 ⁺ En, and has important regulatory effects on perivascular fibrosis (Zhuang et al., 2022) (Figure 4G). Moreover, the SLC9A3R2 ⁺ En subpopulation, whose abundance was significantly reduced in HCC (Figure 4C), was not only significantly enriched in cell adhesion molecule pathways, but was also at an early stage of subpopulation cell differentiation, suggesting that it may gradually differentiate into other endothelial cell stages with the progression of HCC, in turn promoting tumor proliferation and migration (Figure 4H). Therefore, targeting KLF10 to regulate the abundance of the SLC9A3R2 ⁺ En subpopulation may help regulate HCC proliferation and migration. Survival analysis of marker genes showed that patients with low SLC9A3R2 and FCN3 expression had poorer OS (Figure 4I) and RFS (Figure 4J). Overall, these results suggest that SLC9A3R2 ⁺ En cells differentiate into other endothelial cell stages as HCC progresses and contribute to its further development.

TAMs are macrophages that infiltrate tumor tissues and exert tumor-promoting and immunosuppressive effects in various ways (Wium et al., 2018; Gadiyar et al., 2020; Hedrich et al., 2021). We identified nine macrophage subpopulations and subsequently assigned names based on genes specifically highly expressed in these subpopulations (Figures 5A–C). Among them, the APOC1 ⁺ Mac subpopulation was a major component of HCC macrophage ecology compared to levels in adjacent control tissue, whereas the abundance of the FCN1 ⁺ Mac subpopulation was significantly decreased in tumor cells. The APOC1 ⁺ SPP1 ⁺ TAM and HSPA1B ⁺ TAM subpopulations were particularly abundant in the HCC group and were therefore defined as TAMs (Figure 5D). Most macrophage subpopulations were significantly involved in pathways related to antigen processing and expression, suggesting an active immune function of macrophages in the HCC microenvironment (Figure 5E). GRN analysis of macrophages revealed that the target genes were divided into five modules regulated by TFs including NFIC, KLF15, RUNX3, HNF4G, HOXA13, and NFIB (Figures 5F, G). Pseudotime analysis showed similar differentiation trajectories for APOC1 ⁺ Mac and APOC1 ⁺ SPP1 ⁺ TAM cells, which were significantly more abundant in HCC, suggesting increased infiltration by macrophages highly expressing APOC1 with the development of HCC (Figure 5H). Macrophage subpopulations with high APOC1 expression were found at the end of the cell developmental trajectory and specifically expressed the M2 macrophage marker CD163 (Supplementary Figure S2), suggesting that these subpopulations may accelerate HCC development. Survival analysis of the subpopulation-specific markers showed that FCN1 and SPP1 had prognostic potential for HCC (Figures 5I, J). In HCC, the presence of TAM subpopulations correlates with poor outcomes (Graham and Pollard, 2022). In our study, TAM subpopulations were positioned at the end of the cell developmental trajectory and were regulated by the same TFs, implying that macrophage cell states could be changed to prevent the development of tumor cells.

T lymphocytes play an important effector-killing role in antitumor immunity. To explore the differences between lymphocytes in tumors and normal tissues, we identified eight CD8⁺ T cell subpopulations (Figures 6A, B) expressing different specific markers (Figure 6C). Among them, CD27 was specifically expressed in the BTG1 ⁺ RGS1 ⁺ Tcm and HSPA1B ⁺ Tcm subpopulations, thus classifying them as subpopulations of central memory T cells (Tcms) (Figure 6E). The subpopulations with double-positive expression of BTG1 and RGS1 were significantly enriched in tumor lesions (Figure 6D). In addition, the HSPA1B ⁺ Tcm subpopulation was significantly enriched in the HCC group, and HSPA1B was significantly enriched in all cell types in the HCC group (Figure 6C), suggesting its association with the remodeling of the tumor ecological niche in HCC. Enrichment analysis showed that most CD8⁺ T cell subpopulations were significantly involved in chemokine-and natural killer-mediated cytotoxicity-related pathways, especially HSPA1B ⁺ Tcm, which was significantly enriched in HCC. This result suggested that the subpopulation may be immunologically active in the HCC microenvironment, aggregating to lesions and exerting immunocidal effects (Figure 6F).

GRN analysis of CD8⁺ T cells showed that the target genes were divided into three modules regulated by TFs including THAP1, IKZF1, YBX1, FOXA2, NR2F1, and FOXP1 (Figures 6G, H). Pseudotime analysis showed that HSPA1B ⁺ Tcm was positioned at a late stage of development, suggesting that as HCC develops, this subpopulation of cells is stimulated by antigens, proliferates, and differentiates to exert antitumor effects (Figure 6I). Survival analysis revealed that among the subpopulation-specific markers, BTG1 and HSPA1B had prognostic potential for HCC (Figures 6J, K). In conclusion, HSPA1B ⁺ Tcm subpopulation plays an important role in the immune response to HCC, suggesting that increasing the abundance of this subpopulation may help inhibit tumor growth and metastasis.

To explore the relationship between non-tumor and tumor cells, we inferred their receptor–ligand interactions using the iTALK approach and constructed a cellular communication network divided into four modules: cytokines, growth factors, immune checkpoints, and other (Figure 7). We found that the CCL5 ligand secreted by BTG1 ⁺ RGS1 ⁺ Tcm and HSPA1B ⁺ Tcm cells bound to the SDC4/1 receptor on the surfaces of tumor cells (Figure 7B). Interestingly, the checkpoint CD24 ligand secreted by tumor cells bound to the SIGLEC10 receptor on macrophages to transmit inhibitory signals and reduce phagocytosis (Figure 7A). We also observed that tumor and endothelial cells released vascular endothelial growth factors, including VEGFC and VEGFA, as well as platelet-derived growth factors, including PDGFA and PDGFC (Figure 7C). Remarkably, in another module, the TIMP1 ligand secreted by CAFs bound to the receptor CD63 on the surfaces of tumor cells and TAMs, whereas the SPP1 ligand secreted by APOC1 ⁺ SPP1 ⁺ TAM cells bound to the ITGB1 receptor on CAFs (Figure 7D). These cellular communication networks reshape the tumor microenvironment and promote the development of HCC.