The hepatic single-cell RNA sequencing (scRNA-seq) datasets of HCC patients and healthy controls, including GSE112271 (22), GSE149614 (23), GSE151530 (24) and GSE156625 (25), were obtained from the Gene Expression Omnibus (GEO) database. Quality control and pre-processing procedures were performed using Seurat (4.0.5, https://satijalab.org/eurat/) R toolkit (26). To avoid the influence of abnormal cells and technical noise on downstream analysis, low-quality cells such as doublets and empty droplets were removed. In addition, cells with mitochondrial gene expression >10%, or number of detected genes < 200 or > 5000 were also removed. Samples with > 1000 cells were retained for the analysis. Finally, the data of 165,932 cells, including 126,345 cells from 48 HCC patients and 39,587 cells from 14 healthy donors, were used for further analysis.

In order to minimize the technical batch effects among individuals and experiments, the “RunHarmony” function in R package harmony (27) was used to integrate all cells from HCC patients and healthy donors. R package “Seurat” was used for principal component analysis (PCA) and dimensionality reduction. The top 4000 variable genes were used for PCA to reduce dimensionality. The dimensionality of the scaled and integrated data matrix was further reduced to two-dimensional space based on the first 30 principal components (PCs), and visualized by t-Distributed Stochastic Neighbor Embedding (tSNE). The cell clusters were identified based on a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm with a resolution of 1. According to the expression levels of some well-known markers (24), the cells were annotated as B cells, T/NK cells, myeloid cells, fibroblasts, endothelial cells, and epithelial cells.

To identify the metabolic subpopulations of epithelial cells from HCC patients, the cell were first reintegrated using RunHarmony (27) and then classified into 42 clusters in an unsupervised manner using the first 30 PCs and the resolution of 2. The mean expression values of metabolic genes in these cells (28) were downloaded from Kyoto Encyclopedia of Genes and Genomes (KEGG), and was calculated for each cluster. Consensus clustering was then performed to determine the optimal number of stable metabolism-related epithelial subpopulations for HCC according to the expression matrix of metabolic genes for 42 clusters using the “ConsensusClusterPlus” package in R (1000 iterations, 80% resampling). he 42 clusters were divided into two metabolism-related subpopulations (MRS1 and MRS2), and the “FindAllMarkers” function was applied to identify the specific genes for each MRS, with log2 fold-change (avg_log2FC) > 0.25, detectable expression in at least 25% of the cells and percentage ratio > 1 as the criteria. The top 50 genes with higher avg_log2FC for each MRS were designated as the respective MRS1 and MRS2 genes.

The single-sample gene sets enrichment analysis (ssGSEA) scores of 85 KEGG metabolic pathways was calculated for each cluster based on the gene expression level in each cell (29). The differentially activated pathways between MRS1 and MRS2 clusters were then identified by Wilcoxon rank-sum test with p < 0.05 as the cut-off.

The myeloid cells were reintegrated using RunHarmony (27) and classified into 22 unsupervised clusters using the first 20 PCs and resolution of 0.8. The specific genes for each cluster were identified using FindAllMarkers with expression percentage ≥ 0.25 and avg_logFC ≥ 0.25 as the criteria. Based on established biomarkers, the 22 clusters were classified into eight major subpopulations.

The molecular interaction networks between the epithelial MRSs and myeloid subpopulations for HCC patient were identified using CellPhoneDB (30). The ligand-receptor pairs with p value < 0.05 were screened for the different cell clusters.

Single cell proliferation was estimated by predicting the cell cycle distribution using the “CellCycleScoring” function, which is based on the expression levels of ten genes that are upregulated in cycling cells (ASPM, CENPE. CENPF, DLGAP5, MKI67, NUSAP1, PCLAF, STMN1, TOP2A, TUBB) (31). For each of these genes, a background set of 100 genes with the smallest difference in average expression levels was selected. The average expression of the background gene set was then subtracted from each signature gene, and the average of the resulting values was calculated as the proliferation score. To identify metabolic subgroups of HCC patients at bulk levels, the MRS1 and MRS2 scores for HCC patients in TCGA-LIHC and GSE14520 cohort were calculated by GSEA in R package “GSVA” (32) based on the expression levels of the MRS1 and MRS2 genes at single-cell level. The samples were then assigned to the MRS1 and MRS2 subgroups on the basis of the respective scores.

The gene expression profiles of HCC patients were obtained from the GSE14520 (33) dataset of GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA-LIHC cohort of the UCSC Xena website (https://xenabrowser.net/datapages/). The robust multi-array average (RMA) algorithm in the affy package was used to pre-process the array profiles. After background correction, quantile normalization and probe summarization, the gene expression profile was generated based on the platform providing gene and probe mappings. Samples with overall survival (OS) above zero-days were selected for further analysis.

The tumor purity and stromal/immune infiltration were calculated for each sample based on bulk transcriptomic profile using the ESTIMATE algorithm (34). The tumor purity, ESTIMATE score, immune score and stromal score were calculated. In addition, the signature genes of 29 immune cell types and immune-related pathways were obtained from a previous study (35), and the abundance of these signatures in each patient was estimated by ssGSEA using the R package GSVA (32). Since the anti-cancer immune response determines the fate of tumor cells (36), the specific signatures associated with each immune response pathway were obtained from a previous study (37) and ssGSEA was performed.

All statistical analyses were performed using R software (version 4.0.4). The continuous variables between two groups were compared using the Wilcoxon rank-sum test, and Fisher’s exact test was used to compare the categorical variables. The prognostic factors were identified using the log-rank test. All tests were two-tailed and p-value ≤ 0.05 was considered statistically significant.

To explore the inter- and intra-tumoral heterogeneity of HCC at the single-cell level, we analyzed the scRNA-seq data from 48 HCC patients and 14 healthy controls. As shown in Figure 1A , a total of 165,932 individual cells from nine samples were clustered into six major cell types. The clusters were explored by PCA and visualized by t-distributed stochastic neighbor embedding (t-SNE). The cells in each cluster were annotated with the canonical markers ( Figures 1B, C ; Supplementary Figure 1A , and Supplementary Table 1 ) as B cells (6150), endothelial cells (14,926), epithelial cells (44,745), fibroblasts (4871), myeloid cells (22,834) and T/NK cells (72,406). We compared the enrichment of these cell types between HCC patients and healthy controls ( Figures 1D, E ), and found that epithelial cells were significantly enriched in HCC patients compared to the controls (p = 0.011; Figure 1F and Supplementary Figure 1B ), whereas the percentage of endothelial cells was significantly lower in HCC patients (p = 0.053; Supplementary Figure 1B ). In addition, the T/NK and B cells were also considerably reduced in HCC patients (T/NK cells: p = 0.062; B cells: p = 0.025; Figure 1E and Supplementary Figure 1B ), while the percentage of myeloid cells was similar in both groups (p = 0.17; Figure 1E and Supplementary Figure 1B ). These results indicated that epithelial cells and immune cells likely play important roles in the pathogenesis of HCC.

Recent studies have shown that metabolic reprogramming is a hallmark of cancer (7). To evaluate possible metabolic reprogramming in the epithelial cells of HCC tumors, we extracted 1679 metabolism-related genes from the KEGG database (28). Consensus clustering and PCA revealed a distinct metabolic gene expression pattern of the epithelial cells from HCC patients compared to the controls, whereas the epithelial cells from the controls samples showed high similarity ( Supplementary Figures 2A, B ). Furthermore, the global shifts in metabolic gene expression between and within epithelial cells from HCC patients and healthy controls were measured by the correlation distance. The distance between epithelial cells from HCC patients and healthy controls or within the cells from HCC patients was significantly greater than that within the cells from healthy controls ( Supplementary Figure 2C ), indicating considerable metabolic heterogeneity among the epithelial cells from HCC patients.

To confirm this metabolic heterogeneity, we re-clustered the epithelial cells from HCC patients into 42 clusters ( Supplementary Figures 3A, B ), which were then classified into two heterogeneous subpopulations (MRS1 and MRS2) based on metabolism-related genes expression matrix using consensus clustering ( Figure 2A and Supplementary Figure 3C ). We also identified the DEGs between MRS1 and MRS2 cells at the single-cell level ( Figure 2A and Supplementary Table 2 ), and their functional enrichment analyses revealed distinct metabolic patterns of the MRS1 and MRS2 cells. As shown in Figure 2B , pathways related to the metabolism of glycine, serine, threonine and other amino acids were obviously activated in the MRS1 cells, while glycolysis/gluconeogenesis was significantly activated in MRS2 cells. Furthermore, single-cell metabolic pathway analysis also revealed activation of tyrosine, glycine, serine, threonine, and phenylalanine metabolism, and phenylalanine, tyrosine and tryptophan biosynthesis in MRS1 cells. In contrast, glycolysis and gluconeogenesis were up-regulated in MRS2 cells at the single-cell level ( Figure 2C ). Consistent with the role of glycolysis in supporting tumor cell proliferation (39), the MRS2 cells also showed higher proliferation scores ( Figure 2D ), and more than 50% of these cells were either in S or G2M phase ( Figure 2E ). Taken together, the epithelial cells from HCC patients were classified into two distinct metabolic subtypes, of which the MRS2 cells had significant activation of glycolysis/gluconeogenesis and higher proliferation rates.

To further explore the relationship between the metabolism-related epithelial subpopulations and clinical outcomes, we calculated the MRS1 and MRS2 scores for HCC patients in TCGA-LIHC and GSE1450 cohorts. As shown in Figures 2F, G , higher proportion of MRS1 corresponded to favorable overall survival (OS) (log-rank test, p = 0.00018), whereas the predominance of MRS2 correlated to poor OS (log-rank test, p = 0.0019). We validated these results in an independent cohort, and the trend was consistent with the observations in TCGA-LIHC cohort ( Figures 2H, I ). Accordingly, the HCC patients were divided into the MRS1 and MRS2 groups according to the respective enrichment scores ( Figures 3A, B ). As expected, patients in the MRS2 group had significantly wore OS compared to the MRS1 patients (log-rank test, p < 0.001, Figure 3C ). Furthermore, multivariate COX regression analysis revealed that MRS2 was an independent predictor of worse prognosis in HCC after adjusting for stage, grade, HBV and HCV infection, hepatic cirrhosis and alcohol ( Figure 3D ). In another independent cohort as well, the MRS2 patients showed worse OS compared to the MRS1 patients ( Figures 3E, F ). Taken together, HCC can be divided into two metabolic phenotypes, and the predominance of epithelial cells with enhanced glycolysis/gluconeogenesis is linked to worse outcome.

Studies show that the immune cells in the TME can promote tumor development and progression by interacting with tumor cells (40, 41), leading to poor outcomes. Therefore, we next analyzed the differences in various immune-related factors between the two metabolic subtypes. Patients with the MRS2 phenotype exhibited higher ESTIMATE score, immune score and stromal score, which corresponded to lower tumor purity ( Figures 4A–D and Supplementary Figures 4A–D ). Furthermore, most immune cell populations, and the activity of immune-related pathways were significantly enriched in HCC patients with the MRS2 phenotype ( Figure 4E and Supplementary Figure 4E ). The regulatory T cells (Tregs) and myeloid-derived suppressor cells (MDSCs) that contribute to the immunosuppressive TME (42, 43) were particularly enriched in the MRS2 patients ( Figure 4F and Supplementary Figure 4F ). These results suggested that the worse prognosis associated with the MRS2 subtype can be attributed to greater infiltration of Tregs and MDSCs.

Given the positive correlation observed between the MRS2 epithelial cells and immune cells in HCC ( Supplementary Figures 5A–D ), we validated this result at the single-cell level. While MRS1 cells were not significantly correlated to the different immune cell populations, the MRS2 cells were co-enriched with myeloid cells ( Supplementary Figures 5E–G ). In addition, the percentage of MRS2 cells was positively correlated with that of myeloid cells (R² = 0.29, p = 0.00011; Supplementary Figure 5G ), indicating that these co-enriched myeloid cells may contribute to tumor progression. Thus, we dissected the myeloid subpopulations by integrating their single-cell transcriptomic data from HCC patients ( Supplementary Figure 6A ), and found that the myeloid cells were re-clustered into 22 clusters ( Supplementary Figure 6B ). According to the cluster-specific genes and established myeloid markers, the 22 clusters were divided into eight major subpopulations ( Figure 5A ; Supplementary Figures 6C, D and Supplementary Table 3 ).

We also explored the ligand-receptor interactions between the epithelial and myeloid subpopulations using CellPhoneDB (30), and detected significantly more interactions between the MRS2 cells and myeloid subpopulations compared to that among MRS1 cells and myeloid subsets ( Figure 5B ). The MRS2 cells interacted the most with the IL32+, FCN1+ and STMN1+ mono/macrophages ( Figure 5B ). Since these myeloid subpopulations contribute to the immunosuppressive TME (44, 45), we analyzed the type of interactions between MRS2 cells and these myeloid populations, and detected significant enrichment of CD47-SIRPA, TGFB1-TGFβ and other inhibitory ligand-receptor interactions ( Figures 5C–E ). The CD47-SIRPA interaction between tumor cells and myeloid cells is critical to phagocytosis blockage and immune escape (46), which further supports the immunosuppressive role of MRS2 cells in HCC.

We also explored the potential myeloid subpopulations co-enriched with MRS2 cells. The STMN1+ mono/macrophages were positively correlated with MRS2 cells at the single cell and bulk levels ( Figures 5F, G and Supplementary Figure 7A, B ). Furthermore, patients with high levels of STMN1+ mono/macrophages showed worse OS (TCGA-LIHC: log-rank test, p = 0.004; GSE14520: log-rank test, p = 0.011; Figures 5H, I ). Thus, the interaction between STMN1+ mono/macrophages and MRS2 cells via SPP1 and CD44 ( Figure 5E ) may result in a persistent immunosuppressive M2 state of myeloid cells (47), leading to the inhibition of immune surveillance. Taken together, the communication between MRS2 cells and myeloid cells shapes the immunosuppressive microenvironment of HCC, resulting in poor clinical outcomes in HCC patients with the MRS2 phenotype.