ScRNA-seq data were collected from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database and Mendeley Data (https://data.mendeley.com/). The inclusion criteria for scRNA-seq datasets were as follows: (1) liver tissue samples with clearly annotated disease status relevant to MASH, liver cirrhosis, or HCC; (2) the availability of raw or processed gene expression matrices and corresponding cell-type annotations; (3) the scRNA-seq analysis results have been confirmed by FISH, flow cytometry, immunostaining, or experimental studies in original article. We selected 5 cases of healthy controls and 5 cases of liver cirrhotic patients from GSE136103, (5) 3 cases of MASH patients from GSE159977. (6) Because the annual incidence of hepatocellular carcinoma arising from NASH-related cirrhosis was relatively low, we were unable to collect the dataset. Viral hepatitis-induced and NASH-induced HCC have been reported to share overlapping inflammatory and immune signaling pathways that drive malignant transformation (11). Therefore, we collected 6 cases of HBV-related cirrhosis with HCC from skrx2fz79n (12) as a surrogate. The data accession numbers and information were provided in Supplementary Table 1. Additionally, the RNA expression data and clinical information of HCC used in this study were obtained from The Cancer Genome Atlas public access web portal (TCGA-LIHC; https://portal.gdc.cancer.gov/).

Cases from GSE136103 and GSE159977 were obtained raw fastq files. They were processed with the Cell Ranger pipeline (version 6.0.2, 10X Genomics) and mapped to GRCh38 reference genome to generate count matrices. The count matrix of skrx2fz79n dataset was obtained from rds files (https://data.mendeley.com/datasets/skrx2fz79n/1). We used the “emptyDrops” function of the R package DropletUtils (version 1.14.2) (13) to remove barcode-swapped pseudo-cells. We applied “doubletFinder_v3” function of R package DoubletFinder (version 2.0.3) (14) to identify doublets. For quality control, cells with mitochondrial gene percentages less than 10%, genes expressed in over 3 cells, and detected genes between 200 to 4000 were kept. The R package Seurat (version 4.2.1) (15) was utilized to perform data normalization and dimensionality reduction. Gene expression counts were normalized and scaled by “SCTransform” function with glmGamPoi model, and we calculated a PCA matrix by “RunPCA” function. After PCA, we used the “RunHarmony” function with SCT assay for batch effect correction and datasets integration in the R package harmony (version 0.1.0) (16).

The top 40 harmony dimensions were used to carry out the uniform manifold approximation and projection (UMAP) dimensional reduction. We then constructed the nearest-neighbour graph by the “FindNeighbors” function with the reduction as ‘harmony’. The “FindClusters” function was then used to identify clusters with the resolution parameter of 0.09 of the whole object, which resulted in 15 clusters; the resolution of 0.3 for T cells and 0.4 for mononuclear phagocytes. To identify cluster-specific marker genes, we used “FindAllMarkers” function to select those detect in a minimum of 25% of cells within the cluster, displaying a p value < 0.05 in the Wilcoxon rank-sum test, and demonstrating a differential expression threshold of 0.25 log fold change (log2FC). The main cell types were annotated with known cell-type marker genes based on CellTypist (https://www.celltypist.org/) (17) and CellMarker 2.0 (http://117.50.127.228/CellMarker/) (18). We applied type-specific markers to annotate mononuclear phagocyte and T cell subtypes.

Trajectory and diffusion-pseudotime (DPT) analysis were generated using the R package destiny (version 3.22.0) (19). The number of nearest neighbors, k, was set to 10. The principal components were set to 30. The R package CellChat (version 1.6.1) (20) was utilized to explore the communication of macrophages and T cells into functionally relevant signaling pathways. We used the processed expression matrices from the “SCT” data slot of the Seurat object with corresponding annotations to create a CellChat object. CellChatDB.human was set as the ligand-receptor interaction database. We applied the minimum cell count criterion of 10.

We used the R package AUCell (21) to calculate the signature score of curated gene sets relate to macrophage and T cell functional states (22, 23) (Supplementary Tables 2, 3). AUCell could provide a relative measure of gene importance in each cell to evaluate the degree of gene set enrichment. The ranked gene expression matrix was built by the AUCell_buildRankings function, and then we calculated the AUC value using the AUCell_calcAUC function.

We performed GO enrichment analysis with the R package clusterProfiler (version 4.2.2) (24). The top 8 GO annotations were chosen for C1QA⁺ Mac and SPP1⁺ Mac for visualization. Enrichment scores for the eight selected GO annotations were calculated by a hypergeometric statistical test with a significance threshold of 0.01.

SCENIC (21) was a tool that utilizes scRNA-seq data to reconstruct gene regulatory networks. We used the pySCENIC (version 0.12.1) package in Python (version 3.10) to assess transcription factor enrichment and regulator activity of monocyte and nine macrophage subtypes. The grn function was used to infer gene co-expression relationships between transcription factors and their potential target. Then, ctx function was used to refine the regulons and separate direct and indirect target. Finally, the regulon activity was calculated by aucell function. The R package ComplexHeatmap (25) was used to visualize the regulon activity.

Based on cluster-specific marker genes, we identified the signature genes from FindAllMarkers function of Seurat (15). For macrophages, genes with adjusted p value < 0.01 by Wilcoxon rank-sum test, in a minimum of 50% of cells within the cluster and in a maximum of 50% of cells within the others cluster (pct.1 > 0.5 and pct.2 < 0.5), and demonstrating a differential expression threshold of 1 log fold change (avg_log2FC > 1) were defined as macrophage subtypes specific signature genes. For T cells, genes with adjusted p value < 0.01 by Wilcoxon rank-sum test, and demonstrating a differential expression threshold of 0.5 log fold change (avg_log2FC > 0.5) were defined as T cell subtypes specific signature genes. Then, we manually checked the cell subtypes specific signature genes. We removed ALB and ATP5F1E to confirm the signature genes were specific to cell subtypes. Furthermore, we performed SelectGene function of R package (26) (version 1.0) to calculate total entropy difference. The genes with greater total entropy difference tend to be more specific and would be retained as cell subtypes specific signature genes.

We collected 424 samples from the TCGA-LIHC cohort, including 50 normal liver samples and 374 HCC samples. Clinical information and gene expression matrix were retrieved using the R package TCGAbiolinks. The clinical information was retrieved using the GDCquery_clinic function. Bulk RNA-seq data were obtained using the GDCquery function with project set to “TCGA-LIHC”, data_category set to “Transcriptome Profiling”, data_type set to “Gene Expression Quantification”, and workflow_type set to “STAR–Counts. The TCGA-LIHC data was used to evaluate the prognostic performance of signature gene sets derived from macrophage and T cell subtypes clusters. We performed survival analysis using the Cox proportional hazards model implemented in the R package survival (version 3.3-1). And the ggsurvplot function was employed to correct patient age and plot Kaplan–Meier survival curves. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was applied to construct the prognostic model of the signature genes of SPP1⁺ Mac and Tc. Univariate Cox proportional hazards regression analysis was performed on the signature genes of SPP1⁺ Mac and Tc to screen genes significantly associated with OS in HCC (27).

We employed ImmuCellAI (28) to calculate immune cell enrichment score (IS) of each sample to correct the bias of immune cell abundance. Hepatocytes can impair the function of T cell (29), so we calculated the enrichment score of Tc (CA) to assign the deviation. Next, we used the expression of signature gene set as input of ssGSEA to obtain the enrichment score (ES) of SPP1⁺ Mac and Tc. Finally, the gene set enrichment score are as follows:

Gene set enrichment score of SPP1⁺ Mac (Mscore):

Gene set enrichment score of Tc (Tscore):

Experiments were conducted on 7–8 weeks old C57BL/6J male mice. Mice were purchased from the Hunan SJA Laboratory Animal Co., Ltd (Changsha, Hunan, China) and were bred in a specific pathogen‐free barrier facility at the Animal Center of Wuhan National Laboratory for Optoelectronics. 5 healthy mice as control. For Diet-induced MASH mice models, 5 mice were fed with MCD diet for up to 3 weeks. Liver fibrosis models were induced in 5 mice through an intraperitoneal injection of 25% CCL₄ in corn oil (Sinopharm Chemical Reagent Co., Ltd., China) twice a week for 6 weeks. 5 HCC mice models generated by hydrodynamic tail vein injection of plasmids carrying the Akt and N-ras genes were established over six weeks and purchased from Shouzheng Pharma (Wuhan) Biotechnology Co., Ltd. Finally, all mice were anesthetized and transcardially perfused with RNase-free PBS. Liver samples of mice were immediately collected and subjected to further experiments.

H&E staining and Masson’s trichrome staining were used to evaluate and examine the histopathologic changes in liver structure. The sections were imaged with PanoBrain (Tinyphoton, Hubei, China). For each specimen, 5 fields per tissue section were randomly chosen and quantified by ImageJ software (National Institutes of Health, USA).

For FISH, the livers were maintained and fixed with 4% paraformaldehyde (PFA) for 6 hours, then incubated with 30% sucrose with 4% PFA. The liver lobes were immersed in optimal cutting temperature (OCT) tissue blocks, and stored at -80 °C. Cryosections 30 μm thick were used for in situ amplification and hybridization. Padlock probes and primers were designed as target sequence of genes (C1qa, Cd3, Clec4f, Cxcl3, Cxcl10, Klf2, Trem2, Spp1, Akt and Nras). In brief, the sections were incubated with buffer (2×SSC, 35% formamide, 1% TritonX-100 in RNase-free PBS) contain probes overnight at 37 °C. Then, 5.7 U/μl T4 DNA ligase (TaKaRa, Japan) was added to ligate padlock probes at room temperature for 4 hours. After ligation, the sections were incubated with rolling circle amplification mixture (250 μM dNTP, 4 μM Dithiothreitol, 0.4 U/μl Recombinant RNase Inhibitor in ddH₂O) by Phi29 DNA polymerase (New England Biolabs, America) at 30°C for 2 hours. The sections were then incubated with monomer buffer (4% acrylamide, 0.2% bis-acrylamide, 0.1% ammonium persulfate, 0.2% tetramethylethylenediamine) at room temperature for 2 hours. Next, the sections were digested with Proteinase K (Biofrox, Germany) overnight in 37 °C. Finally, the sections were incubated with fluorescent oligo (Alexa Fluor 488, TAMRA, Alexa Fluor 647) complementary to DNA amplicon at 37 °C for 30 minutes. The sections were imaged on Olympus FLUOVIEW FV3000 confocal laser scanning microscopes for gene expression validations or cell subtypes analysis. The detailed probe information was listed on Supplementary Table 4.

The quantification of spatial distance was computed using cv2 python module. We used DAPI to stain the cell nucleus as the center point of the cell. We defined macrophages based on the fluorescent points of cell marker within the diameter of 20 μm, T cells were 10 μm. The FISH images were converted to the HSV color space, and which binary masks of valid regions were generated by applying specific thresholds. Next, the connected component analysis was performed on the masks to filter out small connected components and eliminate noise. The effective contours of fluorescent cells were extracted, and the distribution of their center points was calculated and visualized. For SPP1⁺ Mac and Tc, the shortest distance to a tumor cell was extracted from the pairwise cell-cell distance matrix of all cell distances, which defined as the distance to the closest tumor cell.

To uncovered the single-cell landscape during different processes of CLDs, we applied scRNA-seq data to characterize the dynamics from health liver, MASH, liver cirrhosis and HCC. In total, 99, 593 high-quality cells from 29 samples of 19 patients were analyzed after quality control in our study (Supplementary Table 1). After performing quality control, sctransform normalization, dimensionality reduction, batch effect removal, and clustering, the 13 distinct main cell types were annotated by CellTypist (17) and cell marker genes from CellMarker (30) (Figure 1a). We computed the spearman correlation of the main cell, and we observed a distinct correlation between adaptive immune cells (NKT and T cell) and innate immune cells (MP) and wanted to explore this further (Figure 1b). Each dataset with no obvious batch effects after batch correction (Supplementary Figure 1a). All cells with the number of genes per cell between 200 and 4000, and the proportion of mitochondrial gene counts less than 10% were selected (Supplementary Figure 1b). Of note, MP, T cell and NKT proportion were obviously dominant and distinct across different processes of CLDs, suggesting they might play a role in the immune microenvironment (Figure 1c). T cell expressed the T-cell receptor (TCR) signaling mediators CD3D (31), Mononuclear phagocyte (MP) was identified by the expression of C1QA, Natural killer T cell (NKT) was marked by XCL1 (32), Natural killer cell (NK) was identified by natural killer cell granule protein 7 (NKG7) (33), Endothelial cell was marked by known liver endothelial cell marker VWF (34), Cholangiocyte marker gene was KRT19 (35), B cell was identified by the canonical marker gene CD79A, Epithelial cell was defined by TNFAIP3 (36), Mesenchymal cell was marked by ACTA2 (37), Plasma cell was identified by the expression of JCHAIN (17), Plasmacytoid dendritic cell (pDC) was marked by IRF7 (38), Hepatocyte was defined by their classical marker ALB, Mast cell was positive for expression TPSB2 (39) (Figure 1d, Supplementary Figure 1c). To characterize the heterogeneous cellular compositions among the four processes of CLDs, we calculated the proportion of cell types by a bias-corrected and accelerated (BCa) bootstrap algorithm (40). Interestingly, we observed MPs and T cells (including T cells and NKT cells) proportion were distinct across different processes, which suggesting the progression heterogeneity of CLDs (Figure 1e). This indicated the complexity immune-specific microenvironment of CLDs (41). Hence, we proceeded to further annotated MPs and T cell subtypes with manual marker-based annotation method for in-depth analysis.

To pinpoint the transcriptional diversity in four processes of CLDs, we performed unsupervised clustering analyses of MPs. In addition to conventional dendritic cells (cDC1 and cDC2) and neutrophils, we annotated 10 other myeloid cell subtypes, including monocytes and 9 macrophage subtypes (Figure 2a). We found C1QA⁺ Mac highly enriched in healthy liver. SAMac and SPP1⁺ Mac was predominant in liver cirrhosis and HCC (Figure 2b). CXCL3⁺ Mac and CXCL10⁺ Mac highly expressed the genes of CXCL subfamily. CXCL3⁺ Mac exhibited the highest expression of CXCL8, CXCL3, CXCL2 and CCL20, which implied CXCL3⁺ Mac might be associated with liver disease progression and survival time (42, 43). CXCL10⁺ Mac was characterized by the inflammatory response chemokines of CXCL10 and CXCL9. It has been reported CXCL10 and CXCL9 were often localized with CXCL13-expressing T cells, which suggested they could participate in inflammation and antitumor reactivity (44). C1QA⁺ Mac highly expressed C1QA and C1QB, indicating that C1QA⁺ Mac might be involved in inhibiting tumor progression (45). SAMac specifically expressed TREM2, which was known to regulate scar-producing myofibroblasts (5). SPP1⁺ Mac was strikingly enriched in HCC and highest expressed immunosuppressive gene of SPP1 (46) (Figure 2c).

We next evaluated the macrophage polarization states of M1, M2, angiogenesis and phagocytosis properties to infer their functional states (Supplementary Figure 3a; Supplementary Table 2). CXCL10⁺ Mac showed the highest M1 score, whereas CXCL3⁺ Mac showed the higher M2 and angiogenesis score, suggesting CXCL10⁺ Mac take part in chronic inflammation, CXCL3⁺ Mac was an immunosuppressive phenotype (22). C1QA⁺ Mac display the highest M2 and phagocytosis score, revealing that C1QA⁺ Mac might promote tolerance and diminished pro-inflammation (47). Notably, C1QA⁺ Mac also expressed the highest level of CD5L and MARCO, reported to be associated with anti-inflammation (48, 49), which implicating that C1QA⁺ Mac could be a beneficial cell population of CLDs (Figures 2c, d). SPP1⁺ Mac showed the lowest M1 score, supporting SPP1⁺ Mac contributors to pro-tumor (50). Our observations elucidated the inflammatory response of CXCL3⁺ Mac, CXCL10⁺ Mac and C1QA⁺ Mac, and revealing SPP1⁺ Mac predominant present in HCC, implicating their important roles during CLDs processes.

We identified 14 T cell subtypes by unsupervised clustering analysis of T cells and NK cells (Supplementary Figure 2a), including naive T cell (Tn), NKT, exhausted T cell (Tex), effector T cell (Teff), proliferating T cell (Tproli), T helper cell (Th), T regulatory cell (Treg), stress response T cell (Tstr), effector memory T cell (Tem), cytotoxic T cell (Tc), Cycling NK & T, gamma delta T cell (γδT), mucosal-associated invariant T cell (MAIT) and tissue-resident memory T cell (Trm). But we didn’t observe changing patterns in T cell compositions during CLDs processes (Supplementary Figure 2b). Tn expressed high levels of naive markers. Tex demonstrated the highest expression of CXCL13. Tstr characterized by high expression of heat shock genes. Tc markedly expressed cytolytic activity-related genes. Tem showed high expression of immunoglobulin-related genes (Supplementary Figure 2c). To understand the state of T cells, we calculated the scores of naive, activation, cytotoxicity, and exhaustion onto the UMAP, which was following their expected functions (Supplementary Figure 2d; Supplementary Table 3).

The CLDs processes were dynamic (51). To reveal novel cellular and molecular mechanisms driving MASH to HCC, we further explored the differentiation trajectories and transcriptional regulation of macrophages and T cells subtypes. We applied Destiny (19) to infer the differentiation trajectory of macrophages and T cells subtypes. Macrophages subtypes displayed a trajectory that started with monocytes. FCN1⁺ Mac, S100A8⁺ Mac, and VSIR⁺ Mac were primary site in initial state same as monocyte. CXCL3⁺ Mac and CXCL10⁺ Mac mainly in intermediate state, which might drive the unique transcriptomic identities of MASH microenvironment (52). C1QA⁺ Mac, SAMac and SPP1⁺ Mac were dominating located in terminal state (Figures 3a, b), SAMac has been proved accumulation within the fibrotic niche (5). These findings supported that the HCC was correlated with C1QA⁺ Mac and SPP1⁺ Mac. T cells displayed two paths started with Tn. The path 1 ending in a terminally differentiation of Tex, Tc, γδT and Trm. The path 2 contained Tn, Treg and MAIT (Supplementary Figures 4b, c). These suggested the path 1 of T cells might be associated with HCC (23).

To investigate the transcriptional regulation dynamics of different macrophage subtypes during CLDs progression, we applied SCENIC (21) to calculate the activity of cell-type-specific regulators across CLDs processes. Interestingly, we found various transcription factors were highly expressed within the corresponding macrophage subtypes. The expression of ETS2 and IRF1 were consisted with transcription factors activities, which suggested that ETS2/IRF1 was specific targeting to CXCL3⁺ Mac/CXCL10⁺ Mac (Figure 3c; Supplementary Figure 4a). In addition, we explored gene ontology (GO) analyses to reveal CXCL3⁺ Mac was enriched in chemokine−mediated signaling pathway and inflammatory response (Figure 3d), CXCL10⁺ Mac was significantly enriched in response to biotic stimulus and immune response (Figure 3e). These implied that CXCL3⁺ Mac and CXCL10⁺ Mac might play important roles in MASH.

To validate the compositions of macrophage subtypes across CLDs progresses from scRNA-seq analysis, we constructed the mouse models of MASH, liver fibrosis and HCC (Figure 4a). Firstly, we evaluated the inflammation and fibrosis at the histological level by H&E and Masson’s trichrome staining to confirm the results in the scRNA-seq analysis (Figure 4b). Next, we utilized FISH to visualize those macrophage subtypes and directly calculated their cellular abundances. We used Clec4f as canonical liver macrophage marker. The well-conserved genes of C1qa, Cxcl3, Cxcl10, Trem2 and Spp1 were used to mark the signature genes (Figure 4c; Supplementary Table 4). If the signal of macrophage subtype signature genes were within 20 microns in diameter of Clec4f signal at 400 µm × 400 µm region, we recognized they were the corresponding macrophage subtype. As a result, C1QA⁺ Mac was enriched in healthy liver, CXCL3⁺ Mac and CXCL10⁺ Mac were commonly present in MASH, SAMac was enriched in liver fibrosis, and SPP1⁺ Mac was predominantly present in HCC (Figures 4d, e). Collectively, we affirmed the 5 macrophage subtypes presented at different processes of CLDs progresses.

We next aimed to understand the complex cellular interactions among diseases-associated macrophages and T cells. CellChat (20) was used to compare macrophages-T cells interactions in the previous versus subsequent processes of CLDs. The comparisons of macrophages-T cells interactions were distinct differences during the progression of CLDs. We observed CXCL3⁺ Mac and CXCL10⁺ Mac showed stronger crosstalks with Tex, Trm and γδT between health and MASH processes. The interaction strength of SAMac and SPP1⁺ Mac were increased in both MASH versus cirrhosis and cirrhosis versus HCC. C1QA⁺ Mac showed strong interaction strength in healthy liver and from MASH to cirrhosis (Figures 5a, b). MHC class II (MHC-II) were constitutively expressed on the surface of macrophages. MHC-II complexes depended on the distinct macrophage subtypes (53). We inferred cell-cell communication at the signaling pathway level. The interactions between CXCL3⁺ Mac and CXCL10⁺ Mac were mediated mainly by MHC-II signaling pathway. The interactions of MHC-II signaling pathway increased significantly from C1QA⁺ Mac, SAMac and SPP1⁺ Mac in healthy liver and HCC (Figures 5c-f). To further elucidated the molecular characteristics of diseases-associated macrophages and  T cells, we utilized entropy test to detect the signature gene sets from the differentially expressed genes of each cell subtype (Figures 5g, h). Notably, we found SPP1, IFI27, FOLR2 and SELENOP were the signature gene sets of SPP1⁺ Mac. These genes were reported to participate in immunosuppression and influence liver metabolic activites (12, 54, 55). The signature gene sets of Tc contained cytotoxicity-related molecules (FGFBP2, GNLY, GZMB, PRF1 and GZMH) and the transcription factor of suppress exhaustion (KLF2) (56, 57). These results indicated that SPP1⁺ Mac and Tc could be potential targets for HCC immunotherapy. We want to further identify the potential therapeutic targets from the signature gene sets. In brief, our results suggested a differentiated macrophages-T cells interactions and molecular phenotypes during the progression of CLDs.

ScRNA-seq provided an insight on the cellular heterogeneity in the TME. However, the scRNA-seq studies lacked survival data. We employed TCGA-LIHC data to verify the expression patterns of signature gene sets in HCC and reveal the key genes contributing to HCC. We developed a method of gene sets enrichment score to validate the signature gene sets from scRNA-seq (Figures 5g, h) based on TCGA-LIHC dataset. We calculated gene sets enrichment score of macrophage subtypes and defined as Mscore. The Tscore was defined for gene sets enrichment score of T cell subtypes. The Mscore and Tscore were extensions of the ImmuCellAI (28), designed to provide a simple and intuitive summary of predefined immune gene sets derived from scRNA-seq data, enabling the comparison of estimated cell populations in HCC. The details of the method calculating gene sets enrichment score were described in Methods. We found Mscore of SPP1⁺ Mac and Tscore of Tc were significant difference in the HCC group than the normal group (Figures 6a, b), which further validating SPP1⁺ Mac and Tc were associated with HCC. Next, basing on the clinical data from the TCGA-LIHC project, we confirmed that the higher Tscore of Tc was associated with a survival advantage in HCC, while an opposite effect was observed in SPP1⁺ Mac (Figure 6c), consistent with our hypothesis. Furthermore, we employed LASSO regression model and Cox regression analysis to reveal KLF2 was contributed to favorable prognosis and SPP1 were associated with poor outcomes in HCC (Figure 6d; Supplementary Figures 6a, b). As reported, KLF2 was favored effector differentiation and suppressed exhaustion (58). SPP1 overexpression was identified in tumor-associated macrophages across several cancer types and associated with poor prognosis (59).

In addition, we revealed the spatial distribution and cellular interaction of SPP1⁺ Mac and Tc. We quantified distance from SPP1⁺ Mac and Tc to its closest HCC cell by FISH micrographs. Figures 6e, f showed the FISH micrographs of the markers (top left, Supplementary Table 4), the corresponding dot plots of Zoom 3 representing the spatial distribution of SPP1⁺ Mac and Tc (bottom left), and the corresponding distance quantification from SPP1⁺ Mac and Tc to the closest HCC cell in Zoom 3 (bottom right). We found that Tc localized further away from the HCC cell (Figure 6e and Supplementary Figure 6c, with a median distance of 12.6 μm) compared to SPP1⁺ Mac that infiltrated and tightly surrounded the HCC cell (Figure 6f and Supplementary Figure 6c, with a median distance of 6.1 μm). The FISH micrographs also demonstrated that the number of SPP1⁺ Mac and Tc were significantly higher in HCC, and the double-positive cell pairs of Spp1 with Clec4f and Klf2 with Cd3 were close proximity to that in healthy liver (Supplementary Figures 6d-g, with a median distance of 0 μm). It suggested the potential crosstalk between SPP1⁺ Mac and Tc. Taken together, these findings suggested SPP1⁺ Mac contributed to the immunosuppressive microenvironment in HCC, Tc played anti-tumorigenic roles of HCC. We also determined that KLF2 and SPP1 involved in the progression of HCC and could be further clinical investigation.