Figure 1 shows the flowchart of this study. Six mice NASH datasets were retrieved and downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE83596 consists of 32 samples: 3 steatosis stage control, 3 NASH stage control, 4 fibrosis stage control, 4 tumor stage control, 3 steatosis, 3 NASH, 4 fibrosis, 4 non-tumor, 4 tumor. GSE189066 consists of 6 samples: 3 control, 3 NASH. GSE233767 consists of 8 samples: 4 control, 4 NASH. GSE207281 contains 3 control, 3 NASH. GSE205846 consists of 8 samples: 4 control, 4 NASH. GSE242881 contains 3 control, 3 NASH. GSE246221 contains 5 control, 14 NASH-associated HCC. Human NASH dataset GSE135251 was also downloaded from GEO database and it contains 10 control, 11 NAFLD patients with NAS score 1, 21 NAFLD patients with NAS score 2, 26 NAFLD patients with NAS score 3, 38 NAFLD patients with NAS score 4, 47 NAFLD patients with NAS score 5, 37 NAFLD patients with NAS score 6, 18 NAFLD patients with NAS score 7, 8 NAFLD patients with NAS score 8. All bulk-seq data were processed using R language and the fold change between control and NASH were calculated. Identification of differentially expressed genes was performed using the R package, “DESeq2”, and a p value <0.05 was used to identify the differentially expressed genes (DEGs).

The analysis of single-cell sequencing data utilized publicly available samples from GEO (GSE216836), including two controls and six NASH samples, and was processed in R. To identify cellular subpopulations in single-cell RNA sequencing (scRNA-seq) data, we used the Seurat package (version 4.3.0.1) for data processing. The workflow involved the following steps:1) Data Conversion and Quality Control: We converted 10x Genomics scRNA-seq data into Seurat objects using the Seurat package. Quality control was performed before specific analyses to ensure data integrity. We retained genes that were expressed in at least three cells and kept cells with at least 200 detected genes. Using the PercentageFeatureSet function, we calculated the percentage of mitochondrial content, ensuring that it was below 25%. 2) Normalization and Variable Gene Selection: We used the NormalizeData function to normalize gene expression, reducing technical noise and enhancing downstream analysis accuracy. We then identified the top 2000 highly variable genes using the FindVariableFeatures function. 3) Principal Component Analysis (PCA): We performed PCA to reduce data dimensionality and capture dominant data signals for downstream analysis. 4) Batch Effect Correction and Data Integration: To correct batch effects, we used the RunHarmony function from the Harmony R package, integrating data based on PCA results. 5) Further Dimensionality Reduction and Clustering: Using the Harmony-corrected dimensionality reduction results, we performed cluster analysis and visualized the data. UMAP was applied to further reduce dimensionality and visualize the integrated data, revealing distinct cell clusters. 6) Cell Type Annotation: Finally, cell clusters were annotated using the CellMarker database (http://xteam.xbio.top/CellMarker/) to assign cell types based on known marker genes.

Robust Rank Aggregation (RRA) was used to further identify robust DEGs from different datasets in an unbiased manner using a comprehensive ranking list algorithm, and a p value of <0.01 and log| FC| > 1 were considered to indicate significance (25).

Gene Set Enrichment Analysis (GSEA) and Gene Ontology (GO) pathway enrichment analysis were conducted using the “ClusterProfiler” R package (26).

“ConsensusClusterPlus” R package was applied to calculate how frequently HCC samples were grouped by RRA-selected hub genes. We used the proportion of ambiguously clustered pairs (PAC) to accurately estimate the optimal cluster number (K). Two clusters were identified, and further survival analysis was conducted by the Kaplan-Meier curve with the log-rank test. Principal component analysis (PCA) was performed by the “ggplot2” R package.

Gene mutation was conducted through “maftools” package, based on the somatic mutation data from TCGA-LIHC. And then we calculated tumor mutation burden (TMB) of each patient and compared TMB between the high- and low-risk groups. Survival analysis was also performed according to TMB score.

Hydrodynamic tail vein injection was performed as described in the literature (27). Plasmids were mixed in 2 mL of Normal Saline(NS) and injected into mice via the tail vein within 7 seconds. 12 mice were acclimated and randomly divided into a control group (n=6) and a NASH group (n=6). The model group received a rapid tail vein injection of a 2 mL solution containing 15 μg of ΔN90-beta-catenin, 15 μg of myr-AKT, and 5 μg of pCMV/SB plasmid mix, while the control group was administered an equivalent volume of saline. Mice were anesthetized and euthanized at 4 days, 4 weeks, and 2 months post-injection.

Liver tissue were collected, and protein was extracted using a RIPA lysis buffer, and protein concentrations were measured using BCA (Beyotime, China). The proteins were loaded onto 10–12% SDS-PAGE and transferred onto PVDF membranes, which were blocked using 5% non-fat milk at room temperature for 1 h, and incubated with primary antibodies overnight at 4°C. After washing the membrane with TBST three times, the membrane was incubated with the secondary antibody at room temperature for 1 h. Finally, a chemiluminescence reagent imaging system was used to detect the bands by using the Tanon 4800 system.

Liver tissues from wild-type mice and NASH models were thoroughly homogenized, and RNA was extracted using the TRIzol method. Reverse transcription and PCR were performed using Vazyme’s reverse transcription kit (catalog number R323) and PCR kit (catalog number Q341), respectively. The reverse transcription took place on the GeneAmp PCR System 9700 from Applied Biosystems, and PCR amplification was executed on the LightCycler 480 II system from Roche. All primers were sourced from Tsingke Biotechnology. The primer sequences for RT-PCR are as follows: Rplp0 Forward: GAAACTGCTGCCTCACATCCG, Reverse: GCTGGCACAGTGACCTCACACG; Trem2 Forward: CAGCACCTCCAGGAATCAAGA, Reverse: AGGATCTGAAGTTGGTGCCC; Gdf15 Forward: CTGGCAATGCCTGAACAACG, Reverse: GGTCGGGACTTGGTTCTGAG; Anxa2 Forward: GTGCCTACGGGTCAGTCAAA, Reverse: CACATTGCTGCGGTTTGTCA; Ttc39a Forward: CAGAAGGGCCACAAGGACTC, Reverse: AATCCTGGTGGGAAGCATGG.

All experiments were performed in triplicate. Melting curve analysis confirmed PCR specificity with single peaks. Ct values were analyzed using the 2^(-ΔΔCt) method, with Rplp0 as the reference, to calculate relative RNA expression levels.

In this study, we utilized the R package “oncoPredict” to assess drug sensitivity in LIHC based on gene expression data from patients. This tool enabled us to evaluate the potential effectiveness of various chemotherapeutic agents by calculating drug sensitivities. Additionally, it allowed for the identification of potential biomarkers that could predict the clinical response to specific treatments. We further explored the associations between these predicted drug responses and patient clinical characteristics to enhance the understanding of therapeutic outcomes in HCC.

The Tumor Immune Dysfunction and Exclusion (TIDE) method evaluates the responsiveness of cancer to immunotherapy using pre-treatment tumor expression profiles. It assesses various transcriptomic biomarkers, including TIDE, Dysfunction, Exclusion, MSI.score, TMB, CD274 and CD8 to predict patient responses to treatments. This approach helps in comparing the effectiveness of different biomarkers based on their predictive power regarding treatment outcomes, aiding in the optimization of immunotherapy strategies.

We conducted a RRA analysis on liver gene expression from six murine models of NASH, identifying 62 upregulated and 24 downregulated genes ( Figure 2A ). The distribution of these genes is illustrated in the volcano plots across the six analyzed GEO datasets ( Supplementary Figure 1 ). Gene Ontology (GO) enrichment analysis revealed significant enrichment of these genes in terms related to lipid metabolism and inflammatory processes, including lipid localization, fatty acid metabolic process, collagen-containing extracellular matrix, regulation of lipid metabolic process, triglyceride metabolic process, lipid transporter activity, cytokine activity, steroid metabolic process, lipid transport, regulation of inflammatory response, epoxygenase P450 pathway, and negative regulation of immune system processes ( Figure 2B ). The enriched pathways reveal significant involvement of metabolic processes that are integral to the development and progression of liver diseases leading to HCC. Critical pathways such as lipid localization, lipid transport, and the metabolic processes of fatty acids and triglycerides are highlighted. Dysregulation in these pathways can lead to excessive lipid accumulation, contributing to steatosis, a hallmark of NASH. Furthermore, abnormalities in lipid metabolism are directly linked to oxidative stress and lipotoxicity, which can exacerbate liver damage and promote fibrogenesis. Chronic inflammation is a key driver in the transition from NASH to HCC, and the pathways identified underscore its pivotal role. The activity of cytokines as mediators of inflammation is crucial. Abnormal cytokine profiles can lead to an unresolved inflammatory state, driving the progression of liver diseases. The genes identified through Robust Rank Aggregation (RRA) may influence the progression of disease by impacting these signaling pathways. This suggests a potential mechanism by which these genes could modulate key metabolic and inflammatory processes, thereby affecting the development and progression of liver diseases such as NASH leading to HCC. Understanding these interactions offers insights into how specific gene alterations contribute to disease dynamics and opens avenues for targeted therapeutic interventions.

Since NAFLD eventually develops into liver cancer, we therefore examined the role of the hub genes in liver cancer. The consensus clustering analysis ( Figure 3A ) demonstrates that the division into two clusters is the most stable and distinct configuration, as reflected by the sharp increase in CDF values up to k=2. This finding is visually supported by the consensus matrix ( Figure 3B ), which shows a clear segregation between the two clusters. The heatmap provides a detailed view of the differences in central gene expression between clusters, with genes in cluster 2 having higher hub genes expression levels ( Figure 3C ). The PCA plot further confirms the separation between clusters, emphasizing substantial molecular differences that are potentially clinically relevant ( Figure 3D ). Kaplan-Meier survival analysis ( Figure 3E ) showed that liver cancer patients of cluster 2 who highly expressed hub genes had a lower survival rate. Finally, the analysis of immune cell infiltration reveals significant variations in the immune landscapes of the two clusters, which indicates differences in tumor microenvironment and response to immunotherapies ( Figure 3F ). Cluster 2 shows a significant reduction in naive B cells, gamma delta T cells, resting NK cells, activated K cells, monocytes, M2 macrophages, and resting mast cells compared to Cluster 1 (Blue Bars). On the other hand, this cluster exhibits an increase in regulatory T cells (Tregs), M0 macrophages, and resting dendritic cells. The increase in Tregs and M0 macrophages, which are often associated with immunosuppressive activities, suggests an environment that may promote tumor growth and inhibit effective anti-tumor immune responses.

To further explore key genes involved in the progression of NASH, we conducted an in-depth analysis using liver sequencing data from patients with NAFLD. We aimed to identify key genes involved in the progression from NASH to HCC based on the following criteria: 1. These genes are positively correlated with the NAS score. 2. Compared to adjacent non-tumorous tissues, these genes are upregulated in cancerous tissues. 3. Patients with high expression of these genes have poorer survival outcomes than those with low expression. 4. These genes are upregulated in NASH-HCC mouse models. GSE135251 dataset includes liver biopsy NAS score (NAFLD Activity Score) ranging from 0 to 8, with 0 referring normal liver histology. We examined the correlation between all genes and NAS scores, especially those obtained by RRA ( Figure 4A ). We identified 11 genes with correlations greater than 0.4 with NAS scores, namely GPNMB, COL1A2, EPHB2, CDKNN1A, LGALS3, SDCBP2, COL3A1, CIDEC, FABP4, ANXA2, and TREM2 ( Figures 4B-L ); 7 genes with correlations between 0.3 and 0.4, namely CPXM1, LCN2, ADRA2A, SERPINE2, GPRC5B, UBD, and COL1A1 ( Supplementary Figures 2J–P ); 6 genes with correlations between 0.2 and 0.3, namely OSBPL3, ITGAX, GDF15, CX3CR1, TTC39A, and SLC15F2 ( Supplementary Figures 2D–I ); and 3 genes with correlations between -0.2 to -0.4, namely C8B, HAO2, and AQP7 ( Supplementary Figures 2A–C ).

Forest plots from TCGA hepatocellular carcinoma data show the prognosis of NAS score-related genes whose expression correlates with poor prognosis ( Supplementary Table S1 ). Based on the previous criteria for screening key genes for the transition of Non-alcoholic steatohepatitis to HCC, we further screened the 27 genes found to have a positive correlation with the NAS score. We first examined the expression of each of these 27 genes in NASH-induced HCC, and the volcano plot showed that 22 of these genes had elevated expression, which were Cdkn1a, Anxa2, Gpnmb, Cidec, Gdf15, Mmp12, Ttc39a, Itgax, Fabp4, Sdcbp2, Adra2a, Slc35f2, Col1a1, Lgals3, Cx3cr1, Ubd, Gprc5b, Col3a1, Col1a2, Osbpl3, Trem2 and Cpxm1 ( Figure 5A ).We next screened these 22 genes further and found that 13 of them were upregulated in cancer tissues compared to the paraneoplastic ones, which were GDF15, TTC39A, TREM2, ANXA2, UBD, FABP4, SERPINE2, COL1A1, COL1A2, OSBPL3, ITGAX, MMP12 and GPNMB ( Figures 5C–F , 6 , Supplementary Figure 3 ). After further screening of these 13 genes, we found that only 4 genes met all criteria for key genes in the transition from NASH to HCC, and they were GDF15, TTC39A, TREM2, and ANXA2. These 4 genes were highly expressed in patients with lower survival rates ( Figures 5G-J ). Moreover, we conducted a detailed examination of gene expression in liver tissue from mice with NASH-associated HCC. This model was induced by streptozotocin (STZ) combined with a high-fat diet (HFD), a well-established method for simulating the progression from NASH to HCC, which mimics the human disease pathology closely. Our analysis revealed a significant upregulation in the expression levels of four critical genes: Trem2, Anxa2, Ttc39a, and Gdf15 ( Figure 5B ). We also analyzed the expression levels of Trem2, Anxa2, Ttc39a, and Gdf15 across various stages of the disease, including steatosis ( Supplementary Figure 4A ), NASH ( Supplementary Figure 4B ), and fibrosis ( Supplementary Figure 4C ). Remarkably, we observed that these genes were not only upregulated in the NASH-associated HCC stage but also showed increased expression at earlier disease stages. The consistent upregulation of Trem2, Anxa2, Ttc39a, and Gdf15 across all stages of liver disease—from steatosis and NASH to fibrosis—suggests that these genes play a pivotal role in the progression of NAFLD.

To better simulate the transition from NASH to HCC, we administered oncogenic plasmids AKT and N90 and the sleeping beauty transposon (SB) via tail vein injection. Four days post-injection, fat accumulation was observed in the mouse liver ( Figure 5K ); by four weeks, severe liver damage occurred, and ACC1 and FASN protein levels significantly increased ( Figures 5L, M ). This confirms the successful creation of the NASH mouse model. Two months after injecting the AKT and N90 plasmids, there were significant morphological changes and severe damage in the liver, with evident fatty liver ( Figure 5N ). The mRNA levels of the four genes, GDF15, TTC39A, TREM2, and ANXA2, were significantly elevated, as shown in Figure 5O , indicating that these genes may play a crucial role in the transition from NASH to HCC.

Previous results showed significant differences in immune infiltration among liver cancer patients grouped by hub genes ( Figure 3F ). We further analyzed the role of key genes in immune infiltration in NASH-associated liver cancer patients using GSEA on the GSE164760 dataset. Specifically, there is a downregulation of Th1 and Th2 cell differentiation pathways in cancerous tissues, which may weaken immune surveillance of tumor cells, as Th1 cells promote the attack and elimination of tumor cells ( Supplementary Figure 5A ). Concurrently, key pathways such as Endocytosis, Necroptosis, Chemical Carcinogenesis – Reactive Oxygen Species, Cell Cycle, and Non-Alcoholic Fatty Liver Disease are upregulated, reflecting their pivotal roles in promoting disease progression and tumor development ( Supplementary Figures 5B–F ). Violin plots using the ESTIMATE algorithm confirmed higher Tumor Microenvironment (TME) scores in non-cancerous than in cancerous tissues of NASH-associated HCC patients ( Supplementary Figure 5G ), indicating more pronounced stromal and immune components in the former. Additionally, scatter plots ( Supplementary Figure 6 ) revealed a positive correlation between TREM2 and ANXA2 with both the StromalScore and ImmuneScore, suggesting their roles in modulating the tumor microenvironment. ESTIMATEScore, ImmuneScore, and StromalScore were higher in patients with high expression of TREM2 and ANXA2, whereas these scores did not differ in patients with high expression of GDF15 and TTC39A ( Supplementary Figure 7 ).

To further explore the possible role of GDF15, TTC39A, TREM2, and ANXA2 in the progression of NAFLD to HCC, we used the single-cell dataset GSE216836 to examine the expression of these genes. The UMAP plot effectively segregates various cell types, providing a comprehensive view of the cellular heterogeneity ( Figures 7A, B ). We found that the proportions of various cell types in control and NASH liver tissue samples. Notably, the proportion of macrophages significantly increases in the NASH samples compared to the control, while T cells and B cells show a relative decrease in their proportions( Figure 7C ). We annotated each cell cluster identified in our study using the CellMarker database, a comprehensive resource for cell markers across various tissues and organisms( Supplementary Figure 8 ).

Overall, the feature plots and dot plots indicate an increased expression of the four genes Trem2, Anxa2, Ttc39a, and Gdf15 in NASH ( Figures 7D–G , 8A-D ). Specifically, there is a noticeable increase in the expression levels of these genes in macrophages within NASH ( Figures 8E-L ). The analysis of cellular communication within the control and NASH conditions revealed distinct patterns of intercellular interactions. In the control condition, the network diagrams ( Figures 9A, B ) displayed a balanced communication pattern among various cell types, with endothelial cells, hepatocytes, and macrophages showing prominent interactions. In contrast, the NASH condition ( Figures 9C, D ) exhibited a significant increase in communication activities, especially involving inflammatory cells such as dendritic cells and macrophages, indicating heightened immune response and cellular stress. The Sema6 signaling pathway analysis ( Figure 9E ) identified key interactions, particularly between Trem2 and Plxna1, suggesting their pivotal roles in cellular navigation and immune regulation. Similarly, the GDF signaling pathway ( Figure 9F ) analysis verified the interaction between Gdf15 and Tgfbr2, crucial for tissue repair and fibrosis, which was often exacerbated in NASH. The complex connections of SEMA6-Trem2 in endothelial cells and macrophage networks ( Figure 9G ) and the distribution of the GDF pathway in hepatic stellate cells ( Figure 9H ) highlight the complex network of signaling events that coordinate cellular responses in liver disease. This analysis emphasizes the intricate interplay between gene expression and the immune landscape, providing insights into potential therapeutic targets and prognostic indicators in NASH-associated HCC.

Clinical prognostic analysis of the core genes TREM2, TTC39A, GDF15, and ANXA2 demonstrated that their expression levels were significantly correlated with the survival outcomes of patients diagnosed with HCC. Nomogram of the core genes showed varying expression levels across different patient samples, with higher expression of these genes correlating with poorer predicted survival probabilities ( Figure 10A ). Figure 10B illustrates the survival status, risk score distribution, and expression levels of the core genes. The mortality rate in the high-risk group is significantly higher than in the low-risk group, and the expression levels of the core genes are elevated in the high-risk group compared to the low-risk group. Decision curve analysis revealed a significant net benefit across various practical threshold probabilities for clinical decision-making, highlighting the robust potential of ANXA2, TREM2, GDF15, and TTC39A as biomarkers for predicting overall survival in patients with HCC ( Figure 10C ). The year-specific DCA curves further validated the efficacy of these biomarkers in forecasting short-term and long-term survival, essential for personalized treatment planning ( Figure 10D ). Patients were stratified into high-risk and low-risk groups based on gene expression profiles. High-risk patients, identified by the higher expression of the gene signature, demonstrated a more extensive and varied mutation pattern compared to the low-risk group ( Figure 10E ). The mutations of key genes in 371 LIHC samples were visualized using oncoplot, and the results showed that the high-risk group exhibited higher TP53 mutations ( Figure 10E ). This stratification validated the potential link between mutation burden and poorer clinical outcomes, suggesting that patients in the high-risk group were prone to experience more aggressive disease progression. Overall, these results suggest that the signature composed of these four core genes may serve as an important prognostic indicator for NASH-associated HCC.

We used GDSC database to predict the chemotherapy response of the common chemotherapy agents between the two groups. The drug sensitivity analysis demonstrates that differential drug responsiveness is associated with risk classification based on these four core genes ( Figure 11A ). We found that high-risk patients were insensitive to the following compounds: 615590, 667880, AZD1208, AZD5991, BEN, CHIR-99021, Dihydrorotenone, GSK2110183B, GSK2256098C, GSK2830371A. GSK626616AC, IAP_5620, JAK1_8709, LCL161, LMB_AB2, LY2109761, N-acetyl cysteine, OF-1, TAF1_5496, VTP-B ( Figure 11A ). Additionally, we assessed the sensitivity to some common anticancer drugs. Generally speaking, HCC patients with NASH exhibit higher sensitivity to the following drugs, as their IC50 values are close to 0: Camptothecin, Epirubicin, MG-132, Mitoxantrone, and Mycophenolic acid. However, for most other drugs, there are higher IC50 values, such as: Sorafenib, 5-Fluorouracil, Afuresertib, AGK2, alpha-lipoic acid, ascorbate (vitamin C), Cisplatin, Cyclophosphamide, glutathione, MIRA-1, N-acetyl cysteine, PRIMA-1MET, and Refametinib ( Figures 11B–E ). Moreover, we further analysis the drug sensitivity individually based on the four key genes expression level ( Figures 11B–E ). The groups were categorized based on the median expression level of the genes. Patients with gene expression levels above the median are classified into the high group, while those with expression levels below the median are categorized into the low group. Patients with high expression of ANXA2 are insensitive to 5-Fluorouracil, Cyclophosphamide, Epirubicin, MG-132, and Refametinib but are sensitive to Sorafenib, Afuresertib, AGK2, Ascorbate (Vitamin C), and N-acetyl cysteine (NAC) ( Figure 11B ). Patients with high expression of GDF15 show insensitivity to Afuresertib, whereas they are sensitive to 5-Fluorouracil, Glutathione, MIRA-1, NAC, and Refametinib ( Figure 11C ). Patients with high expression of TREM2 are sensitive to Sorafenib, 5-Fluorouracil, Afuresertib, AGK2, Alpha-lipoic acid, Vitamin C, Camptothecin, Cisplatin, Cyclophosphamide, Epirubicin, MG-132, PRIMA-1MET and Refametinib ( Figure 11D ). Patients with high expression of TTC39A are insensitive to Mitoxantrone and NAC, but are sensitive to Sorafenib 5-Fluorouracil, AGK2, Vitamin C, Cyclophosphamide and Epirubicin ( Figure 11E ). This finding suggests that these genes not only serve as prognostic biomarkers but may also guide therapeutic decision-making, particularly in selecting more effective personalized treatment regimens for HCC patients.

We then analyzed the relationship between the expression of key genes and various immune metrics in hepatocellular carcinoma associated with NASH. Patients with high ANXA2 expression exhibit higher TIDE and Exclusion scores ( Figures 12A, C ). Patients with high TREM2 and TTC39A expression show higher Dysfunction scores ( Figure 12B ). Additionally, patients with high TREM2 expression have lower MSI scores ( Figure 12D ), suggesting these genes may contribute to immune evasion mechanisms in the tumor environment. Notably, patients who respond to immunotherapy have lower ANXA2 expression levels ( Supplementary Figure 9A ). Patients with high TREM2 and TTC39A expression show higher levels of CD274, while those with high ANXA2 expression have lower CD274 levels ( Supplementary Figure 9B ). Additionally, patients with high GDF15 expression have lower CD8 expression levels ( Supplementary Figure 9C ). These findings emphasize the potential negative impact of the expression levels of these genes on the efficacy of immunotherapy in NASH-associated HCC. High expression of ANXA2, TREM2, TTC39A, and GDF15 may lead to a more suppressive tumor microenvironment and enhanced immune escape, thereby reducing the efficacy of immunotherapy.