Single‐Cell Transcriptomic Analysis Reveals a Hepatic Stellate Cell–Activation Roadmap and Myofibroblast Origin During Liver Fibrosis in Mice

Abstract

Background and Aims 

HSCs and portal fibroblasts (PFs) are the major sources of collagen‐producing myofibroblasts during liver fibrosis, depending on different etiologies. However, the mechanisms by which their dynamic gene expression directs the transition from the quiescent to the activated state—as well as their contributions to fibrotic myofibroblasts—remain unclear. Here, we analyze the activation of HSCs and PFs in CCL₄‐induced and bile duct ligation–induced fibrosis mouse models, using single‐cell RNA sequencing and lineage tracing.

Approach and Results 

We demonstrate that HSCs, rather than PFs, undergo dramatic transcriptomic changes, with the sequential activation of inflammatory, migrative, and extracellular matrix–producing programs. The data also reveal that HSCs are the exclusive source of myofibroblasts in CCL₄‐treated liver, while PFs are the major source of myofibroblasts in early cholestatic liver fibrosis. Single‐cell and lineage‐tracing analysis also uncovers differential gene‐expression features between HSCs and PFs; for example, nitric oxide receptor soluble guanylate cyclase is exclusively expressed in HSCs, but not in PFs. The soluble guanylate cyclase stimulator Riociguat potently reduced liver fibrosis in CCL₄‐treated livers but showed no therapeutic efficacy in bile duct ligation livers.

Conclusions 

This study provides a transcriptional roadmap for the activation of HSCs during liver fibrosis and yields comprehensive evidence that the differential transcriptomic features of HSCs and PFs, along with their relative contributions to liver fibrosis of different etiologies, should be considered in developing effective antifibrotic therapeutic strategies.

Supported by the National Science Foundation of China (31872797 and 92042304) and the Shanghai Municipal Science and Technology Major Project (2019SHZDZX02).

Potential conflict of interest: Nothing to report.

Liver fibrosis and liver cirrhosis (the advanced stage of liver fibrosis) are leading causes of morbidity and mortality worldwide.^(1,2) Liver fibrosis results from sustained liver injury, which is caused by intrahepatic cholestasis, uncontrolled alcohol consumption, infection with hepatitis viruses, and abnormal metabolic conditions associated with obesity and diabetes. Liver fibrosis is characterized by the emergence of myofibroblasts, which produce excessive extracellular matrix (ECM), resulting in the progressive loss of liver microstructure and metabolic function and eventual liver failure.⁽³⁾ However, effective therapeutics to stop or reverse liver fibrosis are not yet available. The growing number of patients with liver fibrosis highlights the urgent need to develop mechanistic‐based therapies for treating fibrosis.

HSCs, which reside between liver sinusoidal endothelial cells and hepatocytes, are recognized as one of the major origins of myofibroblasts in liver fibrosis of various etiologies.⁽⁴⁾ Quiescent HSCs (qHSCs) store vitamin A–containing lipid droplets and express specific markers, such as Desmin and lecithin retinol acyltransferase (Lrat). In response to liver injury, HSCs gradually lose vitamin A–containing lipid droplets, migrate to the injury sites, and transdifferentiate into myofibroblasts expressing α‐smooth muscle actin (α‐SMA) and collagen.⁽⁴⁾ However, the contribution of portal fibroblasts (PFs) (which reside in the portal triad) in liver fibrosis remains controversial.⁽⁵⁾Lrat‐cre‐based lineage tracing suggests that HSCs give rise to about 90% of myofibroblasts in mouse models of toxic, cholestatic, and fatty liver disease,⁽⁶⁾ while other studies using collagen type I alpha 1 (Col1a1)–enhanced green fluorescent protein (EGFP) (Col‐GFP) mice and multidrug resistance protein 2 (Mdr2) knockout mice have found that both HSCs and PFs contribute to the myofibroblasts reservoir, with more than 70% of myofibroblasts originating from PFs in early cholestatic fibrosis.^(7,8) Previous studies have shown that myofibroblasts originating from HSCs and PFs express different markers,⁽⁷⁾ although it is unclear whether these myofibroblasts of different origins share similar gene signatures or respond similarly to the same antifibrotic therapies.

Single‐cell RNA sequencing (scRNA‐seq) enables genome‐wide gene‐expression analysis at single‐cell resolution, providing unprecedented benefits for identifying cellular heterogeneity, the transition of cellular states, and intercellular communications in complex tissue in health and disease conditions.⁽⁹⁾ Recent studies using scRNA‐seq reveal the emergence of liver parenchymal and nonparenchymal cells during embryonic development,⁽¹⁰⁾ as well as the zonated gene‐expression patterns of hepatocytes, liver sinusoidal endothelial cells (LSECs), and HSCs in adult mouse livers.^(11‐13) Analysis of human and mouse fibrotic livers with scRNA‐seq has uncovered complex intrahepatic communications and the emergence of TREM2+ macrophages.^(11,14) Although different HSC subsets during activation have been identified,^(15‐17) a detailed roadmap of the transcriptional dynamics of the transition from qHSCs to activated collagen‐producing myofibroblasts remains unavailable.

In this study, we perform scRNA‐seq to examine the transcriptional dynamics of HSC activation in hepatotoxin‐induced and cholestatic liver fibrosis mouse models. Our analysis reveals that, during HSC fibrotic activation, the expression of genes associated with the quiescent state is quickly down‐regulated, followed by sequential activation of genes involved in the inflammatory response, cell migration, and eventual ECM production. Furthermore, by combining scRNA‐seq with lineage tracing, we demonstrate that HSCs are the predominant origin for collagen‐producing myofibroblasts in CCL₄‐induced liver fibrosis, whereas PFs are the major source of collagen‐producing myofibroblasts in early bile duct ligation (BDL)–induced cholestatic liver fibrosis. Consequently, Riociguat, which activates soluble guanylate cyclase (sGC) in HSCs, effectively suppresses liver fibrosis in the CCL₄ model, while it shows no obvious therapeutic effect in the BDL model. Overall, our data suggest that different cell origins and the heterogeneity of profibrotic cells should be accounted for in the development of effective anti‐fibrosis treatments.

Materials and Methods

Mice and Genotyping

Guanylate cyclase 1 soluble subunit alpha 1 (Gucy1a1)‐CreER^T2 knock‐in mice were generated by inserting a CreER^T2‐WPRE‐polyA cassette directly behind the start codon of Gucy1a1 (Biocytogen Co., Beijing, China). Gucy1a1‐EGFP transgenic mice were acquired from GENSAT.org. Rosa26‐LSL‐tdTomato (Stock No. 007900) and chemokine (C‐C motif) ligand 2 (CCL2)‐RFP^flox (Stock No. 016849) mice were acquired from the Jackson Laboratory (Bar Harbor, ME). Wild‐type C57BL/6J mice were purchased from Charles River (Beijing, China). Only male mice at the age of 8‐12 weeks were used in this study. All mice were housed in a barrier specific‐pathogen free facility at the Interdisciplinary Research Center on Biology and Chemistry (IRCBC) with free access to food and water. All animal experiments were performed according to the protocol (IRCBC‐2017‐002) that was approved by the Institutional Animal Care and Use Committee of the IRCBC. All genotyping primers are listed in Supporting Table S1.

Liver Injury Models and Riociguat Treatment Regimen

To model hepatotoxin‐induced liver fibrosis, mice were intraperitoneally injected with carbon tetrachloride (CCL₄, 1 mL/kg body weight and 1:5 diluted in corn oil) twice a week for 3 weeks. Livers were harvested 24 hours after the final injection. To model cholestatic liver fibrosis, the common bile duct was ligated under isoflurane anesthesia. The livers were harvested 10 days after the operation. Livers were harvest 24 hours after the second injection. For treatment, mice received Riociguat (Selleck, S8135; 10 mg/kg body weight) or vehicle control twice a day via oral gavage.

Cell Isolation

Liver nonparenchymal cells were isolated following a two‐step perfusion protocol as previously described.⁽¹⁸⁾ The livers were briefly perfused in situ with 37ºC ethylene glycol tetraacetic acid solution, followed by sequential perfusion with Pronase E (1.07433.0005; 0.4 mg/mL; Merck Millipore, Burlington, Ma) for 5 minutes and Collagenase IV (C5138; 0.48 mg/mL; Sigma‐Aldrich, St. Louis, MO) for 7 minutes. Next, the livers were dissected and transferred into 50‐mL falcon tubes and further digested with buffer containing 0.5 mg/mL Pronase E, 0.5 mg/mL Collagenase IV, and 0.02 mg/mL DNase I (10104159001; Roche, Basel, Switzerland) in a 37°C water bath with gentle agitation for 25 minutes. Digested livers were passed through a 70‐μm cell strainer and centrifuged at 580g for 10 minutes. The liver cells were washed twice with Gey’s balanced salt solution (GBSS)/B buffer and resuspended in 32 mL GBSS/B buffer. Finally, the liver cells were subjected to density gradient separation with 9.69% Nycodenz (Histodenz) solution (cat. no. D2158; Sigma‐Aldrich) and centrifuged at 1,380g for 17 minutes at 4ºC without brake to enrich liver nonparenchymal cells. Cell viability (>90%) was determined using trypan blue staining.

Single‐Cell RNA Sequencing and Data Analysis

The liver nonparenchymal cells of three mice per group were pooled and loaded into the 10$\times$ Genomics Chromium Single Cell chips. Libraries were prepared using the Chromium Single Cell 3’ GEM Library & Gel Bead Kit v3 (10$\times$ Genomics, cat. no. PN‐1000075) according to the manufacturer’s instructions, and sequenced on an Illumina Novaseq (San Diego, CA). Sequencing data were processed with Seurat (v3) single‐cell analysis pipeline with modifications. Pseudo‐time trajectory was inferred using Slingshot. (For more details, see Supporting Materials.)

Histology and Immunohistochemistry

Paraffin‐embedded paraformaldehyde (PFA)–fixed livers sections were stained with sirius red. For immunofluorescent staining, all liver samples were fixed with PFA and cryopreserved in optimal cutting temperature compound. The following primary antibodies were used in this study: anti‐Desmin antibody (Ab15200‐1; Abcam, Cambridge, UK), anti‐CD31 antibody (MA3105; Thermo Fisher Scientific, Waltham, MA), anti‐GUCY1A1 antibody (12605‐1‐AP; Proteintech Group, Rosemont, IL), anti‐aSMA antibody (C6198; Sigma‐Aldrich), anti‐THY1 antibody (105302; BioLegend, San Diego, CA), anti‐COLLAGEN 1 antibody (2150‐1410; Bio‐Rad Laboratories, Hercules, CA), anti‐MCP1 antibody (66272‐1‐lg; Proteintech), anti‐Ki67 antibody (ab16667; Abcam), and anti‐CD45 antibody (20103‐1‐AP; Proteintech). (For more details, see Supporting Materials.)

Statistical Analysis

All statistical analyses were performed using GraphPad Prism (v8) or functions in R/Bioconductor packages. Student t test was performed for two‐group comparison, and one‐way ANOVA with Tukey’s Honest Significant Difference post‐hoc test was performed for multiple‐group comparison. Differences between groups were considered to be significant at a P value < 0.05. All data are presented as mean ± SD unless specifically explained.

Results

Single‐Cell RNA‐Sequencing Analysis of Nonparenchymal Hepatic Cells

To understand HSC activation during liver fibrosis, we performed a droplet‐based single‐cell transcriptomic analysis (10$\times$ Chromium). Because liver fibrosis can occur in the pericentral or periportal regions due to different etiologies, we adopted two different experimental models—the intraperitoneal CCL₄ injection–induced hepatotoxic fibrosis model and the BDL‐induced cholestatic fibrosis model—to induce pericentral and periportal injury, respectively (Fig. 1A). We isolated HSCs from the livers of oil‐treated control, CCL₄‐treated, or BDL‐treated mice, using sequential pronase‐collagenase digestion and the Nycodenz gradient separation protocol; we then performed scRNA‐seq.⁽⁶⁾ After quality filtering, 47,752 cells (10,636 cells from oil control, 18,185 cells from CCL₄‐treated, and 18,931 cells from BDL mice), with an average sequencing depth of 2,353 genes per cell, were subjected to further analysis. Unsupervised clustering classified these cells into eight distinct clusters with HSCs as the largest cell cluster (Fig. 1B,C and Supporting Fig. S1A,B). The other cell clusters were classified as endothelial cells, Kupffer cells, PFs, leukocytes, hepatocytes, and cholangiocytes, based on the expression of cell type–specific marker genes (Fig. 1B‐D). Interestingly, we found that cholangiocytes were divided into two clusters—Sox9+ cholangiocytes and Sox9‐ cholangiocytes (Fig. 1B‐D)—which aligns with Tulasi’s finding based on SOX9 (SRY [sex determining region Y]‐box 9) immunofluorescent staining, reflecting the heterogeneity of their intrahepatic localization and regeneration capacity.⁽¹⁹⁾

FIG. 1: ScRNA‐seq of nonparenchymal liver cells. (A) Schematic overview depicting the approach for the experiments. (B) UMAP visualization of isolated hepatic cells, based on 47,752 single‐cell transcriptomes pooled from control (10,636), CCL₄‐treated (18,185), and bile duct‐ligated mice (18,931). Inset: Three‐dimensional UMAP showing the spatial distribution of HSCs and portal fibroblasts. (C) Heatmap depicting the top 25 marker genes (ordered by adjusted P‐values) for each cell cluster. (D) Distinct expression of the cell‐type‐specific genes overlaid on the UMAP of (B). (E) The expression patterns for Gucy1a1 and Gucy1b1. (F) The EGFP expression pattern in the livers of Gucy1a1‐EGFP mice, showing typical HSC morphology. Liver sections were co‐stained with the HSC‐specific marker Desmin and the endothelial cell‐specific marker CD31.

HSCs represent the largest cell clusters, accounting for about 60% of total cells analyzed in this study (Supporting Fig. S1A). Interestingly, we observe a dramatic shift of the HSC clusters from the livers with CCL₄ treatment or BDL toward the PF clusters, while the HSC cluster of the control mouse was distinctly separated from the PFs (Fig. 1B). In contrast, this transcriptomic shift was not observed in PFs, which undergo marked expansion in BDL livers, indicating that HSCs and PFs undergo different activation mechanisms during liver fibrosis.

Single‐cell RNA‐seq provides tremendous power in identifying novel cell type–specific genes.^(9,20) In addition to the well‐known HSC markers Desmin, Lrat, and Pdgfrb (platelet‐derived growth factor receptor beta), our single‐cell analysis revealed specific expression of regulator of G protein signaling (Rgs5), reelin (Reln), parathyroid hormone 1 receptor (Pth1r), transmembrane protein 56 (Tmem56), vasoactive intestinal peptide receptor 1 (Vipr1), angiopoietin like 6 (Angptl6), and ficolin A (Fcna) in HSCs (Supporting Fig. S1B), similar to other recently published single‐cell studies on HSCs.^{(11,14,17,21)} Interestingly, our study found that Gucy1a1 and Gucy1b1, which are the α1 and β1 subunits of the NO receptor sGC, was highly expressed in HSCs (Fig. 1E). To validate the expression specificity of Gucy1a1 and Gucy1b1, we isolated different hepatic cells, including hepatocytes, LSECs, Kupffer cells, and HSCs, and performed bulk RNA‐seq and quantitative PCR analysis. Indeed, the expression of Gucy1a1 and Gucy1b1 was detected only in HSCs, but not in other hepatic cell populations (Supporting Fig. S2A,B). We further confirmed the expression of sGC in the liver using a GUCY1A1 antibody and Gucy1a1‐EGFP reporter mice. Immunofluorescent staining shows that the GUCY1A1 antibody completely colocalized with the EGFP signal, indicating that EGFP expression faithfully reflected the endogenous Gucy1a1 expression (Supporting Fig. S2C). Additional immunostaining and 3D reconstruction revealed that EGFP was colocalized with the HSC‐specific Desmin‐positive cells, and that EGFP+ cells showed typical HSC morphology, with extended cytoplasmic processes wrapping the sinusoidal endothelium (Fig. 1F and Supporting Fig. S2D‐F). Flow cytometry analysis demonstrates that EGFP‐labeled cells were positive for vitamin A–containing lipid droplets, a unique feature of HSCs, further substantiating their HSC identity (Supporting Fig. S2G).

Gene‐Expression Dynamics Define HSC Stages During Fibrotic Activation

To elucidate HSC activation during liver fibrosis, we extracted transcriptomic information on all HSCs for further analysis. Uniform Manifold Approximation and Projection (UMAP) visualization showed that HSCs could be divided into two distinct clusters: qHSCs and active HSCs (aHSCs) (Fig. 2A). In accordance with the established dogma that aHSCs up‐regulate migration and ECM‐associated genes during liver fibrosis,⁽⁴⁾ the aHSCs from both the CCL₄‐treated and BDL livers showed increased actin alpha 2, smooth muscle (Acta2), transgelin (Tagln), Col1a1, Col1a2, and Timp1 expression (Fig. 2B,C). On the contrary, a group of genes including Lrat and Rgs5, Ecm1, Angptl6, Vipr1, Gucy1a1, and Gucy1b1—which were highly expressed in qHSCs—were markedly down‐regulated in aHSCs (Fig. 2B,C); this indicates that a dramatic transcriptional reprogramming occurs during HSC activation.

FIG. 2: Transcriptomic roadmap of HSC activation. (A) UMAP visualization of the qHSCs and aHSCs from pooled HSCs of control, CCL₄‐treated, and bile duct‐ligated mice. (B) The expression of the qHSC marker genes Lrat and Ecm1 and aHSC marker genes Col1a1 and Acta2 overlaid on the UMAP of (A). (C) Heatmap depicting the top 50 representative genes (ordered by adjusted P‐values) of qHSC and aHSC. (D) Pseudotime trajectory indicating the activation of HSC during liver fibrosis. (E) Heatmap representing four gene clusters based on the dynamic gene expression pattern over pseudotime from the qHSC state to the aHSC state. Cluster 1: 202 genes; cluster 2: 82 genes; cluster 3: 288 genes, and cluster 4: 401 genes. Plot showing average expression of the four gene clusters along pseudotime. GO biological pathway analysis of the four gene clusters. Expression kinetics of representative genes of four gene clusters.

To understand the transcriptional dynamics of HSCs during activation, we reconstructed the gene‐expression dynamics of HSC activation using Slingshot; we identified four distinct gene‐expression modules along a single pseudo‐time trajectory (Fig. 2D,E). Genes of the first module were highly expressed in qHSCs, with their expression decreasing immediately following HSC activation. Interestingly, we found that the first module contained well‐known HSC markers, such as Lrat1 and Rgs5, as well as genes that are associated with vascular tone and relaxation, such as Gucy1a1, Gucy1b1, and Vipr1. This indicates that HSCs not only quickly lose their quiescent identity but also lose their vascular relaxation capacity, contributing to the increase in intrahepatic vascular resistance and the development of portal hypertension during liver fibrosis (cluster 1 of Fig. 2E). The genes of the second module were transiently up‐regulated at the early stages but decreased rapidly afterward. Most genes of the second module encode for inflammatory cytokines, such as Ccl2 and chemokine (C‐X‐C motif) ligand 10 (Cxcl10), which are involved in cytokine–chemokine–mediated signaling pathways and leukocyte recruitment (cluster 2 of Fig. 2E). The third module contains genes such as Acta2, Tagln, and other cell migration and contractility‐associated genes; their expression levels slowly increased, reaching their peaks at midstage, followed by a mild decrease at the end stage during liver fibrosis (cluster 3 of Fig. 2E). The fourth module, which primarily encompasses genes associated with ECM components such as Col1a1 and Col1a2, exhibited delayed up‐regulation and reached their peaks at the end stage during HSC activation (cluster 4 of Fig. 2E).

To analyze whether aHSCs could be further divided into different subclusters, we integrated the differential gene‐expression patterns and the relative spatial distribution of the HSCs along the pseudo‐time trajectory (Supporting Fig. S3). Based on the different gene‐expression kinetics of the four gene clusters, we divided aHSCs into three distinct subclusters, defined as stage‐1 aHSCs, stage‐2 aHSCs, and stage‐3 aHSCs (Fig. 3A). In line with the four dynamic gene‐expression modules, a heatmap of the highly expressed genes of each HSC subcluster—as well as Gene Expression Omnibus analysis—demonstrated that stage‐1 aHSCs expressed high levels of inflammatory cytokines, including Ccl2, Cxcl10, and Ccl7. Stage‐2 aHSCs demonstrated the expression of cell mobility and contractility‐associated genes such as Acta2, Tpm1, vimentin (Vim), Tagln, and tenascin C (Tnc). Stage‐3 aHSCs expressed extremely high levels of ECM deposition and organization‐related genes, including Col1a1, Lox, and Lum (Fig. 3B,C). Because activated HSCs undergo proliferation to expand the collagen‐producing cell reservoir during liver fibrosis, we analyzed aHSC proliferation based on the expression of proliferation markers (marker of proliferation Ki‐67 [Mki67] and minichromosome maintenance complex component 6 [Mcm6]). Overlay of the Mki67‐expressing and Mcm6‐expressing cells with an HSC UMAP indicated that only a small fraction of aHSCs underwent proliferation, with most of them being stage‐2 aHSCs (Fig. 3D).

FIG. 3: Subclusters of activated HSCs. (A) UMAP displaying the qHSC and the three subclusters of aHSCs (stage‐1 aHSCs, stage‐2 aHSCs, stage‐3 aHSCs). (B) Heatmap depicting the top 50 differentially expressed genes for qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, and stage‐3 aHSCs (sort by average logFC in descending order). (C) GO biological analysis of the top 50 differentially expressed genes of (B). (D) UMAP visualization of the expression pattern of the cell proliferation markers Mki67 and Mcm6. (E) UMAP displaying the qHSC and the three subclusters of aHSCs in the individual livers of control, CCL₄, and BDL mice (left panel). Proportion of qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, and stage‐3 aHSCs in the individual livers of control, CCL₄, and BDL mice (right panel).

We further compared the difference in HSC activation in fibrotic livers induced by either CCL₄ or BDL. UMAP projection of HSCs from individual livers indicates that aHSCs of the CCL₄‐treated liver showed a different activation pattern than aHSCs from BDL liver. Most of the aHSCs in CCL₄‐treated liver were stage‐2 aHSCs and stage‐3 aHSCs, accounting for 38.1% and 43.2% of total HSCs, respectively. On the contrary, most of the aHSCs in BDL liver were stage‐1 aHSCs (44.5%), whereas stage‐2 and stage‐3 aHSCs only accounted for 11.3% and 11.0%, respectively, of total HSCs in BDL liver (Fig. 3E). This result suggests that, unlike the widespread activation of HSCs in the CCL₄ model, HSCs are less activated in the BDL liver at day 10, probably due to the confined portal injury in BDL livers.

To systematically delineate the gene regulatory networks that control HSC homeostasis and activation, we subjected the transcriptomes of the HSCs to unbiased single‐cell regulatory networks inference and clustering (SCENIC) analysis.^(22,23) We identified 284 orthologous transcription factors grouped into four regulatory modules, indicating differential transcription factor activity during HSC activation. Consistent with previous studies based on HSC‐specific gene deletion,^(24,25) SCENIC analysis showed that transcription factors of the Ets family, Gata family, interferon regulatory factor 1 (Irf1), Irf2, forkhead box o 1 (Foxo1), and peroxisome proliferator‐activated receptor (Ppar) family were involved in maintaining the HSC quiescent status. Notably, as with the gene‐expression patterns in HSCs at different stages (Supporting Fig. S4A,B), Nfkb1/2, Relb, and signal transducer and activator of transcription 2 (Stat2/3) were involved in controlling inflammatory gene expression in HSCs at stage 1; the serum response factor (SRF), Jun, and Fos transcription factors were responsible for the enhanced cell migration phenotype, whereas Wilms tumor protein 1 (Wt1) and runt‐related transcription factor 1 (Runx1) participated in increased collagen deposition.⁽²⁶⁾ Interestingly, SCENIC analysis also identified transcription factors not previously implicated in controlling HSC function, such as hes related family bHLH transcription factor with YRPW motif 2 (Hey2) and transformation‐related protein 53 (Trp53)—which may play roles in maintaining HSC quiescence—as well as Activating Transcription Factor 3/6 (Atf3/6), and meis homeobox 1/2 (Meis1/2), which may be involved in regulating HSC activation.

Morphological Dynamics of HSCs During Fibrotic Activation

It has been shown that, along with the dramatic changes in gene expression during in vitro activation, HSCs lose their unique stellate shape and acquire a distinct spindle‐like myofibroblast morphology.⁽²⁷⁾ However, it is unclear whether this morphological change is associated with HSC activation in vivo. To investigate, we used Gucy1a1‐EGFP mice to track the change in HSC morphology. Gucy1a1‐EGFP mice were treated with CCL₄ or BDL to induce liver fibrosis. Next, the whole liver sections of control, CCL₄‐treated, and BDL‐treated Gucy1a1‐EGFP mice were scanned in z‐stack mode using spinning disk confocal microscopy. The morphology of EGFP‐positive HSCs was first 3D‐reconstructed, followed by skeletonization, to show the cell body and cytoplasmic processes. The analysis demonstrated that HSCs of CCL₄‐treated livers showed a marked reduction in total filament volume, branch numbers per cell, total branch length per cell, and mean branch length per cell, compared with HSCs from the control liver (Fig. 4A,B). On the contrary, HSCs from BDL livers showed a reduced branch number per cell compared with HSCs from control livers; their total filament volume, branch number per cell, total branch length, and mean branch length per cell were still significantly higher than in HSCs from CCL₄‐treated livers. This result indicates that HSCs are more activated in CCL₄‐treated liver than in BDL‐treated liver, which aligns with the results of the transcriptomic‐based single‐cell analysis, showing that most aHSCs in CCL₄‐treated livers were at stage 2 and stage 3, while most aHSCs in BDL liver were at stage 1 (Fig. 3E).

FIG. 4: Morphological analysis of HSCs in liver fibrosis. (A) The liver sections of sham‐operated, CCL₄‐treated, or bile duct‐ligated Gucy1a1‐EGFP mice were stained with collagen‐1 and CD31 antibodies. The morphology of EGFP+ HSCs were skeletonized; cell bodies are in blue and cellular processes are in red. (B) Morphological analysis of EGFP+ HSCs in the liver sections of sham‐operated, CCL₄‐treated, or bile duct‐ligated Gucy1a1‐EGFP mice. Statistical analyses were performed using one‐way ANOVA with Tukey HSD test. ns: not significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. (C) Liver sections of CCL₄‐treated Gucy1a1‐EGFP mice were stained with collagen1 and CD31. Stage‐1, stage‐2, and stage‐3 HSCs were defined based on their relative distance to the collagen‐1‐positive stage‐3 aHSCs. Zoomed‐in image showing the morphology of qHSCs and aHSCs at stage 1, stage 2, and stage 3. Arrow indicate the skeletonized morphology of individual aHSCs. (D) Morphological comparison of qHSCs and HSCs at stage 1, stage 2, and stage 3. Statistical analyses were performed using one‐way ANOVA with Tukey HSD test. ns: not significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

We next examined whether aHSCs at different activated stages have distinct morphological features. To this aim, we characterized the morphology of qHSCs from the livers of control Gucy1a1‐EGFP mice to serve as a reference (Fig. 4C, first and second rows). To identify aHSCs at different stages, we first defined the terminally activated collagen‐producing stage‐3 aHSCs in the liver sections of CCL₄‐treated Gucy1a1‐EGFP mice, based on their collagen‐1 expression. Next, we defined stage‐1 and stage‐2 aHSCs based on their relative distance to stage‐3 aHSCs (Fig. 4C, third through sixth rows). While the results show that the stage‐1 aHSCs still maintained a stellate morphology similar to qHSCs, stage‐2 and stage‐3 aHSCs acquired a spindle‐like morphology with reduced cell volume and branch counts (Fig. 4C,D). Nevertheless, HSCs (even those adjacent to the collagen‐positive fibrotic scar) in the BDL livers still maintained a stellate morphology, indicating that most HSCs were not terminally transdifferentiated into myofibroblasts in the livers 10 days after BDL (Supporting Fig. S5).

Portal Fibroblasts—Not HSCs—Are the Major Source of Collagen‐Producing Cells in Cholestatic Liver Fibrosis

Both HSCs and PFs are activated during cholestatic liver fibrosis⁽⁵⁾; however, their contributions to collagen‐producing myofibroblasts in cholestatic liver fibrosis remain controversial.^(6‐8) In this study, we found that the low percentage of collagen‐producing stage‐3 aHSCs in the BDL liver (Fig. 3E) could not explain the massive collagen deposition in the portal area (Supporting Fig. S6), suggesting that the collagen in the BDL liver may be deposited primarily by PFs. Therefore, we compared the single‐cell transcriptomes of PFs and HSCs, finding that the gene‐expression patterns of PFs closely resembled those of stage‐3 aHSCs (Fig. 5A). Nevertheless, PFs could be easily distinguished from stage‐3 aHSCs in 2D and 3D UMAP presentations, due to the expression of PF‐specific makers, including Thy1, fibulin‐1 (Fbln1), elastin (Eln), dermatopontin (Dpt), microfibril‐associated glycoprotein 4 (Mfap4), and growth arrest‐specific 6 (Gas6) (Fig. 1B,C and Fig. 5B,C). Furthermore, the PFs of CCL₄‐treated liver expressed similar levels of Col1a1 and Col1a2 compared with aHSCs, whereas the PFs of BDL liver expressed significantly higher levels of Col1a1 and Col1a2 than aHSCs (Fig. 5C), indicating that PFs are the principal collagen‐producing cells in cholestatic liver fibrosis. As in previous studies,⁽²⁸⁾ the percentage of PFs was very low in the healthy liver (Fig. 5D). The number of PFs was increased in the CCL₄‐treated liver. However, the number of PFs was dramatically elevated in the BDL liver–a 57‐fold increase compared with the number of PFs in healthy liver and a 12‐fold increase compared with the number of PFs in CCL₄‐treated livers (Fig. 5D‐F).

FIG. 5: PFs are the major source of collagen‐producing myofibroblasts in BDL liver. (A) Heatmap depicting the top 50 differentially expressed genes in the qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, stage‐3 aHSCs, and PFs of CCL₄‐treated and BDL livers (sort by average logFC in descending order). The right panel shows the representative genes. (B) UMAP visualization of qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, stage‐3 aHSCs, and portal fibroblasts pooled from sham, CCL₄, and BDL mice. The right panel shows the 3‐D UMAP representation of the five cell clusters. (C) Violin plot showing the expression levels of Thy1, Fbln1, Mfap4, Eln, Dpt, Gas6, Col1a1, and Col1a2 in aHSCs and PFs in the livers of CCL₄ or BDL mice. (D) UMAP visualization of qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, stage‐3 aHSCs, and portal fibroblasts in the individual livers of control, CCL₄, and BDL mice. (E) Proportion of qHSCs, stage‐1 aHSCs, stage‐2 aHSCs, stage‐3 aHSCs, and portal fibroblasts in the individual livers of control, CCL₄, and BDL mice. (F) Normalized percentage of PFs in sham, CCL₄, and BDL livers.

We further verified the contribution of HSCs and PFs to collagen‐producing myofibroblasts in CCL₄‐induced and BDL‐induced liver fibrosis, using Gucy1a1‐EGFP mice. Confocal imaging analysis demonstrates that, in the CCL₄‐treated liver, most of the collagen‐producing cells were derived from HSCs, as they were EGFP‐positive. On the contrary, the collagen‐producing cells in the BDL liver were negative for EGFP but positive for the PF marker Thy1 (Fig. 6A). Notably, our analysis showed that the expression of Gucy1a1 and Gucy1b1 declined during HSC activation (Fig. 2F). Thus, to exclude the possibility that a dramatic reduction in EGFP expression occurred in terminally activated HSCs—and to precisely determine the contribution of HSCs to the collagen‐producing myofibroblasts in cholestatic liver fibrosis—we generated Gucy1a1‐CreER^T2 mice to permanently label HSCs during liver fibrosis (Supporting Fig. S7A). We crossed Gucy1a1‐EGFP with Gucy1a1‐CreERT2 and Rosa26‐LSL‐tdTomato mice; we injected five doses of Tamoxifen at the age of 4 weeks and determined the labeling efficiency at the age of 8 weeks. Immunofluorescent imaging and flow cytometry analysis revealed that approximately 90% of EGFP+ HSCs were labeled with permanently expressed tdTomato following Tamoxifen induction (Supporting Fig. S7B,C). Next, we performed BDL to induce cholestatic liver fibrosis in Gucy1a1‐EGFP::Gucy1a1‐CreERT2::Rosa26‐LSL‐tdTomato mice at the age of 8 weeks. Lineage tracing results showed that the collagen‐positive cells in the BDL liver were predominantly positive for Thy1 but negative for both EGFP and tdTomato (Fig. 6B,C). Therefore, the lineage tracing data and scRNA‐seq analysis demonstrate that, during cholestatic liver injury, PFs (which are activated and undergo dramatic expansion) are the major source of collagen deposition, whereas HSCs make a minimal contribution to the excessive collagen deposition.

FIG. 6: Lineage tracing reveals HSCs contribute to the collagen‐producing myofibroblasts in CCL₄ liver instead of BDL liver. (A) Liver sections of control, CCL₄, or BDL‐treated Gucy1a1‐EGFP mice were stained with Collagen‐1 and Thy1. (B) Liver sections from bile duct‐ligated Gucy1a1‐EGFP::Gucy1a1‐CreERT2::ROSA26‐LSL‐tdTomato mice were stained with Thy1. (C) Liver sections from bile duct‐ligated Gucy1a1‐EGFP::Gucy1a1‐CreERT2::ROSA26‐LSL‐tdTomato mice were stained with collagen‐1.

Activating sGC Inhibits HSC Activation and Ameliorates Hepatotoxin‐Induced Liver Fibrosis, but Not Cholestatic Liver Fibrosis

Several studies have shown that sCG‐activating reagents inhibit liver fibrosis in different preclinical rodent models of pig serum, CCL₄ treatment, and high‐fat diet–induced NASH, or BDL.^(29‐32) In this study, we found that expression of sGC subunits Gucy1a1 and Gucy1b1 were decreased following HSC activation, with the lowest expression levels in stage‐3 aHSCs (Fig. 2F and Fig. 7A,B). More importantly, our data show that collagen‐producing cells in BDL model originated primarily from PFs, which did not express sGC (Fig. 6A,B). These results prompted us to re‐evaluate the therapeutic effects of sGC activation in CCL₄‐induced and BDL‐induced liver fibrosis models. Considering that later‐stage aHSCs may impede their response to sGC‐activating reagents due to the reduction of sGC expression, Riociguat, an sGC stimulator that has been Food and Drug Administration–approved for treating pulmonary arterial hypertension,⁽³³⁾ was administered 1 day after BDL or the first dose of CCL₄ (Fig. 7C). The analysis shows that Riociguat treatment significantly reduced collagen deposition in CCL₄‐treated livers, as evidenced by sirius red staining (Fig. 7C,D). The reduction in collagen deposition following Riociguat treatment was further confirmed by a remarkable decrease in the hydroxyproline content in the liver (Fig. 7D). Furthermore, Riociguat treatment also reduced alanine aminotransferase and aspartate aminotransferase levels in the serum, indicating that activating sGC protected the liver from CCL₄‐induced liver injury (Fig. 7D). On the contrary, Riociguat did not reduce liver fibrosis and liver injury in BDL mice at day 10 (Fig. 7E‐F). These data substantiating the notion that sGC‐activating agents should only be used to combat liver fibrosis that aHSCs, which express sGC, are the predominant source of collagen deposition, but not cholestatic liver fibrosis, in which PFs are the major source of accumulated collagen.

FIG. 7: Differential responses of Riociguat treatment in CCL₄‐ and BDL‐induced liver fibrosis. (A) The correlation of the expression of Gucy1a1 and Gucy1b1 in HSCs and PFs. (B) The expression of Gucy1a1 and Gucy1b1 in individual HSCs and portal fibroblasts were overlaid on the UMAP. Dot lines outline the different subclusters according to Figure 4B. (C) Schematic overview of the experimental design. CCL₄ was injected twice a week for 3 weeks and Riociguat (10 mg/kg body weight) was administered twice a day by oral gavage. Representative images of liver sections from control, CCL₄‐, and CCL₄+ Riociguat‐treated mice were stained with sirius red to show the collagen deposit. (D) Sirius red positive area, hydroxyproline concentration, and serum ALT and AST enzymatic activity in the livers of control, CCL₄‐, and CCL₄+ Riociguat‐treated mice were quantitated (Corn oil + vehicle, n = 4; CCL₄ + vehicle, n = 6; CCL₄ + Riociguat, n = 6). Statistical analyses were performed using one‐way ANOVA with Tukey HSD test. ns: not significant; *P < 0.05; **P < 0.01; ****P < 0.0001. (E) Schematic overview of the experimental design of the BDL experiment. Vehicle or Riociguat (10 mg/kg body weight) was administered twice a day for 10 days. Representative images of liver sections of vehicle‐ or Riociguat‐treated BDL mice were stained with sirius red. (F) Quantification of sirius red positive area, hydroxyproline concentration, ALT and AST enzymatic activity, and total bilirubin (T‐BIL) concentration in the serum of vehicle‐ or Riociguat‐treated BDL mice (sham + vehicle, n = 4; BDL + vehicle, n = 6; BDL + Riociguat, n = 5). Statistical analyses were performed using one‐way ANOVA with Tukey HSD test. ns: not significant; *P < 0.05; **P < 0.01; ****P < 0.0001.

To gain further insight into the mechanisms by which sGC activation inhibits HSCs activation, we stimulated isolated mouse HSCs with Riociguat and performed bulk transcriptomic analysis, which provides higher sequencing depth than scRNA‐seq. Riociguat treatment significantly elevated cyclic guanosine monophosphate concentration in cultured HSCs (Supporting Fig. S8A). Transcriptomic profiling reveals that Riociguat treatment markedly suppressed the expression of 822 genes in HSCs (Supporting Fig. S8B). Gene‐set enrichment analysis indicates that Riociguat treatment potently inhibited gene expression related to cell migration and actin cytoskeleton rearrangement, inflammatory cytokines response, and collagen production (Supporting Fig. S8C‐F). Interestingly, we also identified a small set of genes belonging to the Hippo pathway–including Ankyrin repeat domain 1 (Ankrd1), Cysteine‐rich angiogenic inducer 61 (Cyr61), connective tissue growth factor (Ctgf), and Thrombospondin 1 (Thbs1)–which were also significantly suppressed by Riociguat treatment (Supporting Fig. S8G,H). This agrees with a previous finding that inhibition of YAP impedes liver fibrosis.⁽³⁴⁾ These results suggest that activating sGC could potently inhibit the expression of inflammatory cytokines, as well as HSC migration, thereby blocking HSC activation at early stages and preventing HSCs from differentiating into later‐stage collagen‐producing myofibroblasts.

Discussion

Chronic fibrosis caused by various etiologies, including sustained hepatitis virus infection and alcohol‐associated and nonalcoholic liver diseases, is a global health care burden, with total mortality of two million deaths per year.⁽¹⁾ Many strategies targeting TGFβ, PDGF, CTGF, liver X receptor (LXR), and NOX have been developed to inhibit the differentiation and proliferation of collagen‐producing myofibroblasts. However, their clinical deployment has achieved limited success, and no effective antifibrosis treatment options are yet available. Hence, it is crucial to improve our understanding of the molecular mechanisms underlying liver fibrosis and to translate that knowledge into better mechanistic‐based therapeutic strategies. In this study, we describe the HSC activation program and evaluate the origins of collagen‐producing myofibroblasts in CCL₄‐induced and BDL‐induced fibrotic livers, using single‐cell transcriptomic analysis and lineage‐tracing strategies. We also uncover that myofibroblasts of different origins respond differently to antifibrotic therapy.

Cellular heterogeneity within individual cell types—including hepatocytes, LSECs, and HSCs—in healthy livers has recently been investigated using single‐cell analysis. Similar to hepatocytes, HSCs also demonstrate zonated gene‐expression patterns, with periportal HSCs expressing nerve growth factor receptor (Ngfr) and insulin like growth factor binding protein 3 (Igfbp3) and pericentral HSCs expressing Adamts like 2 (Adamtsl2) and lysyl oxidase like 1 (Loxl1).⁽²¹⁾ Single‐cell analysis of HSCs from various fibrosis mouse models reveals that qHSCs and aHSCs displayed distinct gene‐expression profiles. In line with previous findings using bulk‐RNA‐seq, single‐cell analysis data from our study and others consistently show that aHSCs profoundly up‐regulate migration and ECM‐associated genes such as Col1a1 and Acta2.^(14,21) Furthermore, heterogeneity within aHSCs in fibrotic liver has also been recently reported. Rosenthal et al. show that aHSCs in NASH livers can be further divided into a proliferating cluster, an intermediate activated cluster, an immune and inflammatory cluster, and a classic fibrogenic myofibroblast cluster.⁽¹⁷⁾

Active HSCs are known to be inflammatory, chemotactic, migrative, contractile, and characterized by increased collagen production.⁽⁴⁾ However, it is unknown whether HSCs acquire these profibrotic features in a simultaneous or sequential manner during activation. Using scRNA‐seq and pseudo‐time trajectory, we infer a roadmap of cell‐state transition from qHSCs to terminally differentiated collagen‐producing aHSCs during liver fibrosis. Following fibrotic insult, qHSCs sequentially differentiate into stage‐1 aHSCs with an inflammatory gene signature. Subsequently, they progress into stage‐2 aHSCs with a migrative and contractile gene signature, and eventually into ECM‐producing stage‐3 aHSCs. Thus, the analysis indicates that HSC activation follows a sequential activation model orchestrated by waves of gene expression; this requires inactivation of the gene signature of the current stage before the initiation of the gene signature of the next stage. Interestingly, this sequential activation has also been observed in HSCs undergoing in vitro activation.⁽¹⁵⁾ Another surprising finding is that—in contrast to HSCs, which undergo the dramatic transcriptomic shift during activation—PFs undergo dramatic expansion, especially in the BDL livers, as evidenced by the single‐cell and histological analysis, although they did not exhibit transcriptomic changes during liver fibrosis.

Although HSCs from CCL₄‐treated and BDL‐treated livers follow the same activation trajectory, the percentage of aHSCs and the compositions of different subclusters in CCL₄‐treated and BDL‐treated livers are quite different. Our scRNA‐seq reveals that HSCs in CCL₄‐treated liver were highly activated, with the elevated ECM produced by stage‐3 aHSCs. In contrast, the HSCs in BDL liver are less activated, as evidenced by the fact that the dominant subcluster was stage‐1 aHSCs. Meanwhile, the number of PFs was dramatically expanded in BDL liver. We further generated Gucy1a1‐CreER^T2 mice to trace the origin of myofibroblasts in fibrotic livers. Analysis of the livers of Gucy1a1‐CreER^T2::Rosa26‐LSL‐tdTomato mice shows that the predominant myofibroblasts in BDL livers were tdTomato‐negative but Thy1‐positive, whereas the collagen‐producing myofibroblasts in CCL₄‐treated liver were tdTomato‐positive, suggesting that the dominant origin of myofibroblasts in cholestatic liver injury is PFs. Our observations align with those of previous studies that used Col‐GFP and Mdr2 knockout mice^(7,8); our findings contradict those of a previous study based on Lrat‐Cre, which found that myofibroblasts in BDL liver have an HSC origin.⁽⁶⁾ The discrepancy possibly stems from the fact that, in addition to HSCs, the constitutive Lrat‐Cre may also label some PFs during embryonic development, as they share the same mesenchymal origin. Notably, our analyses were performed 10 days after BDL, which only reflects the early or intermediate stage of cholestatic injury. The contribution of HSCs to myofibroblasts in cholestatic liver injury would increase in late‐stage cholestatic models, as Iwaisako et al. demonstrated that the number of aHSCs and the expression of Col1a1 in livers are higher 20 days after BDL than 5 days after BDL.⁽⁷⁾

Single‐cell RNA sequencing shows remarkable power in discovering new cell type–specific makers. In this study, unsupervised single‐cell analysis showed that Lrat, Vipr1, Ecm1, Gucy1a1, and Gucy1b1 were highly and specifically expressed in HSCs. In particular, their expression levels are associated with the maintenance of the quiescent state in HSCs. Among these genes are Gucy1a1 and Gucy1b1, which are the α1 and β1 subunits of the NO receptor sGC.⁽³³⁾ The scRNA‐seq data and Gucy1a1‐CreER^T2::Rosa26‐LSL‐tdTomato lineage‐tracing show that sGC (Gucy1a1/Gucy1b1) is specifically expressed in HSCs and in HSC‐derived myofibroblasts, but not in Thy1‐positive PFs and PF‐derived myofibroblasts. Furthermore, a steady decrease in Gucy1a1 and Gucy1b1 expression in aHSC was observed, which would lead to impaired HSC relaxation and the onset of portal hypertension. Preclinical evidence shows that the activation of sGC can inhibit liver fibrosis in mouse models induced with pig serum, CCL₄, or diet.^(29,31,32) Another study reports that the sGC stimulator Riociguat reduces portal pressure and liver fibrosis in both CCL₄ and BDL rat models.⁽³⁰⁾ However, in our study, Riociguat treatment potently inhibited liver fibrosis in CCL₄‐treated mice, but not in the BDL mouse model. The discrepancy may stem from the different experimental designs, as we performed BDL in mice and analyzed the liver 10 days after operation, whereas Schwabl et al. performed BDL in rats and analyzed the livers 3 weeks and 5 weeks after operation, as the previous study by Iwaisako et al. demonstrates that the contribution of HSCs to collagen‐producing myofibroblasts increases over time after BDL.⁽⁷⁾ Therefore, it is reasonable to expect that a substantial proportion of sGC‐expressing HSCs may differentiate into stage‐3 HSCs in cholestatic livers 3 or 5 weeks after BDL; thus, the administration of the sGC stimulator Riociguat under such conditions may exhibit antifibrosis effect. Nonetheless, our study suggests that the differences between collagen‐producing myofibroblasts of different origins in fibrotic livers should be further investigated and taken into account for designing better antifibrotic therapy strategies in the future.

In fibrotic liver, aHSCs are extensively heterogeneous, with terminally activated aHSCs located at the fibrotic center, where the liver injury happens, and less‐activated aHSCs located away from the fibrotic center. This poses a significant challenge for effective antifibrotic treatment, because most therapies under development are single‐agent treatments that target only one subcluster of aHSCs; therefore, they are not able to effectively stop the progression of liver fibrosis. Currently, many agents targeting specific gene and fibrotic signaling pathways have been developed and tested in various clinical trials. The CCL2‐inhibitor Bindarit⁽³⁵⁾ and the CCR2‐antagonist RS‐504393⁽³⁶⁾ block CCL2 secreted by stage‐1 aHSCs; the Rock‐inhibitor Y‐27632⁽³⁷⁾ and the FAK‐inhibitors PF‐562271⁽³⁸⁾ impede the migration of stage‐2 aHSCs; and LOXL2‐neutralizing antibody interferes in the organization of ECM secreted by stage‐3 aHSCs.⁽³⁹⁾ However, the clinical deployment of these agents as single‐agent therapies has had limited success. Therefore, a combination of agents to target multiple subclusters of aHSCs may provide better clinical outcomes in treating liver fibrosis.

In summary, by using scRNA‐seq with different preclinical fibrosis mouse models, we provide a detailed transcriptional roadmap of HSC activation that is shared by both CCL₄‐induced and BDL‐induced liver fibrosis. Following fibrotic activation, HSCs sequentially acquire inflammatory, migrative, and contractile phenotypes, before eventually differentiating into collagen‐producing myofibroblasts. We also demonstrate that HSCs minimally contribute to collagen‐producing fibroblasts in the early phase of BDL‐induced liver fibrosis, which leads to the differential therapeutic outcome of HSC‐targeting Riociguat treatment. Our discovery provides important insight into the mechanisms of HSC activation, heterogeneity, and differential cellular‐source myofibroblasts during liver fibrosis. Our findings potentially open avenues for developing better therapeutic strategies for the treatment of liver fibrosis.

Author Contributions

W.Y. prepared the fibrosis mouse models and performed therapeutic treatment. W.Y., K.N., and N.L. performed the 10$\times$ single‐cell library preparation. H.H., W.Y., T.W., F.Z., L.W., L.L., W.W., and J.H. analyzed the single‐cell data. H.H. and W.Y. analyzed the bulk RNA‐seq data. W.Y., J.Z., N.S., K.J., X.Y., C.Z., and X.Y. contributed to tissue collection and histology analysis. W.Y., H.H., and J.H. designed the study and wrote the manuscript.

Data and Materials Availability

All sequencing data are deposited in the Gene Expression Omnibus (GEO) with the accession number GSE171904 (single‐cell transcriptomic profiling) and GSE171885 (bulk transcriptomic analysis of Riociguat‐treated HSCs).