High-resolution dissection of two-step reprogramming from fibroblast to induced hepatocytes using scRNA-seq

Our two-step lineage reprogramming consists of a reprogramming step from fibroblasts to hHPLCs (hHPLCs reprogramming step) and a maturation step from hHPLCs to hiHeps (hiHeps maturation step) (Fig. 1A). We observed that critical cellular transformations happened at the late stages of hHPLCs reprogramming step and the following hiHeps maturation step, the critical time points for establishment of hepatocyte-specific GRNs. Therefore, we selected samples from the late stages of hHPLCs reprogramming at day 38 (R-D38), day 54 (R-D54), and day 63 (R-D63), as well as samples from hiHeps maturation step at day 0 (M-D0), day 2 (M-D2), day 4 (M-D4), and day 7 (M-D7) for dissecting cell fate transition (Fig. 1A). Additionally, we included a starting fibroblast sample (R-D0) and a primary human hepatocyte (PHHs) sample as controls for gene expression analysis. These samples were processed to analyze reprogramming trajectories by measuring single-cell gene expression profiles using the Chromium system (10 × Genomics) (Fig. 1A). Sequencing data were collected from 21,304 individual cells across 9 samples, with a median of 6070 genes and 37,291 UMI per cell (Fig. S1A and S1B). Cells with less than 10% mitochondrial content and more than 2000 detected genes (nFeature) were selected for further analysis (Fig. S1C and S1D).

Fig. 1

High-Resolution Dissection of Two-Step Reprogramming from Fibroblast to induced Hepatocytes Using scRNA-Seq. A Schematic diagram of the scRNA-seq analysis strategy during the two-step lineage reprogramming process. In the reprogramming step, samples were collected at day 0 (R-D0), day38 (R-D38), day 54 (R-D54), and day 63 (R-D63). In the maturation step, samples were collected at day 0 (M-D0), day 2 (M-D2), day 4 (M-D4), and day 7 (M-D7). Additionally, a primary human hepatocytes (PHHs) sample was collected. All samples were processed using 10 × Genomics single-cell capture and sequencing. B UMAP of all 21,304 individual cells during the whole reprogramming process, colored by indicated time points. C UMAP of all 21,304 individual cells during the whole reprogramming process, colored by cell clusters. D Dot plots showing the expression of marker genes for 6 cell clusters. E–G Visualizations of scRNA-seq data, projection the expression score of representative fibroblast markers, hepatic progenitor markers and hepatocyte markers onto a UMAP plot. H Heatmap showing 6 groups of DEGs in the 6 cell clusters. Each column represents a single cell. Representative genes from each group of DEGs were listed to the right. I GO enrichment analysis of the six groups of DEGs

We projected all single cells onto a Uniform Manifold Approximation and Projection (UMAP) plot, identifying 6 transcriptionally distinct clusters (Fig. 1B and C). To classify the main cell types, we annotated each cluster based on marker gene expression (Fig. 1D–G). Sample R-D0 fell into cluster 0 identified as fibroblasts, characterized by the expression of mesenchymal genes such as DCN, THY1, FGF5, and STC2. Samples R-D38, R-D54, R-D63 fell into cluster 1, and sample M-D0 into cluster 2. Clusters 1 and 2 were identified as hepatic progenitor-like cells, marked by the expression of DLK1, EPCAM, and CD24. Samples M-D2 and M-D4 in cluster 3 represented immature hepatocytes, showing both progenitor-like markers and low levels of hepatocyte functional genes like ALB, CYP3A4, CYP2C19, UGT2B7, and UGT2B15. Samples M-D7 fell into cluster 4 corresponding to mature hepatocytes, distinguished by high expression of hepatocyte functional genes. Finally, PHHs were grouped into cluster 5.

To dissect the transcriptomic differences among these samples, we further analyzed the differentially expressed genes (DEGs) (Fig. 1H) and their Gene Ontology (GO) enrichment across different cell clusters (Fig. 1I). We identified 6 distinct groups of DEGs that define specific gene signatures for each cluster (Fig. 1H). Group 0 (G0) was enriched with fibroblast genes such as COL1A1, COL1A2, and VIM, predominantly expressed in cluster 0 (Fig. 1H). Groups 1 and 2 (G1 and G2) contained epithelial and hepatic progenitor markers like CD24, KRT8, KRT18, TJP1, and HES1, expressed in clusters 1 and 2, and were associated with processes such as lipid metabolism and aerobic respiration (Fig. 1H and I). Group 3 (G3) was enriched with hepatocyte-related genes including TTR, APOA1, APOC3, HMGCS2, and PPARG, expressed in cluster 3, which involved cell migration and actin filament organization (Fig. 1H and I). Group 4 (G4) included hepatocyte functional genes such as ALB, CYP3A4, UGT2B15, AZGP1, F2, and GLUL, expressed in cluster 4, and was linked to metabolic and catabolic processes like fatty acid metabolic process, small molecule catabolic process and steroid metabolic process (Fig. 1H and 1I). These processes were also similarly represented in Group 5 (G5) genes expressed in PHHs (Fig. 1G). GO enrichment analysis highlighted the significance of metabolic and catabolic pathways, such as fatty acid metabolic processes and small molecule catabolic processes, throughout the reprogramming and maturation stages (Fig. 1I).

Two-step lineage reprogramming recapitulates a trajectory of liver development

Next, we investigated whether the two-step reprogramming process recapitulates a trajectory of liver development in vitro. Wesley et al. recently unveiled a single-cell transcriptome atlas of human liver development, identifying specific liver development stages as hepatoblasts (HB1 and HB2), fetal hepatic cells (FH1 and FH2), and adult hepatic cells (AH), along with key marker genes for each stage [26, 27]. We found that the late stage of our two-step reprogramming process aligns with liver development, as evidenced by corresponding marker gene expression patterns with liver samples across various stages of liver development. Specifically, hHPLCs reprogramming samples expressed hepatoblast (MAP2K2, NPW, BRI3, BAAT, GSTP1, DLK1, APOM, and AMN) and fetal hepatic markers (AFP, CYP3A7, F2, AHSG, APOA1, AGT, ALB, ANGPTL3, APOB, GC, VTN, and AKR1C1). This expression profile indicates that hHPLCs exhibit a mixed characteristic of HB and FH identity. hiHeps maturation samples expressed adult hepatic markers (NNMT, AZGP1, TAT, CES1, ADH1B, HP, and SAA1) (Fig. 2A). Additionally, the expression dynamics of key genes in the late stage of our two-step reprogramming process mirrored those observed in liver development: the reprogramming of hHPLCs from R-D38 to M-D0 resembled the developmental stages from 7 weeks (7W) to 12 weeks (12W) of the human liver, while the maturation of hiHeps from M-D0 to M-D7 paralleled the stages from 13 weeks (13W) to 17 weeks (17W) (Fig. 2B).

Fig. 2

Two-Step Lineage Reprogramming Recapitulates a Trajectory of Liver Development. A Bubble chart displaying marker gene expression in the late stage of two-step reprogramming samples. Genes are colored by their expression stage, namely hepatoblast stages 1 and 2 (HB1, HB2), foetal hepatocyte stages 1 and 2 (FH1, FH2) and adult hepatocytes (AH), as defined in a previous study of human liver development. B Heatmaps displaying marker gene expression in the late stage of two-step reprogramming samples and human liver samples. C, D Pseudotime trajectory analysis of single cells throughout the late stage of two-step reprogramming, colored by indicated time points. Trajectory reconstruction reveals three branches: pre-branch (before bifurcation) and two terminal branches (after bifurcation). E, F Pseudotime trajectory analysis of single cells during the late stage of two-step reprogramming, colored by cell state. G Gene expression heatmap of the top differentially expressed genes (DEGs), cataloged into 3 clusters in a pseudo-temporal order. Two terminal branches of state 2 and state 3 are shown on the top right and left, respectively. H The expression dynamics of the top DEGs categorized into three major clusters in pseudotime trajectory in state 2 cells. Thick lines indicate the average gene expression patterns in each cluster. I Gene signatures and expression dynamics of representative liver developmental genes in each gene cluster as shown in H

To elucidate the trajectory of two-step reprogramming process, we used Monocle 2 to arrange cells in pseudotemporal order (Fig. 2C and D). Samples R-D38, R-D54, R-D63 and M-D0 were positioned in state 1. These cells then branched into two groups, designated as state 2 and 3. State 2 primarily included sample M-D7, while state 3 mainly consisted of samples M-D0 and M-D4 (Fig. 2E and F). According to the pseudotime analysis, state 2 represented the final cell type, illustrating the transition from fibroblasts to hHPLCs and ultimately to mature hepatocytes (Fig. 2F). We then investigated the sequential changes of 3350 DEGs of these cells along pseudotime, which were categorized into 3 clusters (Fig. 2G). We focused on the expression patterns of these 3 clusters as the cells progressed from state 1 to 2 (Fig. 2H and I). Notably, cluster 1 exhibited a stable expression pattern of FH marker genes (APOA1, AGT, GC, and ALB), cluster 2 showed an up-regulation pattern for AH marker genes (NNMT, SAA1, TAT, and HP), and cluster 3 demonstrated a down-regulation pattern for HB marker genes (GSTP1, NPW, and MAP2K2), consistent with their expression during liver development (Fig. 2H and I).

We further performed principal component analysis (PCA) on our samples. The cells progressed along the reprogramming timeline in a sequential order of R-D38, R-D54, R-D63, M-D0, M-D2, M-D4 and M-D7 (Fig. S2A). The major differences in gene expression were captured in PC1, which included genes associated with various metabolic processes, highly similar to the related processes in mouse liver development from E17.5 to P60 [28] (Fig. S2B). These findings support the notion that the late stage of two-step lineage reprogramming mirrors the natural development process of liver.

CD24 and DLK1 marks two different hepatic progenitor populations

To gain a detailed understanding of the two-step reprogramming process in the late stage, we investigated the hHPLCs reprogramming step and the hiHeps maturation step separately. We first explored the molecular and cellular dynamics of hepatic progenitor reprogramming by visualizing the reprogramming of hHPLCs using UMAP, which revealed 6 distinct clusters (Fig. 3A and B). We focused on the transmembrane proteins CD24 and DLK1, both critical for liver development and the maintenance of progenitor cells [29, 30]. Our analysis revealed that cluster 6, comprising M-D0 cells, exhibited enriched expression of CD24 as well as other hepatic progenitor markers such as ITGA6, NQO1, MOK and CDH1. Meanwhile, cluster 2, comprising R-D54 and R-D63 cells, showed enriched expression of DLK1 and other hepatic progenitor markers such as ITGB1, CTNNB1, ITGA6 and ANPEP [31] (Fig. 3C). Other clusters showed low expression of progenitor markers, suggesting a pre-progenitor stage containing R-D38 cells and subpopulation of R-D54 and R-D63 cells.

Fig. 3

CD24 and DLK1 Marks Two Different Hepatic Progenitor Populations. A UMAP of individual cells during the first step of reprogramming process, colored by indicated time points. B UMAP of individual cells during the first step of reprogramming process, colored by cell clusters. C Dot plots showing the expression levels of marker genes for each cluster. D Pseudotime trajectory analysis of reprogramming cells, colored by indicated time points. E, F. Pseudotime trajectory analysis of reprogramming cells, colored according to pseudotime progression. G, H Pseudotime trajectory analysis of reprogramming cells, colored by cell state. I Violin plot displaying the expression of representative hepatic progenitor genes in state 4 and state 5 cells. J Expression pattern of CD24 and DLK1 along pseudotime trajectory. K Schematic diagram of sorting CD24 and DLK1 positive hHPLCs followed by their maturation in 2C medium. L RT-qPCR analysis of hepatocyte gene expression (n = 2). The data are presented as the mean ± SD. Significance was assessed using one-way ANOVA. ***p < 0.001; **p < 0.01; *p < 0.05

To further elucidate the trajectory of progenitor reprogramming, we ordered cells in a pseudotemporal manner using Monocle 2 (Fig. 3D–H). The hHPLCs reprogramming samples were arranged along the pseudotime in the sequence of R-D38, R-D54, R-63, and M-D0 (Fig. 3D–F). These cells were grouped into five states, bifurcating into states 4 and 5, which represented two major progenitor populations at the final progenitor stage (Fig. 3G and H). Cells in the terminal of state 5 expressed key hepatic progenitor genes such as CD24, GSTA1, CKB, and PATJ, indicating successful reprogramming into hepatic progenitors (Fig. 3I). In contrast, cells at the alternate terminal exhibited high expression levels of FABP1, SLC35F1, SNHG5, and KRT19, indicative of unsuccessful reprogramming with a biased program (Fig. 3I). Notably, CD24 was enriched in the successful branch, confirming it as a reliable marker for the final hepatic progenitor population (Fig. 3J). Interestingly, DLK1 was upregulated prior to CD24 upregulation but was not enriched in the terminal branch, suggesting that DLK1-positive cells may represent an earlier stage of progenitor cells (Fig. 3J).

To determine whether CD24 and DLK1 can serve as surface markers for authentic hepatic progenitors during functional hepatocyte maturation, we sorted DLK1-positive and CD24-positive hHPLCs and cultured them in 2C maturation medium for an additional 7 days (Fig. 3K). Both CD24-positive and DLK1-positive cells exhibited significantly enhanced hepatocyte maturation potential compared to unsorted hHPLCs, as evidenced by the upregulation of of key functional hepatocyte markers, including ALB, CYP3A4, CYP2C9, CYP2C19, CYP2C8, CPS1, OTC, ASS, ARG1, NTCP, PXR, HNF1A, CEBPA, and UGT1A1 (Fig. 3L). Notably, while CD24-positive cells exhibited higher expression of ALB, most other functional genes were comparably expressed between the two sorted populations (Fig. 3L). These results suggest that CD24 and DLK1 define different states of progenitor cells, and both markers may be useful for enriching progenitor cells to enhance the efficiency of hepatocyte maturation.

Lipid metabolism mediates reprogramming of hepatic progenitors into functional hepatocytes

To explore the maturation of hepatic progenitors into functional hepatocytes in 2C maturation medium, we next analyzed maturation samples at different maturation stages. Maturing hepatocytes were grouped into 3 clusters: M-D0 formed cluster 0, M-D2 and M-D4 formed cluster 1, and M-D7 formed cluster 2 (Fig. 4A and B). Differential expression analysis revealed that cluster 0 was enriched with progenitor and proliferative markers such as CD24, CDK1, CDK4, PCNA, and KRT18. Cluster 1 expressed basic hepatocyte markers including APOE, TTR, APOA1, and ACSL1, while cluster 2 was enriched with functional hepatocyte genes like ALB, CYP3A4, CYP2C19, and UGT2B15 (Fig. 4C and D), reflecting progressive maturation over the 2C culture timeline. GO enrichment analysis further highlighted fatty acid metabolism as a key process in hepatocyte maturation (Fig. 4E).

Fig. 4

Lipid Metabolism Mediates Reprogramming of Hepatic Progenitors into Functional Hepatocytes. A UMAP of individual cells during the second step of the maturation process, colored by the indicated time points. B UMAP of individual cells during the second step of the maturation process, colored by cell clusters. C Dot plots showing the expression of marker genes for each cluster. D Heatmap showing three groups of DEGs of the 3 cell clusters. Each column represents a single cell. Representative genes for each group of DEGs were listed on the right. E GO enrichment analysis of three groups of DEGs. F Pseudotime trajectory analysis of reprogramming cells, colored by the indicated time points. G, H Pseudotime trajectory analysis of reprogramming cells, colored by cell states. I Violin plot displaying the expression of representative hepatocyte genes in state 2 and 3 cells. J RT-qPCR analysis of functional hepatocyte gene expression (n = 2). The data are presented as the mean ± SD. Significance was assessed using one-way ANOVA. ***p < 0.001; **p < 0.01; *p < 0.05

Furthermore, Monocle pseudotime analysis revealed two distinct cell populations (Fig. 4F–H). M-D0 fell primarily into state 1; M-D2, M-D4 and M-D7 divided into two branches (Fig. 4F–H). One population (state 2) expressed major functional hepatocyte genes such as UGT2B15, SLC7A2, CYP3A4, and CYP3A5, indicating successful hepatocyte maturation (Fig. 4I). The other population (state 3) expressed genes such as AMBP, IGFBP2, MUC13, and AMBRA1, characteristic of immature hepatocytes and intestinal cells, representing the failed maturation branch (Fig. 4I). Analysis of 2626 DEGs post-bifurcation (Fig. S3A) revealed that cells in the successful branch had elevated levels of genes involved in metabolic processes, particularly fatty acid metabolism, consistent with the GO analysis (Fig. S3B-S3D, and Fig. 4E). These gene expression patterns suggest that lipid metabolism plays a critical role in hepatocyte maturation.

To validate the role of lipids in hepatocyte maturation, we supplemented the 2C maturation medium with chemically defined lipid concentrate (CDL) during hiHeps maturation. CDL supplementation significantly enhanced hiHeps maturation, evidenced by the upregulation of functional hepatocyte genes such as CYP2C9, CYP2C19, CYP1A2, UGT2B15, UGT2B7, CPS1, MRP2, and HNF4A (Fig. 4J). These findings underscore the essential role of lipid metabolism in promoting hepatocyte maturation.

HNF4A and HHEX regulate hepatocyte and intestinal cell fate during maturation step in two-step reprogramming

As the maturation process indicated an intestinal cell fate in the hiHeps, we further investigated the cell identity at the final stage (M-D7). M-D7 cells were divided into two sub-populations (Fig. 5A). Cluster 0 expressed functional hepatocyte genes, including CYP3A4, CYP2C19, CYP2C9, CYP2D6, CYP2C8, UGT1A1, UGT2B7, UGT2B15, and PLG (Fig. 5B-5C and S4A), while cluster 1 expressed intestinal cell genes like AGR2, LGALS3, LCN2, S100A6, S100P, MUC5AC, MUC5B, TTF1, and VSTM2L (Fig. 5B, C and S4B). Both clusters expressed common hepatocyte markers, such as ALB, TTR, FABP1, and SERPINA1, indicating a subpopulation of hiHeps acquired a mixed hepatocyte and intestinal cell fate (Fig. 5B, C and S4C). Differential expression and GO enrichment analysis confirmed that cluster 0 was enriched with genes involved in oxidation and respiration processes, while cluster 1 was enriched with metabolic process genes typical of intestinal cells (Fig. 5D and E). CellNet analysis further validated these findings, showing that cluster 0 had a significantly higher liver score, while cluster 1 had a higher intestinal score (Fig. S4D and S4E).

Fig. 5

HNF4A and HHEX Regulate Hepatocyte and Intestinal Cell Fate During Maturation Step in Two-Step Reprogramming. A UMAP of individual cells of the M-D7 sample. B UMAP visualizations showing the expression of representative cluster 0 marker genes (left), cluster 1 marker genes (middle), and shared genes of cluster 0 and 1 (right). C Dot plots showing the expression of marker genes. D Heatmap showing the two groups of DEGs in the 2 cell clusters. Each column represents a single cell. Representative genes for each group of DEGs were listed on the right. E GO enrichment analysis of the two groups of DEGs. F RT-qPCR analysis of hepatocyte gene expression in HNF4A knockdown hiHeps (n = 2). G RT-qPCR analyses of intestine gene expression in HHEX knockdown hiHeps (n = 2). The data are presented as the mean ± SD in figure H and I. Significance was assessed using one-way ANOVA. ***p < 0.001; **p < 0.01; *p < 0.05

To investigate the biased cell fate determination, we examined the correlation between five reprogramming TFs that are expressed across various endoderm-derived tissues and their influence on the determination of hepatocyte and intestinal lineages. Cluster 0 was closely correlated with HNF4A, while cluster 1 was correlated with HHEX (Fig. 5C). HNF4A is essential for liver development and functional hepatocyte differentiation [32, 33], whereas HHEX is a marker of the definitive endoderm, crucial for early liver development [34]. However, the contrasting roles of these TFs in hepatocyte reprogramming remain unclear.

To confirm the roles of HNF4A and HHEX in hiHeps reprogramming, we knocked down these genes in hHPLCs and matured them in 2C medium for an additional 7 days to generate hiHeps. Knockdown of either HNF4A and HHEX reduced the expression of ALB and other key markers (Fig. 5F and G). Specifically, HNF4A knockdown led to a significant reduction in functional hepatocyte markers like ALB, CYP3A4, MRP2, CYP2C9, and CYP2C19, with minimal effect on intestinal genes like MUC5AC, AGR2, and FABP2 (Fig. 5F). Conversely, HHEX knockdown decreased the expression of intestinal genes MUC5AC, AQP8, EPHB2, DLFA, FABP2, CHGA, and CDX2, while only slightly affecting ALB expression (Fig. 5G). These results suggest that HHEX promotes an intestinal cell fate in hiHeps during maturation, whereas HNF4A is crucial for maintaining hepatocyte function.

Gene regulatory networks in the two-step reprogramming process

To gain insight into specific TF regulatory networks involved in the late stage of two-step reprogramming process, we re-mapped the 7 samples representing different stages of reprogramming and maturation, clustering them into four distinct groups (Fig. 6A and B). Single-cell regulatory network inference and clustering (SCENIC) analysis revealed that cluster 0 of early pre-progenitor samples from day 38 to day 63 (R-D38, R-D54 and R-D63) was regulated by TFs such as HES7 and TBX15 (Fig. 6C and D). Cluster 1 of hepatic progenitor samples (M-D0) was regulated by MYBL2 and SMARCA4, etc. (Fig. 6C and D). Cluster 2 representing immature hepatocyte samples (M-D2 and M-D4) was primarily regulated by HOXA10, HOXB7, PPARA, and NR1H4, etc., while cluster 3, which contained mature hepatocyte samples (M-D7) was regulated by RARA and FOXO3, etc. (Fig. 6C and D). The expression pattern of these TFs was shown in the feature plots (Fig. S5A).

Fig. 6

Gene Regulatory Networks in the Two-step Reprogramming Process. A UMAP of individual cells from the late stage of the two-step reprogramming samples, colored by indicated time points. B UMAP of individual cells from the late stage of the two-step reprogramming samples, colored by cell cluster. C Single-Cell Regulatory Network Inference and Clustering (SENIC) analysis of regulatory TFs of each cluster. Top ranked TFs are shown. D Dot plots of representative TFs for each cluster. E Transcriptional regulatory network showing the progression from pre-progenitor cells to mature hepatocytes. This network is divided into four subnetworks, each representing distinct cell states during hepatocyte development: pre-progenitor (red), progenitor (blue), immature hepatocyte (green), and mature hepatocyte (purple)

We also constructed a gene network based on pairwise correlations of the top DEG expressions during progressive cell fate transitions from R-38 to M-D7 (Fig. 6E). This analysis revealed four chronological subnetworks. The pre-progenitor subnetwork was closely linked to progenitor related genes, and the immature hepatocyte genes were tightly connected to the mature hepatocyte subnetwork (Fig. 6E). The sequential switch of transcriptional circuits underscored the close relationship between the pre-progenitor and progenitor subnetworks, which acted as a bridge linking immature hepatocytes subnetwork to that of mature hepatocytes, which consisted of key liver functional transcription factors such as RARA, XBP1, PPARA, CEBPB, and NR5A2 (Fig. 6E).