Publicly available single-cell RNA sequencing (scRNA-seq) datasets were obtained from the GEO and CNGB databases, including the CNGB Nucleotide Sequence Archive (CNP0002310) (Xu et al., 2024), Zenodo repository (Zenodo.6035873) (Ben-Moshe et al., 2022), and GEO accessions (GSE202379) (Gribben et al., 2024). All datasets were generated using the 10x Genomics Chromium platform.

Data analysis was performed using Seurat v5 in R. Highly variable genes (top 2000) were identified using the FindVariableFeatures function, followed by principal component analysis (PCA). Batch effects were corrected using the Harmony package prior to clustering. Cells were clustered with FindNeighbors and FindClusters. Filtering thresholds were set at nFeature_RNA >500 and percent. mt <10% to exclude low-quality cells and potential doublets.

Cell-type annotation was performed using established hepatocyte and biliary marker genes to identify hepatocyte clusters. Within annotated hepatocyte clusters, cells with detectable Sox9 expression (expression value >0) were defined as SOX9⁺ hepatocytes, whereas cells lacking detectable Sox9 expression were classified as SOX9^- hepatocytes. The proportion of SOX9⁺ hepatocytes was calculated for each sample and time point.

Differentially expressed genes (DEGs) between SOX9⁺ and SOX9^- hepatocytes were identified using FindMarkers with the criteria |log₂FC| >0.25 and FDR <0.05. Marker genes of SOX9⁺ hepatocytes were determined based on DEG analysis. Gene Set Enrichment Analysis (GSEA) for KEGG pathways was performed on all SOX9⁺ hepatocytes in each experimental model.

All SOX9⁺ hepatocytes were extracted and re-clustered within each model to identify distinct subpopulations and to characterize their transcriptional changes over time and across disease conditions. KEGG enrichment analysis was performed for each SOX9⁺ subcluster.

Temporal changes in the proportions of SOX9⁺ hepatocytes were quantified and visualized to delineate their activation kinetics during liver regeneration.

To characterize SOX9⁺ hepatocytes during liver regeneration after partial hepatectomy, we first compared the transcriptomes of SOX9⁺ and SOX9^- hepatocytes. Differential gene expression analysis revealed widespread transcriptional differences between the two populations, as visualized by a volcano plot (Figure 1a). A substantial number of genes were significantly upregulated in SOX9⁺ hepatocytes, indicating that these cells adopt a transcriptional program distinct from SOX9^- hepatocytes rather than representing minor transcriptional variation within the hepatocyte compartment.

We next examined the expression of representative lineage and functional markers across hepatocyte populations to define the cellular identity of SOX9⁺ hepatocytes. Dot plot analysis demonstrated that SOX9⁺ hepatocytes continued to express canonical hepatocyte markers, confirming their hepatocyte lineage origin (Figure 1b). Importantly, these cells did not express classical cholangiocyte markers, excluding biliary contamination. In contrast to SOX9^- hepatocytes, SOX9⁺ hepatocytes displayed increased expression of genes associated with progenitor-like and stress-responsive states, suggesting that SOX9⁺ hepatocytes represent a distinct hepatocyte state induced during regeneration rather than a separate epithelial lineage.

To gain functional insight into the SOX9⁺ hepatocyte state, we performed KEGG pathway enrichment analysis on genes upregulated in SOX9⁺ hepatocytes. This analysis revealed significant enrichment of pathways related to inflammatory and immune signaling, including TNF, NF-κB, and Toll-like receptor signaling pathways, as well as pathways associated with integrin signaling and focal adhesion (Figure 1c). These enriched pathways suggest that SOX9⁺ hepatocytes actively engage in inflammatory signaling and cell–matrix interactions, processes known to be critical for effective tissue remodeling and regeneration following liver injury.

We also investigated the temporal dynamics of SOX9⁺ hepatocytes following partial hepatectomy. Quantification of SOX9⁺ hepatocytes across post-operative time points revealed a rapid increase in their proportion during early regenerative stages, followed by a gradual decline as regeneration progressed toward completion (Figure 1d). This transient expansion indicates that SOX9⁺ hepatocytes represent a regeneration-associated hepatocyte population that emerges in response to surgical injury and diminishes as liver homeostasis is restored.

To further resolve the temporal organization of SOX9⁺ hepatocytes during liver regeneration, we visualized SOX9⁺ hepatocytes across post-hepatectomy time points using UMAP embedding (Figure 1e). SOX9⁺ hepatocytes segregated into distinct regions corresponding to different time points, with minimal overlap between early (8 h, D0) and later regenerative stages (D1–D7). This temporal separation indicates pronounced transcriptional remodeling of SOX9⁺ hepatocytes as regeneration progresses, rather than persistence of a static SOX9⁺ cell population. Notably, SOX9⁺ hepatocytes from early time points occupied discrete UMAP regions, whereas cells from later stages showed increased heterogeneity, suggesting diversification of SOX9⁺ hepatocyte states during regenerative progression.

To quantify the temporal dynamics of SOX9⁺ hepatocyte heterogeneity, we examined the composition of SOX9⁺ subclusters across post-hepatectomy time points. The relative abundance of individual SOX9⁺ subclusters changed markedly over time (Figure 1f). Early after hepatectomy, SOX9⁺ hepatocytes were dominated by a limited number of subclusters, whereas later stages exhibited a broader and more complex subcluster composition. This dynamic redistribution indicates that SOX9⁺ hepatocytes undergo progressive state transitions during regeneration, with distinct subpopulations emerging and resolving at different regenerative stages.

To characterize the molecular features underlying SOX9⁺ hepatocyte heterogeneity, we identified marker genes for each SOX9⁺ subcluster. Heatmap visualization revealed clear transcriptional differences among subclusters, with distinct gene expression patterns defining individual states (Figure 1g). Several subclusters were characterized by enrichment of genes associated with stress responses, inflammatory signaling, and cell–matrix interaction, whereas others displayed higher expression of metabolic or housekeeping genes. These findings demonstrate that SOX9⁺ hepatocytes comprise multiple transcriptionally distinct states rather than a homogeneous population.

To link transcriptional heterogeneity to functional differences, we performed KEGG pathway enrichment analysis for SOX9⁺ hepatocyte subclusters across time points (Figure 1h). Early regenerative stages were enriched for metabolic pathways, including bile acid metabolism and xenobiotic metabolism, consistent with retention of hepatocyte functional programs. In contrast, later stages exhibited enrichment of pathways related to cell cycle regulation, p53 signaling, ECM–receptor interaction, and cell adhesion, indicating activation of proliferative and tissue remodeling programs. These pathway shifts suggest a temporal transition of SOX9⁺ hepatocytes from metabolically active states toward proliferative and remodeling-associated states during regeneration.

To further illustrate temporal functional changes, we examined the expression dynamics of representative genes across post-hepatectomy time points. Expression of the hepatocyte metabolic marker Alb decreased during early regenerative stages and partially recovered at later time points (Figure 1i). In contrast, expression of Sox9 increased progressively, consistent with sustained progenitor-associated identity during regeneration. Inflammatory and proliferative markers, including Cxcl1 and Mki67, exhibited transient induction at early to intermediate stages, coinciding with peak regenerative activity. These gene expression dynamics reinforce the concept that SOX9⁺ hepatocytes undergo stepwise functional reprogramming during liver regeneration.

To characterize SOX9⁺ hepatocytes in the context of APAP, we compared the transcriptional profiles of SOX9⁺ and SOX9^- hepatocytes. Differential expression analysis revealed extensive transcriptional differences between the two populations, as visualized by a volcano plot (Figure 2a). A large number of genes were significantly upregulated in SOX9⁺ hepatocytes, indicating that these cells adopt a distinct transcriptional program in response to APAP-induced injury.

To determine the cellular identity of SOX9⁺ hepatocytes after APAP injury, we examined the expression of representative lineage markers across epithelial cell populations. Dot plot analysis demonstrated that SOX9⁺ hepatocytes continued to express canonical hepatocyte-associated genes, while lacking expression of classical cholangiocyte markers such as Krt19 (Figure 2b). These results exclude biliary contamination and confirm that SOX9⁺ hepatocytes represent a hepatocyte-derived population that acquires a distinct transcriptional state following injury.

To gain insight into the functional programs associated with the SOX9⁺ hepatocyte state, we performed KEGG pathway enrichment analysis on genes upregulated in SOX9⁺ hepatocytes. This analysis revealed significant enrichment of pathways related to inflammatory and stress responses, including integrin signaling, TNF signaling, MAPK signaling, and other injury-associated pathways (Figure 2c). These findings indicate that SOX9⁺ hepatocytes engage inflammatory and stress-responsive signaling programs during APAP-induced liver injury.

We next examined the temporal dynamics of SOX9⁺ hepatocytes following APAP administration. Quantification of SOX9⁺ hepatocytes across time points revealed a rapid increase in their proportion during the acute injury phase, followed by a gradual decline during the recovery phase (Figure 2d). This transient expansion suggests that SOX9⁺ hepatocytes represent an injury-induced hepatocyte state that emerges rapidly in response to toxic injury and resolves as liver homeostasis is restored.

To further resolve heterogeneity within the SOX9⁺ hepatocyte population, we analyzed the composition of SOX9⁺ hepatocyte subclusters across time points following APAP injury. The relative abundance of individual subclusters changed markedly over time, with specific subclusters predominating during early injury stages and others emerging during later recovery phases (Figure 2e). These results demonstrate dynamic remodeling of SOX9⁺ hepatocyte states throughout the injury–recovery continuum.

To link transcriptional heterogeneity to functional differences, we performed KEGG pathway enrichment analysis across SOX9⁺ hepatocyte subclusters over time. Early time points following APAP injury were characterized by enrichment of stress and inflammatory pathways, whereas later stages showed gradual reactivation of metabolic pathways, consistent with injury resolution and functional recovery (Figure 2f). These results indicate a temporal shift in functional programs of SOX9⁺ hepatocytes during the course of injury and repair.

Finally, we examined the temporal expression patterns of representative genes to illustrate functional changes in SOX9⁺ hepatocytes following APAP injury. Expression of the hepatocyte metabolic marker Alb decreased during the acute injury phase and partially recovered at later time points. In contrast, Sox9 expression was maintained or modestly increased, while inflammatory and proliferative markers such as Cxcl1 and Mki67 exhibited transient induction during early to intermediate stages (Figure 2g). These expression dynamics highlight stepwise functional reprogramming of SOX9⁺ hepatocytes during acute injury and subsequent repair.

We compared the transcriptomic profiles of SOX9⁺ and SOX9^- hepatocytes to characterize SOX9⁺ hepatocytes in metabolic liver disease. Distinct gene expression patterns between the two groups are shown in the volcano plot (Figure 3a). To clarify the cellular identity of SOX9⁺ hepatocytes, we examined the expression of representative hepatocyte and cholangiocyte marker genes. Dot plot analysis showed that SOX9⁺ hepatocytes expressed hepatocyte-associated markers while lacking expression of classical cholangiocyte markers such as KRT19 (Figure 3b). In contrast, cholangiocytes displayed robust expression of biliary markers.

To gain insight into the functional programs associated with SOX9⁺ hepatocytes, we performed KEGG pathway enrichment analysis on genes upregulated in SOX9⁺ hepatocytes. This analysis revealed significant enrichment of pathways related to inflammatory signaling and stress responses, including MAPK signaling, TNF signaling, bacterial invasion–related pathways, and focal adhesion (Figure 3c). These findings suggest that SOX9⁺ hepatocytes engage stress-responsive and cell–matrix interaction programs, consistent with a non-homeostatic hepatocyte state.

We next quantified the proportion of SOX9⁺ hepatocytes across disease stages. SOX9⁺ hepatocytes were rare in healthy controls and NAFLD, accounting for less than 1% of hepatocytes. Their proportion increased modestly in NASH without cirrhosis and exhibited a pronounced expansion in NASH with cirrhosis and end-stage disease, reaching approximately 2% of the hepatocyte population (Figure 3d). These results indicate a stage-dependent accumulation of SOX9⁺ hepatocytes during liver disease progression, with a sharp increase coinciding with advanced disease stages.

To further resolve heterogeneity within the SOX9⁺ hepatocyte population, we examined the composition of SOX9⁺ hepatocyte subclusters across disease stages. Stacked bar plots revealed substantial shifts in subcluster composition during disease progression (Figure 3e). Early disease stages were dominated by a limited number of subclusters, whereas advanced stages showed increased representation of additional subclusters, indicating progressive remodeling of SOX9⁺ hepatocyte states during disease progression.

We next characterized the transcriptional features underlying SOX9⁺ hepatocyte heterogeneity by identifying marker genes for each subcluster. Heatmap analysis demonstrated clear differences in gene expression patterns among SOX9⁺ hepatocyte subclusters (Figure 3g). Distinct subsets of genes defined individual subclusters, indicating that SOX9⁺ hepatocytes comprise multiple transcriptionally distinct states rather than a homogeneous population.

To investigate functional changes associated with disease progression, we analyzed KEGG pathway enrichment dynamics across disease stages. Early stages were characterized by enrichment of metabolic pathways, including amino acid metabolism, fatty acid metabolism, and xenobiotic metabolism (Figure 3f). In contrast, advanced disease stages showed strong enrichment of pathways related to inflammatory signaling, cytokine–cytokine receptor interaction, cell cycle regulation, ECM–receptor interaction, and cancer-associated pathways, indicating a functional transition toward inflammatory, proliferative, and tissue-remodeling programs.

Finally, we examined the expression dynamics of representative genes across disease stages to illustrate functional reprogramming of SOX9⁺ hepatocytes. Expression of the hepatocyte metabolic marker ALB progressively decreased with disease severity, whereas SOX9 expression was maintained. In parallel, inflammatory and signaling-associated genes such as STAT3 and ECM-related genes such as ITGA6 were upregulated in advanced disease stages (Figure 3h). These expression patterns highlight a progressive shift of SOX9⁺ hepatocytes from a metabolically active hepatocyte state toward an inflammatory and tissue-interactive phenotype during disease progression.