The scRNA-seq data of three human cirrhotic tissue samples and three healthy liver tissue samples, two mouse healthy liver mesenchymal cell samples and two fibrotic (following chronic CCl${}_4$-induced liver injury) liver mesenchymal cell samples, and one Lin-negative cell sample of the bilio-vascular tree from healthy mouse liver were acquired from the Gene Expression Omnibus (GEO) dataset with accession numbers GSE136103 [4], GSE137720 [20], GSE163777 [14]. Raw data (sra files) were obtained from the NCBI Sequence Read Archive (SRA) database with accession numbers SRP218975, SRP222529, SRP299106. We then extracted reads files from sra files using the SRA Tool Kit v3.0.0 (Bethesda, MA, USA) fastq-dump Tool and obtained the fastq files. We aligned reads to the GRCh38 and mm10 reference genomes as appropriate for the input dataset using the Cell Ranger Single-Cell Software Suite v7.0.0 (Pleasanton, CA, USA) from 10X Genomics, with the functions “cellranger mkfastq” and “cellranger count”.

The R package Seurat v4.1.0 (Satija lab, New York, NY, USA) [21] was used to process human liver scRNA-seq data (GSE136103). Genes expressed in fewer than three cells in a sample, cells that expressed fewer than 300 genes, and more than 6000 genes, and cells with mitochondrial gene content $>$20% and ribosomal gene content $>$50% of the total unique molecular identifier (UMI) count were excluded. After quality control, an expression matrix comprising 52,669 cells and 30,947 genes was obtained. Next, data were normalised using the default parameters of the log-normalisation method, and the top 2000 highly variable genes were identified using the FindVariableFeatures function. Z-score transformation of gene expression was performed using the ScaleData function. The following dimensionality reduction analysis was performed using the RunPCA function. Clustering was based on the 20 most significant principal components (PCs). Then, the batch effect was removed by using R package Harmony v0.1.0 (Center for Data Sciences, Brigham and Women’s Hospital, Boston, MA, USA) [22]. Different resolutions were set, and the clustering effect was observed using clustree. Finally, the resolution was set to 0.6, and uniform manifold approximation and projection (UMAP) visualisations were constructed using the same number of PCs as the associated clustering.

We focused on mesenchymal cells, therefore, we extracted subpopulations of mesenchymal cells (cluster 9 and 15) from healthy and cirrhotic samples. Because we were more concerned with the dynamic evolution of mesenchymal cells in cirrhosis, we used the FindClusters function to sub-cluster these mesenchymal cells from healthy and cirrhotic liver separately. Dimensionality reduction analysis was performed using the RunPCA function, and sub-clustering was based on the nine most important key components. The resolution was set to 0.9 (healthy) and 1.2 (cirrhotic). After subclustering, healthy and cirrhotic mesenchymal cells were visualised with t-distributed stochastic neighbor embedding (t-SNE) using the nine most important PCs.

Mouse hepatic mesenchymal cells from two healthy and two fibrotic liver (GSE137720) were analysed using Seurat. Low quality cells (300 genes/cell, $>$4500 genes/cell, 3 cells/gene, $>$7.5% mitochondrion genes, and $>$25% ribosomal genes) were filtered out from the dataset. After quality control, 12,538 and 10,706 mesenchymal cells were obtained from healthy and fibrotic mouse livers, respectively. Normalisation, z-score transformation, dimensionality reduction, and clustering analysis were similar to those of the human liver scRNA-seq process, except that the resolution was set to 0.5. The t-SNE visualisations were constructed using the nine most significant PCs consistent with the clustering.

The scRNA-seq data of mouse Lin-negative cells of bilio-vascular tree (GSE163777) was analysed in a similar manner. Low-quality cells (300 genes/cell, $>$6000 genes/cell, 3 cells/gene, $>$10% mitochondrion genes, and $>$25% ribosomal genes) were filtered out from the dataset. After normalisation, z-score transformation, and dimensionality reduction, clustering was performed using the 20 most significant PCs. Finally, the resolution was set to 0.6, and t-SNE visualisations were constructed using the 20 most significant PCs consistent with the clustering. After cell type annotation, we extracted mesenchymal cells for sub-clustering. We used RunPCA to redimensionalise the extracted mesenchymal cells and clustered the 14 most important PCs. Finally, the resolution was set to 0.5, and t-SNE was used for visualisation.

After obtaining cell clusters, the FindAllMarkers function was used to search for differentially expressed genes in the clusters. We set the log fold-change of the average expression between the two clusters (avg_logFC) $>$0.5. The p value after Bonferroni correction was 0.05. Genes were ranked in ascending p values, differentially expressed in each cluster, and cell types were manually identified using liver cell marker genes from previous studies [4, 14, 20] and marker genes in PanglaoDB, CellMarker, and the Human cell landscape database.

The differentiation trajectory of the mesenchymal subpopulations was constructed using the R package Monocle v2.24.1 (Cole Trapnell’s lab, Seattle, WA, USA) [23]. The newCellDataSetFirst function was used to build a Monocle object. We used dispersionTable (mean_expression $\ge$0.1 & dispersion_empirical $\ge$1 * dispersion_fit) to select highly variable genes as ordering genes and used in Monocle for clustering. We ran reduceDimension (max_components = 2, method = ‘DDRTree’) with t-SNE as the reduction method, and ran orderCells under default parameters. Finally, a pseudo-time differentiation trajectory was plotted using the cell trajectory function.

We verified the Monocle trajectory and its directionality using the velocyto Python package v 0.17.17 (La Manno Lab, Stockholm, Sweden) [24] to estimate the cell velocity from the spliced and unspliced mRNA content. We used velocyto to convert the cellranger result files (output directory) into loom files and merged the loom files of all samples. The merged loom files were used as input for scvelo v0.2.4 (Volker Bergen, Munich, Germany) Python pipeline [25]. The calculation of RNA velocity values for each gene in each cell and embedding of the RNA velocity vector in t-SNE were performed using scvelo and finally visualised with Python v3.8 (Python Software Foundation, Delaware, DE, USA).

The R package CellChat v1.5.0 (Department of Mathematics, University of California, Irvine, CA, USA) [26] was used for cell communication analysis, and the Seurat object was exported as the input to CellChat. After the construction of the CellChat object using the createCellChat function, the ligand-receptor interaction database CellChatDB.human (or CellChatDB.mouse) was set up. The computeCommunProb function was used to calculate the communication probability with default parameters and infer the CellChat network. Function filterCommunication was used to filter groups with fewer than 10 cells. The function subset Communication was then used to save the prediction results of cell communication in the form of a data frame. The computeCommunProbPathway function was used to summarise the communication probability of all ligand-receptor interactions related to each signalling pathway, to calculate the communication probability at the signal pathway level. The integrated cell communication network was evaluated using the aggregateNet function by calculating the number of links or summarising the communication probability. The results of cell-cell communication were displayed using netVisual heatmap function.

The Python package pySCENIC v0.12.0 (VIB Center for Brain Disease Research, Laboratory of Computational Biology, Leuven, Belgium) [27] was used for transcription factor analysis. First, the gene expression matrix was exported from R and saved as “csv” files, and then the “csv” file was converted into a loom file in Python. The pyscenic grn, pyscenic cistarget, and pyscenic AUCell functions were run sequentially for transcription factor analysis and the results were visualised.

The human hepatic stellate cell line (LX-2) was provided without Mycoplasma contamination from Professor Wang Qingqing’s research group at Zhejiang University School of Medicine and was verified as LX-2 cells by the STR technology of Wuhan Punosei Life Technology. To maintain the cells in the same state before dosing, LX-2 cells in the logarithmic growth phase were starved for 24 h (the medium did not contain foetal bovine serum or penicillin/streptomycin). TGF-$\beta$1 at 50 ng/mL was then added for 24 h to construct the liver injury model (stimulate period), and no drug treatment was added to the blank control group (quiescent period). After 24 h, the LX-2 cells were left without treatment for 7 days (recovery period).

The stimulate, control, and recovery period samples were repeated three times, and RNA extraction, quality control, library construction, and RNA-seq were performed on nine samples. Total RNA was extracted using TRIzol reagents, RNA integrity was assessed using 1% agar-gel, and RNA concentration and purity were quantitatively and qualitatively analysed using a nanoophotometer spectrophotometre. RNA volume was 3 µg/ sample, the RNA-seq library was constructed using NEBNext UltraTM RNA Library Prep Kit for Illumina, and index codes were added to the attribute sequence of each sample. Following the manufacturer’s instructions, we used the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, CA, USA) , to cluster samples using index codes on the cBot Cluster Generation System. After clustering, we sequenced the library using the Illumina platform and obtained 150 bp paried-end sequences.

First, FastQC was used to evaluate Illumina reads. The Fastp software was used to discard low-quality reads. Then, we downloaded the human reference genome (GRCh38) from the Ensembl database and used HISAT2 software for read alignment. Finally, the gene levels were quantitatively analysed using the featureCounts tool to obtain the final standardised matrix. The expression level of each gene was quantified as normalised fragments per kilobase of transcript per million mapped reads (FPKM).

Unless otherwise stated, statistical analyses were performed using Graphpad Prism v8.0 (GraphPad Software, San Diego, CA, USA). A paired t-test was used to compare the differences between the two groups. Differences were considered statistically significant at p 0.05.

A graphical overview of the study design is shown in Fig. 1A. To examine the heterogeneity of hepatic NPCs, we analysed the hepatic NPC scRNA-seq data obtained from the GEO database (GSE136103). Hepatic NPCs were isolated from three healthy and three cirrhotic human livers. NAFLD was the cause of liver fibrosis in one female patient, and alcohol was the cause of liver fibrosis in the other two male patients (Supplementary Fig. 1A). Leucocytes (CD45${}^{\mathrm{+}}$) or other NPC (CD45${}^{\mathrm{-}}$) fractions were FACS-sorted before scRNA-seq [4]. After filtering low-quality cells (Supplementary Fig. 1B,C), the batch effect was removed using Harmony, and cells of the same type, such as endothelial cells, clustered more distinctly (Supplementary Fig. 1D,E). After dimensionality reduction and clustering (clustree is shown in Supplementary Fig. 2A), 52669 NPCs revealed 24 populations (Fig. 1B), and the top marker genes from the 24 cell clusters were identified by Seurat (Supplementary Fig. 2B, Supplementary Table 1). NPCs were represented by a total of 11 distinct cell types (Fig. 1C), which correspond to T cell (TRAC${}^{\mathrm{+}}$IL7R${}^{\mathrm{+}}$), innate lymphoid cell (ILC; PRF1${}^{\mathrm{+}}$FGFBP2${}^{\mathrm{+}}$), Endothelial cell (VWF${}^{\mathrm{+}}$STAB2${}^{\mathrm{+}}$), mononuclear phagocyte (MP; C1QB${}^{\mathrm{+}}$MARCO${}^{\mathrm{+}}$), Epithelial cell (ALB${}^{\mathrm{+}}$KRT18${}^{\mathrm{+}}$), mesenchymal cell (COL3A1${}^{\mathrm{+}}$TAGLN${}^{\mathrm{+}}$ACTA2${}^{\mathrm{+}}$PDGFRB${}^{\mathrm{+}}$), B cell (CD79A${}^{\mathrm{+}}$MS4A1${}^{\mathrm{+}}$), plasma cell (FCRL5${}^{\mathrm{+}}$JSRP1${}^{\mathrm{+}}$), cycling cell (TOP2A${}^{\mathrm{+}}$MKI67${}^{\mathrm{+}}$), plasmacytoid dendritic cell (pDC; LILRA4${}^{\mathrm{+}}$PTCRA${}^{\mathrm{+}}$), mast cell (TPSAB1${}^{\mathrm{+}}$CPA3${}^{\mathrm{+}}$) (Fig. 1D). Because the parenchymal cells were not fully removed during centrifugation, a small cluster of epithelial cells remained (Fig. 1C). MPs, Mast cells, ILCs, and T cells displayed high expression levels of CD45 (PTPRC) (Supplementary Fig. 2C). Although all 11 major cell types were present in both healthy and cirrhotic tissues, the proportions of major cell types were different (Fig. 1E). Cell sorting approaches led to loss of fragile cells, particularly hepatocytes and cholangiocytes, and the relative enrichment of more robust NPCs, such as endothelial cells/MPs [28]. As a result, a large number of endothelial cells and MP cells were enriched in the NPCs clustering results, while only a small proportion of mesenchymal cells (the main mesenchymal cells, HSC, accounted for only 6% of the liver).

Fig. 1.

Single-cell transcriptomic profile of human liver. (A) Graphic overview of the study design. (B) Clustering 52,669 cells from three healthy and three cirrhotic human livers. (C) The 11 lineages of 52,669 human liver cells, inferred from the expression of marker genes. Right, annotation by injury condition. ILC, innate lymphoid cell; MP, mononuclear phagocyte; pDC, plasmacytoid dendritic cell. (D) UMAP plots of conserved hepatic cell markers. CD79A, CD79a Molecule; MKI67, Marker Of Proliferation Ki-67; VWF, Von Willebrand Factor; ALB, Albumin; KRT18, Keratin 18; PRF1, Perforin 1; TPSAB1, Tryptase Alpha/Beta 1; COL3A1, Collagen Type III Alpha 1 Chain; ACTA2, Actin Alpha 2, Smooth Muscle; TAGLN, Transgelin; PDGFRB, Platelet Derived Growth Factor Receptor Beta; C1QB, Complement C1q B Chain; FCRL5, Fc Receptor Like 5; TRAC, T Cell Receptor Alpha Constant. (E) Bar plot representing the relative contribution of cells from each donor (three healthy livers, three cirrhosis livers) for each cell type.

We focused on mesenchymal cells, and we isolated clusters 9 and 15 (2918 mesenchymal cells) for further analysis. In healthy liver, mesenchymal cells sub-clustered into eight distinct mesenchymal subpopulations (Fig. 2A). We identified the cells in clusters 1, 2, 3, and 7 as VSMCs that displayed high expression levels of MYH11, PLN, CNN1, and RCAN2. Cells in clusters 0, 4, and 5 were classified as qHSCs based on their high expression levels of markers PPARG and NGFR [1]. Finally, cells in cluster 6 were annotated as FBs because of their high expression levels of markers DPT, GGT5, and MFAP4 (Fig. 2B, Supplementary Table 2). Each cluster was assigned a unique name based on the expression of the predominant marker (Supplementary Fig. 3A).

Fig. 2.

Pseudotime analysis revealed the differentiation trajectories of healthy mesenchymal cells. (A) Mesenchymal cells clustered into eight subpopulations in human healthy livers. FB-c6-GGT5, the marker gene of fibroblasts from cluster 6 is GGT5. VSMC, vascular smooth muscle cell; FB, fibroblast; HSC, hepatic stellate cell; tSNE, t-distributed Stochastic Neighbor Embedding. (B) Dot plots showing the expression of mesenchymal cellular marker genes in healthy human livers. (C) A differentiation trajectory of healthy human liver mesenchymal cells inferred by Monocle 2 (pseudotime along differentiation trajectory in inset). (D) Gene expression dynamics along the trajectory of healthy hepatic mesenchymal subtypes. (E) RNA velocity of mesenchymal subtypes in healthy livers, pseudotime (right) detected (purple to yellow). (F) Heatmap showing the TF activity of mesenchymal subtypes in healthy livers. (G) The relative expression of mesenchymal cellular marker genes and transcription factors (TFs) along the differentiation trajectory in healthy livers. (H) Schematic diagram illustrating the molecular mechanisms of mesenchymal cell differentiation in healthy human livers. (I) Dot plot showing the significance and average expression of ligand-receptor interactions between FB and other cell types in healthy human livers.

To determine the evolutionary trajectory of hepatic mesenchymal cells in healthy livers, we used Monocle 2 for trajectory analysis. This trajectory indicated that FBs differentiated into VSMCs and qHSCs (Fig. 2C). To provide further insights into FB differentiation, we identified three sets of differentially expressed genes (DEGs) along the trajectory within branch 1 (Fig. 2D). The first set, consisting of VSMC markers (ACTA2, TAGLN, PLN, and TPM2), increased towards the end of trajectory. Meanwhile, the second set, consisting of HLA-DRA, ANK3 and PCDH9, were highly expressed in FB and showed an upward trend in the early developmental trajectory. The third set of genes (PDE3A, PLA2G5, and CTNNA3) were highly expressed in HSC and showed an upward trend in the late developmental trajectory (Supplementary Fig. 3B).

To validate the developmental trajectory of Monocle 2, we utilized the scVelo pipeline to compute RNA velocity values for each gene in each cell. The analysis predicted that HSCs and VSMCs arose from FBs (Fig. 2E), which is consistent with Monocle2’s results (Fig. 2C). Moreover, VSMC-c2-GBP2, VSMC-c3-SEMA4A, and HSC-c4-IGFBP5 exhibited faster velocities in the late pseudotime than FB-c6-GGT5.

FB differentiation is regulated by a complex network of transcription factors (TFs) that function in conjunction with one another. We used PySCENIC to explore the top five activities in the TF regulatory network. We observed that NFATC2 and NR1H4 displayed high expression and activity levels in the FB regulatory network (Fig. 2F). However, transcriptomics-based TF exploration provides indirect results. As pseudotime increased, the expression of NFATC2, NR1H4, and ZEB2 gradually decreased, along with the expression of FBs-specific genes (IGFBP3, MFAP4, DPT, GGT5, LUM), and the FBs ultimately differentiated into qHSCs and VSMCs (Fig. 2G,H, Supplementary Fig. 3B). During the differentiation of FBs into qHSCs, the qHSC marker genes RGS5 and PLAT was upregulated. Regarding the differentiation of FBs into VSMC, the expression of ACTA2, TPM2, PLN, and MYH11 were upregulated. Finally, CellChat was used to analyse unbiased interactions of mesenchymal subpopulations with ligands and receptors [26]. CellChat revealed that FB significantly interacted with HSC and VSMC, particularly through the ligand-receptor pair GAS6-AXL (Fig. 2I).

In human cirrhotic livers, mesenchymal cells were sub-clustered into 10 subpopulations (Fig. 3A) based on the expression patterns of marker genes (Fig. 3B, Supplementary Table 3). Cells in cluster 0 and 6 express high levels of MYH11, ACTA2, TAGLN, and CNN1, and were identified as VSMCs. Cells in cluster 1, 5, and 7 expressed high levels of activation markers TAGLN, ACTA2, VIM, A2M, CD36, IGFBP5, PDGFRB, and quiescent marker PPARG, showing characteristics common to both qHSC and MFBs [29], and were annotated as HSCs. Cells in cluster 2, 4, 8, and 9 expressed high levels of PDGFRA, VIM, COL3A1, and COL1A1, were identified as MFBs. The cells in cluster 3 were annotated as mesothelial cells expressing high levels of SLPI and KRT19 (Fig. 3B). Each cluster was renamed based on the specific marker gene expressed (Supplementary Fig. 3C).

Fig. 3.

Pseudotime analysis revealed the differentiation trajectories of mesenchymal cells in cirrhotic liver. (A) Human cirrhotic hepatic mesenchymal cells clustered into 10 subpopulations. MFB, myofibroblast. (B) Dot plots showing the expression of mesenchymal cellular marker genes in cirrhotic human livers. (C) A differentiation trajectory of cirrhotic human liver mesenchymal cells inferred by monocle. (D) Gene expression dynamics along the trajectory of cirrhotic liver mesenchymal subtypes. (E) RNA velocity of mesenchymal subtypes in cirrhotic livers. (F) Heatmap showing the TFs activity of mesenchymal cells in cirrhotic livers. (G) The relative expression of mesenchymal cellular marker genes and TFs along the trajectory in cirrhotic livers. (H) Schematic illustrating molecular mechanisms of mesenchymal subpopulations differentiation in cirrhotic human livers. (I) Dot plot showing the significance and average expression of ligand-receptor interactions between HSC and MFB in cirrhotic livers.

Monocle 2 indicated that VSMC-c0-MYH11 and VSMC-c6-ADRA1A were at the beginning of the trajectory, and HSC-c7-PDGFRB was located at the end stage of branch 1. MFB-c8-IGFBP3, MFB-C9-PDGFRA, MFB-c4-CRABP2, and Meso-c3-SLPI were concentrated at the two terminal stages of branch 2 (Fig. 3C). The top 50 genes related to developmental trajectory in branch 1 were clustered into three groups (Fig. 3D). Genes in cluster 1 were highly expressed during the early stage of trajectory development, such as MYH11 and COX4I2, which were highly expressed in VSMCs. Genes in cluster 2 (S100A6, TIMP1, COL1A1, COL3A1, COL1A2, LUM) were highly expressed in MFBs; it has been reported that S100A6 is highly upregulated on activated MFBs [16]. Genes in cluster 3 showed high expression levels in the transitional stage of the developmental trajectory, such as PDGFRB, CD36, MEF2C, and SOX5, which were highly expressed in HSCs (Fig. 3D).

The RNA velocity results are mostly consistent with those of the Monocle, indicating that VSMC-c0-MYH11 and VSMC-c6-ADRA1A were in an early developmental trajectory and tended to differentiate into HSC-c5-RGS5 cells (Fig. 3E). A small proportion of HSC-c1-CD36 and most HSC-C7-PDGFRBs were in the late stages of development, whereas some cells tended to develop into MFB-c2-LUM and MFB-c8-IGFBP3. Surprisingly, MFB-C9-PDGFRA showed a tendency to transform into HSC-C7-PDGFRB (Fig. 3E); we hypothesised that, as liver damage was interrupted, MFBs underwent apoptosis or reverted to the iHSCs state.

To further determine the primary regulators of mesenchymal subpopulations, we evaluated the top six activities of the TF regulatory network using pySCENIC (Fig. 3F). During VSMC differentiation into HSCs, the expression of the VSMCs marker genes MYH11, ACTA2 and TAGLN decreased gradually, whereas the expression of HSC-specific genes (PDGFRB, IGFBP5, CD36, RGS5, A2M) increased gradually, and the expression of HSC TFs (SOX5, MEF2C, TCF7L2) also increased (Fig. 3G,H, Supplementary Fig. 3C). During the transformation of HSCs to MFBs, the expression of MFB-specific genes (COL1A1, COL3A1, COL1A2, LUM, PDGFRA, TIMP1) and TFs (NR1H4, SOX6, PBX1, ZNF23) increased, whereas the expression of PDGFRB, RGS5, JUNB, and MEF2C was gradually downregulated and almost undetectable in MFBs. Furthermore, we investigated the cell-cell communication mechanisms of HSCs and MFBs. The results indicated that all HSC subtypes, except HSC-c7-PDGFRB, had strong interactions with MFBs through the adhesive ligand-receptor pairs MIF-(CD74 + CD44), MDK-LRP1, and NGF-NGFR (Fig. 3I).

Next, we analysed scRNA-seq (GSE137720) data of mouse livers, including two healthy livers and two livers with fibrosis resulting from chronic CCl${}_4$-induced liver injury (Supplementary Fig. 4A). Following quality control (Supplementary Fig. 4B,C), 12538 and 10706 mesenchymal cells were successfully classified into nine clusters (Fig. 4A) and ten clusters (Fig. 4B), respectively. In healthy livers, clusters 0, 1, 3, and 6 expressing high levels of quiescence markers (Ecm1, Reln, Lrat, Hgf, Rgs5, Ngfr) were annotated as qHSCs, whereas clusters 2, 4, and 7 were annotated as FBs expressing high levels of Pi16, Cd34, Mfap4, Col1a1, Col15a1. Cluster 5 and 8 were annotated as VSMCs expressing high levels of Tagln, Acta2, Cnn1, and Myh11 (Fig. 4C, Supplementary Table 4). In the fibrotic livers, cluster 4, 5, and 9 were annotated as VSMCs expressing high levels of Tagln, Acta2, Cnn1, Myh11, and Tpm2. Cluster 0 and 6 with high expression of Pdgfra, Pi16, Fmo2, Col1a1, Col3a1 were annotated as MFBs, whereas cells in cluster 1, 2, 3, 7, and 8 were annotated as aHSCs expressing high levels of Pdgfrb, Hgf, Pth1r, Ecm1, Ngfr, and Adamtsl2 (Fig. 4D, Supplementary Table 5). We renamed the mesenchymal cell subpopulations based on the markers specifically exp-ressed in each cluster (Supplementary Fig. 4D,E).

Fig. 4.

Dynamic evolution trajectory of mouse liver mesenchymal cells is similar to that of humans. (A,B) t-SNE visualisation of mesenchymal cells clustered into nine and 10 subpopulations in mouse healthy and fibrotic livers, respectively. (C,D) Dot plots showing the expression of mesenchymal cellular marker genes in healthy and cirrhotic mouse livers. (E) Trajectory of mouse healthy hepatic mesenchymal cells inferred by Monocle 2. (F,I) Gene expression dynamics along the healthy and fibrotic mouse hepatic mesenchymal cells trajectory. (G,J) RNA velocity of mesenchymal subtypes in healthy and fibrotic livers, pseudotime (right) detected (purple to yellow). (H) A differentiation trajectory of mouse fibrotic liver mesenchymal cells.

In healthy mouse livers, Monocle 2 analysis suggested a differentiation trajectory from FBs to VSMCs and HSCs, which is similar to that in humans (Fig. 2C) but with a simpler trajectory, with each cell type distributed more centrally in three branches (Fig. 4E). Furthermore, this trajectory is consistent with the previously reported results of trajectory analysis of mouse Lin-negative cell in biliary tree fragments (Supplementary Fig. 5A–G, Supplementary Table 6) [14]. Pi16, Col15a1, Cd34, Thy1, Clec3b, Fbln2 are highly expressed in FB (Supplementary Table 7), and Pi16${}^{\mathrm{+}}$ cells serve as resource fibroblasts [30]; CD34 and Thy1 are both stem cell markers [14], further supporting FB as the source cell in the trajectory. Genes related to trajectory were clustered into three clusters. Genes in cluster 1 (Col3a1, Mfap4, Pi16, Clec3b) were highly expressed during the early trajectory. Genes in cluster 2 (Tagln, Tpm1, Myh11, Acta2) were highly expressed during the transitional stage, whereas genes in cluster 3 (Pth1r, Ecm1, Lrat, Ngfr, Hgf, Reln) were highly expressed during the final stage (Fig. 4F). The results of RNA velocity analysis in healthy livers (Fig. 4G) are consistent with the results of Monocle (Fig. 4E). HSCs and VSMCs originated from FBs, which coincided with the results of pseudotime analysis in human healthy livers (Fig. 2C–E).

In mouse fibrotic liver, Monocle trajectory analysis indicated that HSCs and MFBs originated form VSMCs (Fig. 4H), which is consistent with pseudotime analysis results from human hepatic mesenchymal cells with cirrhosis (Fig. 3C,E). The first 50 key genes related to cell trajectory were divided into three clusters. Genes in cluster 1 (Pth1r, Stat1, Vipr1 and Eng) showed high expression mainly in the final stage of the trajectory, genes in cluster 2 (Dpt, Lama2, Cd34, Clec3b, Pi16 and Thy1) were highly expressed in the transitional stage of trajectory, and genes in cluster 3 (Des, Tagln, Rgs4) were highly expressed in the early trajectory (Fig. 4I). The results of RNA velocity analysis in mouse fibrotic liver (Fig. 4J) are consistent with the Monocle results (Fig. 4H). Interestingly, MFB-c0-Pdgfra showed a tendency to differentiated into HSC-c1-Pdgfrb, which is consistent with the RNA velocity of human hepatic mesenchymal cells with cirrhosis (Fig. 3E). LAMA2, CCBE1, PDGFRA and NAALADL2, which were highly expressed in the MFBs of human liver, were also highly expressed in the MFBs of the mouse liver (Supplementary Fig. 6A,B). PDGFRB, ARHGAP42, ASAP1, and GUCY1A2, which were highly expressed in human HSCs, were also highly expressed in mouse HSCs (Supplementary Fig. 6C,D). Therefore, we suggest that MFBs and HSCs in mouse fibrotic liver are similar subpopulations as MFBs and HSCs in human cirrhotic livers.

Finally, we examined the molecular mechanisms underlying mesenchymal cell differentiation in the mouse liver. In healthy mouse livers, FBs TFs Zeb1 showed higher activity during the transformation stage of FBs to VSMCs and qHSCs, the expressions of FBs marker genes (Cd34, Pi16) and TFs (Zeb, Jund) continuously decreased, and the expressions of VSMCs marker genes (Acta2, Myh11, Tagln) and TFs (Mef2c, Klf2) gradually increased (Fig. 5A,B). During FBs differentiation into qHSCs, the expressions of Acta2, Myh11 and TFs Mef2c gradually decreased, whereas the expression of qHSCs marker genes (Ngfr, Lrat and Hgf), qHSC-related TFs (Sox5, Nr1h4, Hand2 and Cebpb) were gradually increased (Fig. 5B,C, Supplementary Fig. 7A). Cell-cell communication analysis showed that FBs interacted with HSCs and VSMCs through the Angptl1-(Itga1 + Itgb1) ligand-receptor pair (Supplementary Fig. 7B).

Fig. 5.

Molecular mechanism of the differentiation of mesenchymal subpopulations in the mouse liver. (A) Heatmap showing the TF activity of mesenchymal cells in healthy mouse livers. (B) Relative expression of mesenchymal cellular marker genes and TFs along differentiation trajectory in healthy mouse livers. (C) Schematic diagram illustrating the molecular mechanism of the differentiation of mesenchymal subpopulations in healthy mouse livers. (D) Heatmap showing the TF activity of mesenchymal cells in mouse fibrotic livers. (E) Relative expression of mesenchymal cellular marker genes and TFs along the trajectory in mouse fibrotic livers. (F) Schematic diagram illustrating the molecular mechanism of differentiation of mesenchymal subpopulations in mouse fibrotic livers.

In mouse chronic fibrosis liver, pySCENIC combined trajectory analysis showed that, at the stage of trans-formation from VSMCs to MFBs and aHSCs, the expression of VSMCs TFs (Mef2c, Thra, Gata4), markers Myh11 and Tagln continuously decreased, the expression of Tagln was the lowest when differentiation to MFBs, and it gradually increased when differentiation into aHSCs (Fig. 5D–F, Supplementary Fig. 7C). In the differentiation from VSMCs to MFBs, the expressions of the MFBs markers Col1a1, Clec3b, and TF Zeb1 was upregulated, and then downregulated with the differentiation of VSMCs to aHSCs. In the differentiation from VSMCs to aHSCs, the expression of aHSCs marker genes Hgf and Pth1r was upregulated, and the expression of the TF Irf1 and Junb gradually decreased from the beginning of VSMCs differentiation to the lowest level at the MFBs stage, and then gradually increased at the differentiation stage of aHSCs (Fig. 5D–F, Supplementary Fig. 7C). Among the ligand-receptor pairs in which VSMCs interacted with MFBs and HSCs, Fgf1-Fgfr2 was the most significant, followed by Ngf-Ngfr; Gas6-Axl also significant (Supplementary Fig. 7D).

To examine the translatability of the core genes in scRNA-seq result, we used TGF-$\beta$1 to treat human hepatic stellate cells (LX-2 cells) to establish a liver fibrosis model (Fig. 6A). In a previous study, we verified the expression of fibrogenic genes (including COL1A1, ACTA2, and FN1) in unstimulated, stimulated and recovered LX-2 cells using reverse transcription polymerase chain reaction (RT-PCR), and successfully constructed a human liver fibrosis model. Following 24 h of starvation, LX-2 cells were stimulated with TGF-$\beta$1 for a further 24 h to induce the transition of qHSCs into MFB. Subsequently, the TGF-$\beta$1 stimulation was removed, and the cells underwent 7-days recovery period, promoting the reversal of MFB into inactive iHSCs.

Fig. 6.

MFB exhibited a reversal to HSC following removal of hepatic injury stimulation. (A) Schematic diagram of the obtention of a human liver fibrosis model using TGF-$\beta$1 to induce LX-2 cells. (B) Quiescent gene (PPARG, NGFR, RGS5, CEBPA, PLAT, RELN, CD36, and IGFBP5) expression in LX-2 (qHSC), TGF-$\beta$1 stimulate stage (MFB) and TGF-$\beta$1 recovery stage (iHSC) (n = 3). Mean $\pm$ SEM; *p 0.05; **p 0.01 (paired t-test). (C) Fibrogenic gene (COL1A1, TIMP1, TGFB1, FN1, PDGFRA, IGFBP3, ACTA2 and PDGFRB) expression in LX-2 (qHSC), TGF-$\beta$1 stimulate stage (MFB) and TGF-$\beta$1 recovery stage (iHSC) (n = 3). Mean $\pm$ SEM; *p 0.05; **p 0.01 (paired t-test).

In LX-2 cells, after exposure to TGF-$\beta$1, qHSC-specific genes (PPARG, NGFR, RGS5, CEBPA, PLAT, RELN, CD36, IGFBP5, GFAP) and TF NR1H4 were mostly downregulated in the fibrotic compared to the quiescent stage (Fig. 6B), and MFB-related genes (COL1A1, TIMP1, TGFB1, FN1, PDGFRA, IGFBP3, ACTA2, PDGFRB, and COL1A2) and TFs (JUNB, ZEB1, PBX1, TCF7L2, and NFATC2) were mostly upregulated (Fig. 6C). This is in agreement with the RNA expression trends from HSC to MFB in pseudotime analysis of human cirrhotic mesenchymal cells (Fig. 3G,H). Among these, ACTA2 expression was upregulated during the initial stimulus stage. After the 7-day recovery period, immune cells or inflammation may still be present in this model, which causes the expression of this gene to continue to be upregulated.

After stimulation removal, qHSC-specific genes were upregulated in the recovery stage, and fibrogenic genes were downregulated (Fig. 6B,C). This implies that MFB transdifferentiates into iHSC, which further supports the results of the reversal of MFB to iHSCs in human cirrhosis and mouse liver fibrosis in the pseudotime analysis (Figs. 3E,4J).