Single-cell Transcriptome Analysis Reveals the Potential Role of Hepatic Stellate Cells in Liver Fibrosis

Abstract

Upon activation, hepatic stellate cells (HSCs) can convert into fibroblasts and increase the production of extracellular matrix, a major cause of liver fibrosis (LF) and a growing health issue worldwide. Other mechanisms by which HSCs may induce fibrosis remain to be explored, and the role of cell dynamic gene expression in liver fibrogenesis is not well understood. In this study, analysis by single-cell transcriptome sequencing (scRNA-seq) was used to explore the potential effects of HSCs in a bile duct ligation (BDL)-induced mouse model of LF, followed by the identification of novel targets for clinical diagnosis. Methods: Liver tissue collected from BDL and sham-operated C57BL/6J mice was used for scRNA-seq. To systematically dissect the molecular and cellular events following fibrosis, the scRNA-seq data was analyzed for differential gene expression, KEGG, pseudotime trajectory, and cellular communication. Morphological changes in the BDL and sham livers were examined by hematoxylin and eosin (H&E) staining, Masson’s trichrome staining, fiber staining, and Sirius red staining. Results: The scRNA-seq analysis performed on the BDL and sham groups revealed the gene expression of 20,764 cells across 27 cell types. Antioxidant levels declined markedly in HSCs from BDL mice, leading to a more pronounced occurrence of ferroptosis. We also found evidence suggesting that elevated apelin signaling and platelet activation in HSCs contributed to the increased synthesis of extracellular matrix and collagen fibers. The large accumulation of immune cells in the liver of BDL mice induces different outcomes for HSCs. Conclusion: The results of this study provide further insight into the cellular and molecular alterations that occur within a specific subset of HSCs during LF, offering valuable information on potential targets for therapeutic intervention.

Graphical abstract

Keywords

scRNA-seq / liver fibrosis / HSCs / immune cells / ferroptosis / platelet activation

1. Introduction

Reversal of the early, pathological stage of liver fibrosis (LF) offers the greatest potential for restoring liver function and overall hepatic health [1]. Currently, there are no established markers for the diagnosis and prevention of LF. Given the increasing number of patients affected by this disease, there is an urgent need for new diagnostic and therapeutic targets to prevent and treat fibrosis [2, 3].

The activation of hepatic stellate cells (HSCs) results in excessive deposition of extracellular matrix (ECM), which represents the classical mechanism underlying LF. There is increasing evidence that, in addition to HSCs, other cell types within the liver are also crucial in the pathogenesis and progression of this disease [4, 5, 6, 7]. These cell types promote LF by modulating HSCs through intercellular signaling [8]. Currently, research on liver macrophages and myeloid-derived suppressor cells (MDSCs) focuses on two aspects. Firstly, cytokines or chemokines secreted by these cells affect the proliferation or differentiation of HSCs. Secondly, following liver injury, the aggregation and infiltration of cells promotes injury repair and fibrosis, thus aggravating the progression of LF. Bioinformatics analysis combined with in vivo mouse experiments have shown that liver macrophages secrete Tumor Necrosis Factor-\alpha (TNF-\alpha) and Interleukin-1\beta (IL-1\beta), which activate HSCs through the nuclear factor kappa-B (NF-\kappaB) signaling pathway [9]. Matrix metalloproteinase 8 (MMP8) belongs to the matrix metalloproteinase family and has been shown to activate HSCs via the Extracellular-regulated kinase 1/2 (ERK1/2) pathway by regulating collagen1 (Col1a1) expression, as well as the production of \alpha-smooth muscle actin (\alpha-SMA) when secreted by macrophages [10]. The blocking of Transforming growth factor-\beta (TGF-\beta) signaling induces HSCs to downregulate the expression of related genes in MDSCs, thereby slowing the progression of LF [11]. Additionally, MDSC-related cytokines can suppress \alpha-SMA expression and induce apoptosis in HSCs [7, 12].

A better understanding of liver structure and of the function of each liver cell type is essential for elucidating the pathogenesis of liver-related diseases and identifying valuable therapeutic targets. Sequencing technology has garnered increasing attention in genomic and transcriptomic studies due to its advantages of high throughput, short cycle, low cost, and efficient cell capture [13, 14]. The application of this technology can lead to a better understanding of genetic development and of the molecular mechanisms underlying the heterogeneity between healthy individuals and those with LF.

This study utilized single-cell transcriptome sequencing (scRNA-seq) to investigate the specific cell types and cell type-specific differential gene expression patterns associated with LF in a bile duct ligation (BDL)-induced mouse model of LF. The BDL mice showed increased antioxidant levels that were associated with a significant reduction in HSC-specific ferroptosis. Furthermore, increased levels of apelin signaling and platelet activation were observed in HSCs, potentially contributing to increased synthesis of ECM and collagen fiber. These findings provide valuable insights into potential targets for the treatment of LF.

2. Materials and Methods

2.1 Animals and Treatment

Male C57BL/6J mice (8 to 10 weeks of age) were housed in specific pathogen free (SPF) cages in approved facilities under a 12 h light/dark cycle and controlled temperature (22 °C) and humidity (40% to 60%). The mice were randomly divided into sham and BDL groups. BDL and sham operation in mice were performed, under anesthesia with isoflurane. Inhale oxygen containing 2% isoflurane (at a flow rate of 1 L/min) to maintain the anesthesia of mice. For sham operations, all steps were the same, except for common bile duct ligation. Mice were euthanized with 2% pentobarbital sodium (150 mg/kg) intraperitoneal injection at day 21 after BDL for collecting tissue samples, including blood and livers. All experimental procedures were approved by the Medical Ethics Committee of the Shanxi Institute of Traditional Chinese Medicine.

2.2 Measurement of ALT and AST Concentrations

Ocular blood from mice was collected into prepared tubes and an automatic biochemical analyzer (AU480, BECKMAN COULTER, Brea, CA, USA) was employed to examine liver function markers, including the levels of alanine aminotransferase (ALT) and aspartate aminotransferase (AST).

2.3 Histological Analysis

Mouse livers collected 4 weeks after BDL were fixed in cold 4% formalin, dehydrated, embedded in paraffin, and then cut into 4 µm-thick sections. The sections were rehydrated and stained with hematoxylin and eosin (H&E), Sirius red, reticular fiber staining, or Masson trichrome. The degree of fibrosis was assessed histologically using light microscopy by an observer who was blinded to the sample identity.

2.4 Preparation of Liver Cell Suspension

Liver cells were extracted using a two-step perfusion approach. First, the liver was sequentially perfused with Pronase E (0.4 mg/mL; Merck Millipore, Burlington, MA, USA) and Collagenase IV (0.48 mg/mL; Sigma-Aldrich, St. Louis, MO, USA) for 5 minutes and 7 minutes, respectively. Next, the livers were dissected, placed into 50 mL falcon tubes, and digested in a 37 °C water bath with moderate agitation for 25 minutes using a buffer containing 0.5 mg/mL Pronase E, 0.5 mg/mL Collagenase IV, and 0.02 mg/mL DNase I (Roche, Basel, Switzerland). The livers were then passed through a cell strainer with a pore size of 70 µm and centrifuged at 580 g for 10 minutes. Liver cells were resuspended in 32 mL of Gey’s balanced salt solution (GBSS)/B buffer after being washed twice with GBSS/B buffer. To enrich liver nonparenchymal cells, the liver cells were subjected to density gradient separation using 9.69% Nycodenz (Histodenz) solution (Sigma-Aldrich) and centrifuged at 1380 g for 17 minutes at 4 °C. The viability of cells (>90%) was assessed using trypan blue staining.

The above procedure was performed to obtain as many cell types as possible from the mouse livers. To analyze cell characteristics during the process of cholestatic LF, nuclear suspensions were prepared from mouse liver tissues and the nuclei labeled on the machine. Specific surface markers were used to screen for certain types of cells within the liver. Additionally, single-cell transcriptomics was performed on the 10\times Genomics platform, which is based on polydT capture of mRNA. This experimental process excludes mycoplasma, since prokaryote mRNA does not have a polyA tail. The genome alignment rates for the three samples in our project were all >90%, indicating the data were of mouse origin and there was no contamination from other species.

2.5 Single Cell RNA-seq and Data Preprocessing

scRNA-seq was performed at the Stanford Functional Genomics Facility (SFGF) after Fluorescence-activated cell sorting (FACS) sorting and using the 10\times Chromium 3^′ v2 kit (10\times Genomics, Pleasanton, CA, USA) according to the manufacturer’s protocol. The target number of captured cells ranged from 2000 to 10,000. R1 (read 1 primer sequence) was introduced to the molecules during Gel Bead-in-emulsion (GEM) incubation. During construction of the library, P5, P7, a sample index, and R2 (read 2 primer sequence) were inserted using End Repair, A-tailing, Adaptor Ligation, and Polymerase Chain Reaction (PCR). The P5 and P7 primers used in Illumina bridge amplification were included in the final libraries. Sequencing was performed on a HiSeq 4000 (Illumina, Inc., San Diego, CA, USA) at a median depth of 115,216 reads/cell. Raw sequencing data were processed using the Cell Ranger pipeline (version 2.1.0, 10\times Genomics) and mapped to the hg19 mouse reference genome to generate a unique molecular identifier (UMI) count matrix by cell barcoding. For UMI counting, reads that were uniquely mapped to the transcriptome and intersected an exon at least 50% of the time were included.

2.6 Quantification of Gene Expression

UMI sequences were corrected for sequencing mistakes prior to quantification, while valid barcodes were recognized using the EmptyDrops approach. UMI counting and cell barcode calling were used to generate the cell by gene matrices. Each sample’s cell by gene matrices were independently imported using the R package Seurat (version 3.1.1, https://satijalab.org/seurat/) for downstream analysis. Cells with an extremely high number of UMIs (30,000) or mitochondrial gene content (10%) were excluded. Cells with 500 genes or >4000 genes were also removed.

To reduce the impacts of batch effect and behavioral circumstances on clustering, samples were aggregated by employing canonical correlation analysis and mutual nearest neighbor analysis. A variance stabilizing transformation was applied to identify 2000 highly variable genes in each sample. Anchors were established between data points, and correction vectors were generated to create an integrated expression matrix that was subsequently used for clustering. The matrix was scaled and dimensionally reduced using principal component analysis (PCA). Following the jackstraw process, a resampling test was conducted by permuting a portion of the data (1% by default) and rerunning PCA. This process created a ‘null distribution’ of gene score. Significant principal components (PCs) were defined as those genes exhibiting a high enrichment with low p-value for downstream clustering and dimensional reduction. The Harmony method was employed for batch correction of single-cell transcriptome clustering.

2.7 Cell Clustering Analysis

The R package Seurat was used to perform Uniform Manifold Approximation and Projection (UMAP). Seurat employs a graph-based clustering strategy, where cell distances are determined using previously identified PCs. Cells are embedded in a graph of shared k-nearest neighbors (SKNN), where edges are formed between cells with similar gene expression patterns. To partition this graph into strongly connected quasi-cliques or communities, the SKNN graph was first built in PCA space using Euclidean distance. Edge weights between any two cells were then refined based on the shared overlap with their immediate neighbors, as measured by the Jaccard distance. We then applied modularity optimization techniques (Smart Local Moving algorithm) to iteratively combine cells, thereby optimizing the conventional modularity function.

2.8 Cell Type Annotation

Log-normalized matrices were processed using the SingleR R package to annotate cell types [15]. This was achieved by comparing the gene expression profiles of reference cell types to those of single cells. First, a Spearman coefficient was computed for the expression of each individual cell in the reference dataset. Subsequently, multiple correlation coefficients were aggregated per cell type to yield a single value representing each cell type for each individual cell. Following this, the correlation analysis was performed once more with SingeR, focusing solely on the top cell types identified in the previous step. The investigation specifically targeted highly variable genes. Finally, each single cell was assigned the cell type corresponding to the highest value calculated in the preceding analysis.

2.9 Analysis of Cell-type Specific Differential Gene Expression

The Wilcoxon rank sum test was utilized to compare the expression levels of each gene within a specific cluster against those of the remaining cells. A list of criteria was employed to identify significant genes, while minimizing false positive rates. Firstly, genes within the target cluster needed to show at least a 1.28-fold increase or decrease in expression. Secondly, genes had to be expressed in at least 25% of the cells in the target cluster, with a p-value of 0.05. To evaluate inter-group differences, the Model-based Analysis of Single-cell Transcriptomics (MAST) method was utilized to test the significance of differentially expressed genes. The top 20 genes were further selected as marker genes according to the results of differential expression analysis per cluster. The expression patterns of each marker gene were then visualized using heatmap, UMAP, and bubble diagram.

2.10 KEGG Enrichment Analysis

For the gene set generated from the above differential gene expression analysis, clusterProfiler R tools were utilized to conduct pathway enrichment analyses based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/). Statistical analysis was performed using Fisher’s exact test, followed by the Benjamini-Hochberg correction, with an adjusted p-value of 0.05 to determine the significance.

3. Results

3.1 Changes in Liver Tissue and Functional Parameters in Mice With BDL-induced Liver Fibrosis

A BDL-induced LF model of C57BL/6J mice was used to comprehensively investigate the cellular and genetic changes associated with LF by comparison with sham-operated mice. In the BDL group, the proximal bile duct was thickened, the gallbladder gradually expanded, and the liver became smaller, with a harder texture and increased surface granularity (Fig. 1A). The levels of ALT and AST were significantly higher in the BDL group compared to the sham-operated group, peaking 2 weeks after surgery. They subsequently began decreasing 3–4 weeks after surgery, but remained significantly higher in BDL mice (p 0.01) (Fig. 1B).

Histological analysis of liver tissue from BDL mice revealed severe focal necrosis, vacuolated hepatocytes, increased and enlarged reticular fibers, and elevated levels of collagen fibers (Fig. 1C–F). In contrast, the liver tissue morphology, lobular structure, central venules, and bile duct structure of sham mice appeared normal, with no evidence of fibrosis, cellular degeneration or inflammatory changes.

3.2 Single-cell Transcriptomic Analysis of Parenchymal and Nonparenchymal Hepatic Cells

Single-cell transcriptomic analysis (10\times Chromium) was performed to objectively study the mechanism involving HSCs during LF. A total of 20,764 cells (10,789 in the BDL group and 9975 in the sham group) were analyzed following data quality control at the gene and cell levels. Unsupervised clustering of these cells was performed using the canonical correlation analysis function, which generated 27 clusters with different transcript levels (Fig. 2A). A heatmap was used to identify cell type-specific genes in HSCs, hepatocytes and macrophages that showed significant differential expression between the sham and BDL groups (Fig. 2B). Of note, the expression levels of collectin sub-family member 10 (Colec10) and Neurexin 1 (Nrxn1) were higher in HSCs compared to other cell types, while the expression levels of Colec10 and Nrxn1 were much higher in the BDL group than in the sham group (Fig. 2C).

We next investigated genes that were positively regulated and showed increased expression in the BDL group compared to the sham-operated group. The functional implications of up-regulated genes in different cell types was explored, followed by KEGG enrichment analysis of up-regulated genes that were specific to HSCs, hepatocytes, and macrophages (Fig. 2D). Genes associated with focal adhesion, platelet activation, ECM-receptor interaction, and regulation of the actin cytoskeleton were up-regulated in HSCs. Genes related to lysosomes, endocytosis, platelet activation and regulation of the actin cytoskeleton showed higher expression levels in hepatocytes, while genes enriched in pathways related to peroxisome, complement and coagulation pathways, as well as carbon metabolism, were up-regulated in macrophages. Overall, both HSCs and hepatocytes demonstrated significant activation of platelet-related processes and synthesis of the actin cytoskeleton in BDL mice, with HSCs in particular showing pronounced accumulation of ECM.

3.3 Subtype Identification of HSCs From BDL and Sham Mice

HSCs from the BDL and sham mice were divided into 10 distinct subgroups based on their gene expression patterns and distribution (Fig. 3A,B). The distribution of HSCs in sham-operated and BDL mice was significantly clustered, with the majority located in the cluster 2 subgroup (Fig. 3C,D). All HSC subgroups expressed cell type-specific genes (e.g., Ecm1, Colec11, Hgf, Vipr1) (Fig. 3E–R). Interestingly, genes related to high levels of fibrosis, including Col1a1, Col1a2 and Acta2, were specifically identified in clusters 0, 5, and 8 (Supplementary Fig. 1s1–s3). The HSCs were then divided into different sample subgroups that were either unique to the sham group, or shared by both groups (Fig. 3S). Differential analysis identified molecular markers that were shared (Lyz2, Cd74, Scube1, C1qb, Krt8), or sham-specific (Ttc39c, Dpyd, Ppara, Hao1, Egfr) (Fig. 3T). According to KEGG analysis, highly expressed genes in the shared subgroup were related to ribosomes, oxidative phosphorylation, PPAR signaling, and ferroptosis (Fig. 3U), while those in the sham-specific subgroup were related to bile secretion and the peroxisome (Fig. 3V).

A trajectory study of HSCs in both the BDL and sham groups showed distinct expression patterns between resting and activated HSCs (Supplementary Fig. 1s4–s9). In the sham group, HSCs were predominately in a resting state, whereas HSCs in the BDL group mainly showed higher expression levels of pro-fibrosis genes. ScRNA-seq analysis consistently demonstrated that HSCs from BDL mice expressed high levels of migration- and ECM-related genes, including Col1a1 and Acta2. ScRNA-seq and temporal trajectory analysis revealed that resting HSCs, which did not express key pro-fibrosis genes, were only present in sham mice. In BDL mice, HSCs were activated and there was an increased level of ECM1 expression. As the liver damage progressed, these cells transformed into terminal HSCs with high Acta2 and Col1a1 expression.

3.4 The Interaction Between HSCs and Hepatocytes Modulates the Expression of Ferroptosis-related Genes, Thereby Facilitating Liver Fibrosis

Among the co-enrichment results for HSCs and hepatocytes, the ferroptosis signaling pathway stood out in both (Fig. 4A). A role for ferroptosis in LF has been supported by several studies [16, 17, 18]. To identify genes associated with the ferroptosis pathway in BDL mice, the expression levels of key ferroptosis marker genes in hepatocytes and HSCs were comprehensively analyzed. Ferroptosis pathway-related inhibitory genes, such as Ascl1, Glutathione peroxidase 4 (Gpx4), Cp, Fth1, Ftl1, Gclm, Trf, were all up-regulated in the HSCs of BDL mice, while the Ascl1 ferroptosis-specific regulatory gene was down-regulated (Fig. 4B,C). The expression patterns of these genes were identical in both HSCs and hepatocytes (Fig. 4B,C). Gpx4 and Gclm may affect HSC and hepatocyte ferroptosis via synthesis of glutathione (GSH), while Cp and Trf could reduce free Fe by binding Fe ions, thereby inhibiting Fe-dependent cell death. The distribution of Fth1 and Acsl1 expression in BDL and sham mice indicated that Ascl1 was differentially expressed mainly in hepatocytes and HSCs, while the expression of Fth1 varied across multiple cell types (Fig. 4D).

3.5 Changes in Apelin Signaling-related Gene Expression and the Cell Cycle in BDL Mice

Apelin signaling-related genes were significantly upregulated in HSCs and hepatocytes from BDL mice (Fig. 5A,B). Superoxide dismutase (SOD1/SOD2) and Gpx1 also showed elevated expression levels in BDL HSCs (Fig. 5C), indicating these cells have elevated antioxidant capacity that can potentially render them less susceptible to ferroptosis. Conversely, SOD1/SOD2 were down-regulated in BDL hepatocytes, while NOX4/MAPK3 were up-regulated (Fig. 5C). Apelin has been implicated in the proliferation of bile duct cells through Nox4/ROS/ERK-dependent signaling, potentially inducing HSC activation through intracellular ROS pathways. Apelin expression might increase antioxidant capacity by up-regulating the expression of SOD1/SOD2 and consequently reducing ferroptosis in HSCs, while concurrently promoting NOX4/ROS/ERK signaling in hepatocytes. To investigate the association between Apelin signaling and the cell cycle, the cell cycle status of both BDL and sham groups was evaluated. A considerable reduction in the number of non-cycling cells was observed in BDL mice, accompanied by an increase in cells within the G2 and M phases (Fig. 5D). Further examination of cell cycle stage-specific genes revealed significant up-regulation in the expression of G2- and M-phase-related genes in BDL mice, suggesting potential acceleration of the cell cycle by Apelin (Fig. 5E).

3.6 Platelet Activation and Changes in Downstream Gene Expression in HSCs

Up-regulation of the receptor for platelet-derived growth factor (PDGF) is an important metabolic manifestation observed in HSCs during LF [19]. Notably, genes in HSCs and hepatocytes that show significant differential expression between the BDL and sham mice were enriched in pathways related to platelet activation (Fig. 2D), suggesting dysregulation of platelet-related gene expression in BDL-induced LF. The expression profile in HSCs and hepatocytes of genes related to platelet activation was compared between BDL and sham mice. The majority of these genes were up-regulated in the HSCs and hepatocytes of BDL mice (Fig. 6A). Concurrently, PDGF expression was elevated in both HSCs and hepatocytes (Fig. 6B), with significant up-regulation of the PDGF receptor observed specifically in the HSCs of BDL mice (Fig. 6C). Additional investigation into platelet-specific downstream factors revealed significant up-regulation of TGF-\beta1, the TGF-\beta1 receptor, and SMAD in HSCs (Fig. 6D). Platelets are likely to initiate activation of the TGF-\beta1 signaling pathway by phosphorylating SMAD, thereby stimulating trans-differentiation of HSCs into myofibroblasts.

3.7 Accumulation of Immune Cells in the Liver of BDL Mice Induces Different Outcomes for HSCs

To assess the level of immune cell infiltration in BDL, common immune cell markers were used to annotate the immune cell types in 27 cell clusters (Fig. 7A). Eosinophils, red pulp macrophages, and T memory cells were found to be the most important immune cells in the sham group (Fig. 7A), and T memory cells, red pulp macrophages, neutrophils and MDSC in BDL mice (Fig. 7A). Red pulp macrophages are mainly involved in immune surveillance, iron metabolism and cell clearance [20]. The interaction between macrophages and HSCs plays a crucial role in regulating the formation and excessive deposition of ECM in the liver [4, 21]. Specifically, macrophages secrete prostaglandins that bind to receptors on the surface of HSCs, leading to the activation of autophagic signaling in these cells [22]. Various studies have demonstrated that macrophages play a crucial role in promoting the proliferation and differentiation of HSCs into fibroblasts through the release of TGF-\beta, thereby facilitating the development of LF [23, 24]. However, when the proliferation rate of macrophages decreases, the ECM is degraded, which may play a role in reversing LF [25]. Given that macrophages have both pro- and anti-fibrotic effects, we compared the expression of M1 and M2 macrophage marker genes in macrophages from BDL mice and sham mice. Interestingly, M1 macrophage marker genes (e.g., Il1b, Tlr8, Cd80 and Cd86) were down-regulated in BDL mice, whereas M2 macrophage marker genes (e.g., pparg, Arg-1, F13a1 and Slc40a1) were up-regulated in the sham group (Fig. 7B). This suggests that macrophages in BDL mice may undergo a polarity shift from M1 to M2, thereby altering their regulatory role in the liver and further promoting LF.

Other annotated immune cells such as MDSCs were also notable (Fig. 2A). These cells promote normal immune responses in various organs and tissues. During liver injury, MDSCs accumulate in the liver and induce the proliferation of HSCs and promote fibrosis. This has become a new hotspot area for research into LF. In the present study, a significant increase in MDSCs was observed in BDL mice, in parallel with a decrease in eosinophils (Fig. 7C). In contrast, sham mice had very few MDSCs and a high level of eosinophils (Fig. 7D). MDSC marker genes (e.g., Ly6g, Itgam, Fut4, Cd33, Ccr2) were also expressed at much higher levels in BDL mice compared to sham mice, with the expression mainly concentrated in MDSCs and neutrophils (Fig. 7E).

3.8 High SCD1 Expression in HSCs may Influence the Development of Liver Fibrosis Through the Wnt Signaling Pathway

Stearoyl-CoA desaturase 1 (SCD1) is a pivotal enzyme in hepatic lipid anabolism, catalyzing the crucial step of converting saturated fatty acids into monounsaturated fatty acids. The overexpression of SCD1 in activated HSCs is associated with diet-induced steatohepatitis and fibrosis through the regulation of Wnt signaling. The present study found that SCD1 was significantly upregulated in the majority of HSCs in BDL mice, whereas only a small fraction of HSCs in the sham group expressed SCD1 (Fig. 8A,B). Furthermore, BDL-induced HSCs and hepatocytes showed increased expression of many of the Wnt signaling-related genes (Fig. 8C). Additionally, \beta-catenin, a key marker of Wnt signaling activation, was prominently expressed in the majority of HSCs from BDL mice, but showed only limited expression in HSCs from sham mice (Fig. 8D). In summary, elevated expression of SCD1 in BDL mice triggers abnormal activation of the canonical Wnt signaling pathway, likely contributing to the activation of HSCs and the development of LF.

4. Discussion

This study characterized the transcriptomic profile of hepatocytes in the livers of sham-operated and BDL mice, identifying important genes and associated mechanisms that underlie liver fibrogenesis across various cell populations. A total of 27 distinct cell lineages were identified, of which HSCs emerged as one of the predominant cell populations undergoing modulation during LF. Subsequent subgroup and trajectory assessments of HSCs from the BDL mouse model and sham-operated mice revealed distinct gene expression patterns in resting and activated cells. While predominantly quiescent HSCs were observed in the sham group, HSCs derived from BDL mice showed increased expression of profibrotic genes such as Col1a1 and Acta2. Colec10 was specifically expressed in BDL HSCs, consistent with findings from other studies [26]. Furthermore, Nrxn1, a member of the neuroin family associated with the upregulation of vascular smooth muscle contraction [27, 28], showed specific expression patterns in BDL-derived HSCs. KEGG analysis of HSCs revealed the enrichment of pathways related to platelet activation signaling, Ras signaling, and vascular smooth muscle contraction.

Ten subsets of HSCs were identified through the analysis of scRNA-seq clustering. These were present in both the sham and BDL groups, as well as in exclusive subsets unique to the sham group. Previous study has classified HSCs into two subsets, namely central venous-associated HSCs (CaHSC) and portal venous-associated HSCs (PaHSC), considered to be the primary source of differentiation into collagen-producing myofibroblasts [29]. The subset of HSCs in the portal vein region referred to as the 1-HSC region appears not to be associated with LF [30]. In the current study, high expression of HSC subset 2 was observed in the sham group, leading us to hypothesize this subset may represent PaHSC. Unfortunately, this result was not verified by subsequent experiments. Notably, Krt8 (keratin 8) was universally expressed in HSC subsets, albeit at lower levels in healthy mice from the sham group. Remarkably, these findings are consistent with a previous study indicating that Krt8 deficiency confers protection against pulmonary fibrosis and suppresses pro-fibrotic cell activation by modulating oxidative phosphorylation [31]. This was also supported by our analysis of differential gene expression and by the KEGG analysis. In addition, Cd74 is a marker of microglial activation, while TGF-\beta, a cytokine in LF, inhibits Cd74 [32]. Consequently, reduced expression of TGF-\beta promotes downregulation of the Cd74 receptor on microglia, thereby attenuating glial reactivation [15].

Macrophages play a critical role in the reversal of mild LF. The differential regulation of M1 macrophages, which are associated with pro-inflammatory activity, and M2 macrophages, which are associated with anti-inflammatory effects, can lead to the progression or reversal of LF [33, 34]. During active hepatitis, M2 macrophages can also induce fibrosis by promoting fibrinolysis after the inflammation resolves, with the anti-inflammatory effect becoming predominant [35]. Therefore, mice with LF may remain in the stage of inflammation progression, during which time the dominant factors for reversal of LF have not yet emerged.

The KEGG analysis of annotated genes in HSCs and hepatocytes revealed involvement of the ferroptosis, Apelin, and SCD1-related signaling pathways in LF. Abnormal iron metabolism leads to the accumulation of reactive oxygen species (ROS), ultimately inducing unprogrammed cell apoptosis. Lipid peroxide (LPO) plays an important role in hepatitis C virus (HCV) replication, inducing high levels of ferritin that lead to excessive iron deposition, a direct manifestation of ferroptosis [36]. Gpx4 regulates ferroptosis by converting toxic LPO into non-toxic lipoalcohols [37]. Inhibition of HSC ferroptosis improves LF in a ccl4-induced mouse model of LF, highlighting the close relationship between ferroptosis and LF [38]. The research have suggested that Apelin promotes cell proliferation of the HSC line LX2 through activation of ERK signaling during the induction of LF [39]. Moreover, MAPK can regulate Apelin signaling to attenuate LF injury in BDL mice [40]. AMPK-SCD1 inhibits ferroptosis during liver cancer, and the down-regulation of SCD1 expression significantly increases ferroptosis [41]. Hence, it is evident that various factors such as ferroptosis, oxidative stress, Apelin signaling, and SCD1-related signaling are interconnected within the liver, providing novel insights into the pathogenesis of LF.

In this study, we identified some specific genes related to liver fibrosis, such as Colec10 and Nrxn1, which might serve as non-invasive serum or imaging-based biomarkers for the early detection of fibrosis. At the same time, single-cell RNA sequencing of patient liver biopsy samples can further identify subtype-specific characteristics to stratify the severity of the disease. Furthermore, in future research, we will continue to conduct in-depth analysis of the results from single-cell sequencing, mainly focusing on the communication analysis between immune cells and hepatic stellate cells. We will also validate the observed interactions between immune cells and hepatic stellate cells through in vitro co-culture experiments. For instance, sorting macrophages or bone marrow-derived suppressive cells from mice with bile duct ligation and co-culturing them with primary hepatic stellate cells will enable direct assessment of paracrine signals (such as the secretion of transforming growth factor-\beta or platelet-derived growth factor) and their effects on the activation markers of hepatic stellate cells. Additionally, neutralizing antibodies or small molecule inhibitors targeting key cytokines (such as anti-transforming growth factor-\beta and platelet-derived growth factor receptor inhibitors) can be used to analyze the specific pathways mediating the activation of hepatic stellate cells. This research has identified operable targets in hematopoietic stem cell biology, immune regulation, and metabolic pathways, opening up new avenues for precision medicine in liver fibrosis. The application of these research results in clinical practice may be able to meet the needs of non-invasive diagnosis and improvement of disease treatment, thereby enhancing the treatment efficacy of patients.

5. Conclusion

Targeted therapies for LF should address the underlying causes, as well as inhibiting the deposition of ECM. However, current treatments are still unable to reduce the high mortality rates associated with liver disease. scRNA-seq can identify individual cells, accurately analyze the cell composition of a sample, and reveal the gene expression patterns of individual cells during LF. This method can also categorize cells and help to identify new cell types. The present study has revealed changes in BDL-induced fibrotic cells in the liver, cell type-specific gene expression, the immune landscape, and the regulation of signaling pathways including ferroptosis, thus providing novel insight into the pathogenesis of LF (Fig. 9).

Availability of Data and Materials

The data of the present study are available from the first or corresponding author upon reasonable request and with prior permission.

References

[1] Tsuchida T, Friedman SL. Mechanisms of hepatic stellate cell activation. Nature Reviews. Gastroenterology & Hepatology. 2017; 14: 397–411. https://doi.org/10.1038/nrgastro.2017.38.

[2] Parola M, Pinzani M. Liver fibrosis: Pathophysiology, pathogenetic targets and clinical issues. Molecular Aspects of Medicine. 2019; 65: 37–55. https://doi.org/10.1016/j.mam.2018.09.002.

[3] Parola M, Pinzani M. Liver fibrosis in NAFLD/NASH: from pathophysiology towards diagnostic and therapeutic strategies. Molecular Aspects of Medicine. 2024; 95: 101231. https://doi.org/10.1016/j.mam.2023.101231.

[4] Cheng S, Zou Y, Zhang M, Bai S, Tao K, Wu J, et al. Single-cell RNA sequencing reveals the heterogeneity and intercellular communication of hepatic stellate cells and macrophages during liver fibrosis. MedComm. 2023; 4: e378. https://doi.org/10.1002/mco2.378.

[5] Tsuchiya Y, Seki T, Kobayashi K, Komazawa-Sakon S, Shichino S, Nishina T, et al. Fibroblast growth factor 18 stimulates the proliferation of hepatic stellate cells, thereby inducing liver fibrosis. Nature Communications. 2023; 14: 6304. https://doi.org/10.1038/s41467-023-42058-z.

[6] Cheng D, Chai J, Wang H, Fu L, Peng S, Ni X. Hepatic macrophages: Key players in the development and progression of liver fibrosis. Liver International: Official Journal of the International Association for the Study of the Liver. 2021; 41: 2279–2294. https://doi.org/10.1111/liv.14940.

[7] Zhang C, Sui Y, Liu S, Yang M. The Roles of Myeloid-Derived Suppressor Cells in Liver Disease. Biomedicines. 2024; 12: 299. https://doi.org/10.3390/biomedicines12020299.

[8] Wesley BT, Ross ADB, Muraro D, Miao Z, Saxton S, Tomaz RA, et al. Single-cell atlas of human liver development reveals pathways directing hepatic cell fates. Nature Cell Biology. 2022; 24: 1487–1498. https://doi.org/10.1038/s41556-022-00989-7.

[9] Pradere JP, Kluwe J, De Minicis S, Jiao JJ, Gwak GY, Dapito DH, et al. Hepatic macrophages but not dendritic cells contribute to liver fibrosis by promoting the survival of activated hepatic stellate cells in mice. Hepatology (Baltimore, Md.). 2013; 58: 1461–1473. https://doi.org/10.1002/hep.26429.

[10] Naim A, Baig MS. Matrix metalloproteinase-8 (MMP-8) regulates the activation of hepatic stellate cells (HSCs) through the ERK-mediated pathway. Molecular and Cellular Biochemistry. 2020; 467: 107–116. https://doi.org/10.1007/s11010-020-03705-x.

[11] Liu X, Liu M, Wu H, Tang W, Yang W, Chan TTH, et al. PPP1R15A-expressing monocytic MDSCs promote immunosuppressive liver microenvironment in fibrosis-associated hepatocellular carcinoma. JHEP Reports: Innovation in Hepatology. 2024; 6: 101087. https://doi.org/10.1016/j.jhepr.2024.101087.

[12] Gao M, Huang A, Sun Z, Sun Y, Chang B, Zhang JY, et al. Granulocytic myeloid-derived suppressor cell population increases with the severity of alcoholic liver disease. Journal of Cellular and Molecular Medicine. 2019; 23: 2032–2041. https://doi.org/10.1111/jcmm.14109.

[13] Xiong X, Kuang H, Ansari S, Liu T, Gong J, Wang S, et al. Landscape of Intercellular Crosstalk in Healthy and NASH Liver Revealed by Single-Cell Secretome Gene Analysis. Molecular Cell. 2019; 75: 644–660.e5. https://doi.org/10.1016/j.molcel.2019.07.028.

[14] Li H, Zhou Y, Wang H, Zhang M, Qiu P, Zhang M, et al. Crosstalk Between Liver Macrophages and Surrounding Cells in Nonalcoholic Steatohepatitis. Frontiers in Immunology. 2020; 11: 1169. https://doi.org/10.3389/fimmu.2020.01169.

[15] Lee JS, Hsu YH, Chiu YS, Jou IM, Chang MS. Anti-IL-20 antibody improved motor function and reduced glial scar formation after traumatic spinal cord injury in rats. Journal of Neuroinflammation. 2020; 17: 156. https://doi.org/10.1186/s12974-020-01814-4.

[16] Horn P, Tacke F. Metabolic reprogramming in liver fibrosis. Cell Metabolism. 2024; 36: 1439–1455. https://doi.org/10.1016/j.cmet.2024.05.003.

[17] Cui S, Ghai A, Deng Y, Li S, Zhang R, Egbulefu C, et al. Identification of hyperoxidized PRDX3 as a ferroptosis marker reveals ferroptotic damage in chronic liver diseases. Molecular Cell. 2023; 83: 3931–3939.e5. https://doi.org/10.1016/j.molcel.2023.09.025.

[18] Bi Y, Liu S, Qin X, Abudureyimu M, Wang L, Zou R, et al. FUNDC1 interacts with GPx4 to govern hepatic ferroptosis and fibrotic injury through a mitophagy-dependent manner. Journal of Advanced Research. 2024; 55: 45–60. https://doi.org/10.1016/j.jare.2023.02.012.

[19] Khomich O, Ivanov AV, Bartosch B. Metabolic Hallmarks of Hepatic Stellate Cells in Liver Fibrosis. Cells. 2019; 9: 24. https://doi.org/10.3390/cells9010024.

[20] Wculek SK, Dunphy G, Heras-Murillo I, Mastrangelo A, Sancho D. Metabolism of tissue macrophages in homeostasis and pathology. Cellular & Molecular Immunology. 2022; 19: 384–408. https://doi.org/10.1038/s41423-021-00791-9.

[21] Liang W, Huang X, Shi J. Macrophages Serve as Bidirectional Regulators and Potential Therapeutic Targets for Liver Fibrosis. Cell Biochemistry and Biophysics. 2023; 81: 659–671. https://doi.org/10.1007/s12013-023-01173-w.

[22] Cao Y, Mai W, Li R, Deng S, Li L, Zhou Y, et al. Macrophages evoke autophagy of hepatic stellate cells to promote liver fibrosis in NAFLD mice via the PGE2/EP4 pathway. Cellular and Molecular Life Sciences: CMLS. 2022; 79: 303. https://doi.org/10.1007/s00018-022-04319-w.

[23] Chen S, Saeed AFUH, Liu Q, Jiang Q, Xu H, Xiao GG, et al. Macrophages in immunoregulation and therapeutics. Signal Transduction and Targeted Therapy. 2023; 8: 207. https://doi.org/10.1038/s41392-023-01452-1.

[24] Wen Y, Lambrecht J, Ju C, Tacke F. Hepatic macrophages in liver homeostasis and diseases-diversity, plasticity and therapeutic opportunities. Cellular & Molecular Immunology. 2021; 18: 45–56. https://doi.org/10.1038/s41423-020-00558-8.

[25] Yang H, Cheng H, Dai R, Shang L, Zhang X, Wen H. Macrophage polarization in tissue fibrosis. PeerJ. 2023; 11: e16092. https://doi.org/10.7717/peerj.16092.

[26] Zhang M, Jing Y, Xu W, Shi X, Zhang W, Chen P, et al. The C-type lectin COLEC10 is predominantly produced by hepatic stellate cells and involved in the pathogenesis of liver fibrosis. Cell Death & Disease. 2023; 14: 785. https://doi.org/10.1038/s41419-023-06324-8.

[27] Deshpande K, Martirosian V, Nakamura BN, Iyer M, Julian A, Eisenbarth R, et al. Neuronal exposure induces neurotransmitter signaling and synaptic mediators in tumors early in brain metastasis. Neuro-oncology. 2022; 24: 914–924. https://doi.org/10.1093/neuonc/noab290.

[28] Ding D, Lou X, Hua D, Yu W, Li L, Wang J, et al. Recurrent targeted genes of hepatitis B virus in the liver cancer genomes identified by a next-generation sequencing-based approach. PLoS Genetics. 2012; 8: e1003065. https://doi.org/10.1371/journal.pgen.1003065.

[29] Dobie R, Wilson-Kanamori JR, Henderson BEP, Smith JR, Matchett KP, Portman JR, et al. Single-Cell Transcriptomics Uncovers Zonation of Function in the Mesenchyme during Liver Fibrosis. Cell Reports. 2019; 29: 1832–1847.e8. https://doi.org/10.1016/j.celrep.2019.10.024.

[30] Khan MA, Fischer J, Harrer L, Schwiering F, Groneberg D, Friebe A. Hepatic stellate cells in zone 1 engage in capillarization rather than myofibroblast formation in murine liver fibrosis. Scientific Reports. 2024; 14: 18840. https://doi.org/10.1038/s41598-024-69898-z.

[31] Wang F, Ting C, Riemondy KA, Douglas M, Foster K, Patel N, et al. Regulation of epithelial transitional states in murine and human pulmonary fibrosis. The Journal of Clinical Investigation. 2023; 133: e165612. https://doi.org/10.1172/JCI165612.

[32] Jahn J, Bollensdorf A, Kalischer C, Piecha R, Weiß-Müller J, Potru PS, et al. Microglial CD74 Expression Is Regulated by TGFβ Signaling. International Journal of Molecular Sciences. 2022; 23: 10247. https://doi.org/10.3390/ijms231810247.

[33] Casari M, Siegl D, Deppermann C, Schuppan D. Macrophages and platelets in liver fibrosis and hepatocellular carcinoma. Frontiers in Immunology. 2023; 14: 1277808. https://doi.org/10.3389/fimmu.2023.1277808.

[34] Orecchioni M, Ghosheh Y, Pramod AB, Ley K. Macrophage Polarization: Different Gene Signatures in M1(LPS+) vs. Classically and M2(LPS-) vs. Alternatively Activated Macrophages. Frontiers in Immunology. 2019; 10: 1084. https://doi.org/10.3389/fimmu.2019.01084.

[35] Weng SY, Wang X, Vijayan S, Tang Y, Kim YO, Padberg K, et al. IL-4 Receptor Alpha Signaling through Macrophages Differentially Regulates Liver Fibrosis Progression and Reversal. EBioMedicine. 2018; 29: 92–103. https://doi.org/10.1016/j.ebiom.2018.01.028.

[36] Yamane D, Hayashi Y, Matsumoto M, Nakanishi H, Imagawa H, Kohara M, et al. FADS2-dependent fatty acid desaturation dictates cellular sensitivity to ferroptosis and permissiveness for hepatitis C virus replication. Cell Chemical Biology. 2022; 29: 799–810.e4. https://doi.org/10.1016/j.chembiol.2021.07.022.

[37] Bersuker K, Hendricks JM, Li Z, Magtanong L, Ford B, Tang PH, et al. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature. 2019; 575: 688–692. https://doi.org/10.1038/s41586-019-1705-2.

[38] Chen J, Li X, Ge C, Min J, Wang F. The multifaceted role of ferroptosis in liver disease. Cell Death and Differentiation. 2022; 29: 467–480. https://doi.org/10.1038/s41418-022-00941-0.

[39] Wang Y, Song J, Bian H, Bo J, Lv S, Pan W, et al. Apelin promotes hepatic fibrosis through ERK signaling in LX-2 cells. Molecular and Cellular Biochemistry. 2019; 460: 205–215. https://doi.org/10.1007/s11010-019-03581-0.

[40] Mohamed MR, Haybaeck J, Wu H, Su H, Bartneck M, Lin C, et al. JNKs protect from cholestatic liver disease progression by modulating Apelin signalling. JHEP Reports: Innovation in Hepatology. 2023; 5: 100854. https://doi.org/10.1016/j.jhepr.2023.100854.

[41] Zhao Y, Li M, Yao X, Fei Y, Lin Z, Li Z, et al. HCAR1/MCT1 Regulates Tumor Ferroptosis through the Lactate-Mediated AMPK-SCD1 Activity and Its Therapeutic Implications. Cell Reports. 2020; 33: 108487. https://doi.org/10.1016/j.celrep.2020.108487.

Funding

Traditional Chinese Medicine Administration of Shanxi Province(2023ZYYDA2001)

Health Commission of Shanxi Province(2021081)

Youth Project of Science and Technology Department of Shanxi Province(202303021222369)

Shanxi Academy of Traditional Chinese Medicine(202404)