Single-cell atlas of human liver and blood immune cells across fatty liver disease stages reveals distinct signatures linked to liver dysfunction and fibrogenesis

Abstract

Immune cells play a central yet poorly understood role in metabolic dysfunction-associated steatotic liver disease and metabolic dysfunction-associated steatohepatitis (MASLD/MASH), a global cause of liver disease with limited treatment. Limited access to human livers and lack of studies across MASLD/MASH stages thwart identification of stage-specific immunological targets. Here we provide a unique single-cell RNA sequencing atlas of paired peripheral blood and liver fine-needle aspirates from a full-spectrum MASLD/MASH human cohort. Our findings included heightened immunoregulatory programs with MASH progression, such as enriched hepatic regulatory T cells, monocytic myeloid-derived suppressor cells, TREM2⁺S100A9⁺ macrophages and S100^hiHLA^lo type 2 conventional dendritic cells. Hepatic cytotoxic T cell functions increased with inflammation, but decreased with fibrosis, while acquiring an exhausted signature, whereas natural killer cell-driven toxicity intensified. Our dataset proposes immunological mechanisms for increased fibrogenesis and vulnerability to liver cancer and infections in MASH and provides a basis for a deeper understanding of human immunological dysfunction in chronic liver disease and a roadmap to new targeted therapies.

Extended Data Fig. 1 |. Serum levels of monocyte and macrophage activation marker, soluble (s)CD163.

Levels of sCD163 significantly correlated with (a) non-invasive markers of fibrosis, liver stiffness (kPa) and APRI score; and (b) liver injury enzymes, ALT and AST, two-sided Spearman correlation. (c) Levels of sCD163 were significantly higher in patients with histologically defined MASH, compared to healthy controls (n = 10), and in advanced MASH (n = 9), compared to NS-SS (n = 6); and (d) in patients with high NAS score (>4, n = 11) compared to low NAS score (<5, n = 14) (the box plots represent the median and interquartile ranges, and the whiskers depict the minimum and maximum of the data set), two-sided Mann–Whitney unpaired U test. *P < .05, **P < .01, ***P < .001, ****P < .0001. Abbreviations: kPa, kilopascal; APRI, aspartate aminotransferase–to-platelet ratio index; ALT, Alanine aminotransferase; AST, Aspartate aminotransferase; NAS; MASLD (former NAFLD) activity score. NS, no steatosis/inflammation/fibrosis.

Extended Data Fig. 2 |. Additional information on the quality of scRNA-seq dataset.

(a) Ridge plot of the housekeeping gene B2M expression per sample, per compartment. (b) Stacked bar chart of sample proportion in each cell type, per compartment. (c) Separate clustering UMAPs for liver FNAs (n = 55,234) and PBMCs (n = 123,753), color-coded by cell type (clustering settings: 20 PCs, Louvain Resolution=0.05).

Extended Data Fig. 3 |. Broad characterization of cell type enrichment per compartment and per liver disease stage.

(a) Box plots showing per-sample (n = 25) frequencies of CD16⁻ NK cells, MAIT cells, CD8⁺ T cells and B cells, and the ratio of CD4⁺/CD8⁺ T cells, per compartment: data are presented as median values (horizontal line), 25th/75th quartile values (bounds of boxes), and non-outlier minimum/maximum (whiskers); two-sided Wilcoxon. (b) Box plot of the proportion of each cell type as a percentage of each sample in FNA (left) and PBMC (right): data are presented as median values (horizontal line), 25th/75th quartile values (bounds of boxes), and non-outlier minimum/maximum (whiskers); NS-SS (n = 6), early MASH (n = 10), advanced MASH (n = 9); two-sided permutation test; significance (*) determined by false-discovery rate (FDR) < 0.05 and log2 false discovery (log2FD)>0.58. (c) Scatter plot of average pathway enrichment score per patient (n = 25) of (left) ROS/RNS production in macrophages by controlled attenuation parameter (CAP reflects steatosis), (middle) collagen formation in hepatic stellate cells (HSCs) by MASLD/NAFLD activity score (NAS), and (right) insulin-growth factor binding protein (IGF) and IGF-binding protein (IGFBP) in hepatocytes by liver enzyme ALT: error bands represent 95% confidence interval; two-sided Pearson.

Extended Data Fig. 4 |. Interferon-induced protein validation and liver sources of interferons.

Box plot by cell type per patient of average (a) ISG score per disease stage (NS-SS (n = 6), early MASH (n = 10), advanced MASH (n = 9)) in liver FNA (left) and PBMCs (right), (b) IFNG expression in liver FNAs (n = 25), and (c) upstream gene expression for interferon-alpha (left) and interferon-beta (right) in liver FNAs (n = 25): data are presented as median values (horizontal line), 25th/75th quartile values (bounds of boxes), and non-outlier minimum/maximum (whiskers); two-sided Wilcoxon; upstream genes are shown in Extended Data Table 2. (d) Magnification of Extended Data Fig. 3a: box plot of the percentage of pDCs out of each sample in FNA (left) and PBMC (right): data presented as in 5a; two-sided permutation test; significance (*) determined by false-discovery rate (FDR) < 0.05 and log2 false discovery (log2FD)>0.58. (e) Comparative ISG expression in liver monocytes (left) and macrophages (right) between patients without steatosis or with steatosis (NS-SS), those with MASH and those with untreated HCV: data are presented as median values (horizontal line); two-sided Wilcoxon; ***P < .001. (f) Immuno-fluorescent (IF) validation of MX1 and IFIT3 in macrophages in human liver tissue biopsies. Representative image merged and per channel (left) and expression fractions of MX1 + , IFIT3 + , and MX1 + IFIT3+ out of macrophages (CD68 +) for three disease stages (right): simple steatosis (SS, n = 1), early MASH (F1, n = 1), and advanced MASH (F2, n = 1).

Extended Data Fig. 5 |. Additional data on monocyte association with MASLD/MASH progression.

(a) Monocyte chemotaxis receptors: box plots of average CCRL2 (left) and CMKLR1 (right) expression in FNA by patient: data are presented as median values (horizontal line), 25th/75th quartile values (bounds of boxes), and non-outlier minimum/maximum (whiskers); NS-SS (n = 6), early MASH (n = 10), advanced MASH (n = 9); two-sided Wilcoxon (b) Monocyte chemotaxis ligands: scatter plot of average RARRES2 expression in hepatocytes per patient (n = 25) by CAP score: error bands represent 95% confidence interval; two-sided Pearson.

Extended Data Fig. 6 |. Investigation of the monocyte and monocyte-derived cell paths without specifying a starting point of cell progression.

(a) Velocity analysis of monocyte and monocyte-derived cells (macrophages and cDCs) shown by UMAP in early MASH (left) and advanced MASH (right) and (b) unbiased pseudotime via slingshot analysis, color-coded per cell type.

Extended Data Fig. 7 |. Investigation of T-cell activation marker ICOS, its ligand, and B-cell light chain kappa versus lambda expression.

(a) Immunofluorescent staining of ICOS in human liver biopsies for three disease stages: simple steatosis (SS, steatosis grade 1, n = 1), early MASH (F1, steatosis grade 1, n = 1), and advanced MASH (F2, steatosis grade 2, n = 1). Box plots of (b) average ICOS expression by T-cell subpopulation per patient FNA (n = 25), (c) average T_reg ICOS expression by histological steatosis grade per patient FNA (0–1: n = 10; 2: n = 11, 3 n = 3), two-sided Wilcoxon, and (d) average ICOSLG expression by cell type per patient FNA (n = 25): data are presented as median values (horizontal line), 25th/75th quartile values (bounds of boxes), and non-outlier minimum/maximum (whiskers). (e) Bar chart showing the fraction of B cells that express kappa (IGKC) or lambda (IGLC2–7) light chain genes for liver FNA (left) and PBMC (right).

Extended Data Fig. 8 |. Additional data on immune cell crosstalk.

Changes in predicted individual ligand-receptor interactions between (a) SS versus early MASH and (b) early versus advanced MASH, from HSCs to immune cells (left) and from immune cells to HSCs (right). (c) Radar plot of the total number of cell-cell interactions (both incoming and outgoing) between macrophages and other cell types (left) and between pDCs and other cell types (right).

Fig. 1 |. The landscape of liver and peripheral blood immune cells across the MASLD/MASH spectrum.

a, Graph showing the number of patients (n = 25) spanning MASLD/MASH disease stages, including NS (n = 2, male sex n = 1, median age = 50 (46–54) years), SS (n = 4, male sex n = 1, median age = 47.5 (25–67) years), early MASH (MASH F0, n = 5, male sex n = 2, median age = 60 (56–73) years; MASH F1, n = 5, male sex n = 1, median age = 34 (26–60) years) and advanced MASH (MASH F2, n = 5, male sex n = 2, median age = 53 (31–58) years; MASH F3 and F4, n = 4, male sex n = 1, median age = 61.5 (50–65) years) who provide liver FNA, clinical liver biopsy and peripheral blood. b, Merged Uniform Manifold Approximation and Projection (UMAP) of single cells color coded by cell type (clustering settings: 20 principal components (PCs); Louvain resolution, 0.05) based on analysis of scRNA-seq of 168,634 cells from the liver FNAs and PBMCs as in a. c, Dot-plot of expression of marker genes in each cell type as in b. Shade represents the average scaled expression and size the percentage of cells that express the gene. d, Stacked bar chart showing the proportion of cell types as in b in patients with NS–SS, early MASH and advanced MASH in liver and PBMCs. e, Box plot showing the proportion of CD14⁺ monocytes, CD14⁺CD16⁺ monocytes, CD16⁺ monocytes, cDCs, B cells and macrophages as a percentage of a patient’s liver FNA or PBMC sample (nonparenchymal cells for FNAs). Data are presented as median values (horizontal line), 25th or 75th quartile values (bounds of boxes) and nonoutlier minimum or maximum (whiskers) (NS–SS (n = 6), early MASH (n = 10) and advanced MASH (n = 9); permutation test). *Significance determined by false discovery rate (FDR) <0.05 and log₂(false discovery (FD)) > 0.58.

Fig. 2 |. Monocytes, monocyte-derived cells and their connectivity across the MASLD/MASH spectrum.

a,b, UMAPs of monocytes, macrophages and cDCs color coded by cell type (a) and FNAs versus PBMCs (b) in cells as in Fig. 1b (clustering settings: 15 PCs; Louvain resolution, 0.6). c, Heatmap showing the top differentially expressed genes of each monocyte subcluster as in a. d, Pathway enrichment analysis for each monocyte subcluster as in a. The one-sided hypergeometric test with Bonferroni’s correction was used. DDX58, retinoic acid-inducible gene I (RIG-I); IFIH1, melanoma differentiation-associated protein 5 (MDA5). e, Box plot showing M-MDSC scores in each cell per monocyte subcluster for liver FNA cells as in a. Data are presented as median values (horizontal line), 25th or 75th quartile values (bounds of boxes) and nonoutlier minimum or maximum (whiskers): CD16⁺ mono (n = 3,617), CD14⁺CD16⁺ mono (n = 3,608), immature mono (n = 5,862), S100^hi mono (n = 3,299), ISG^hi mono (n = 2,575), inflammatory mono (n = 299) and differentiating mono (n = 4,116). Two-sided Wilcoxon’s rank-sum test relative to the reference group (S100^hi mono) was used: ^****P < 0.0001. f, Enrichment of monocyte subsets as in a in liver FNA cells versus PBMCs. Data are presented as observed log₂(fold-changes) ± 95% bootstrap confidence intervals (CIs), n values as in f, two-sided permutation test, significance determined as FDR < 0.05 and log₂(FD) > 0.58. g, Box plot showing the proportion of CD14⁺ monocyte subsets (immature mono, S100^hi mono, ISG^hi mono and differentiation mono) in NS–SS (n = 6), early MASH (n = 10) and advanced MASH (n = 9), shown as a percentage of liver FNA (nonparenchymal cells) or PBMC sample per patient. Data are presented as median values (horizontal line), 25th or 75th quartile values (bounds of boxes) and nonoutlier minimum or maximum (whiskers); permutation test: *Significance determined by FDR < 0.05 and log₂(FD) > 0.58. h, Cluster connectivity by PAGA with immature monocytes as the root node in a Reingold–Tilford tree plot layout for monocyte clusters as in a, with CD14⁺ monocyte subsets highlighted in gray. i, Slingshot analysis showing pseudotime progression of the monocyte or monocyte-derived cell lineage, shown by UMAP (left) and strip plot (right). NS, not significant; ZAP-70, ζ-chain-associated protein kinase 70.

Fig. 3 |. Macrophage and pre-macrophage subsets and association with MASLD/MASH progression.

a, UMAPs of macrophages (MPs) and pre-macrophages (PreMPs) color coded by subcluster (top), compartment (FNAs or PBMCs, middle) and disease stage (NS–SS, early MASH and advanced MASH, bottom) in cells as in Fig. 1b (clustering settings: 15 PCs; Louvain resolution, 1.4). b,c, Marker gene expression dot-plot for all PreMP or MP subclusters (PreMP_1, PreMP_2, PreMP_3, MP_0, MP_1, MP_2 and MP_3) (b) and only MP_1, MP_2 and MP_3 (c) as in a. Color represents the average scaled expression of a gene in a cluster and size the percentage of cells that express a gene in a cluster. d, Matrix plot of pairwise cluster modularity index for each PreMP or MP subcluster as in a. e, Subclustering of MP_0 macrophages (clustering settings: 7 PCs; Louvain resolution, 0.4) as UMAP (left) and heatmap of differentially expressed genes (right). f, Differential abundance graph for MPs and PreMPs as in a depicting closely related cells (nodes), the layout of which is determined by the UMAP position of the neighborhood index cell. Node size represents the number of cells in the neighborhood. Lines represent the number of cells shared between adjacent neighborhoods. Color represents the log(fold-change) (log(FC)) differential abundance of neighborhoods in advanced versus early MASH. g, Volcano plot of differentially expressed genes in pooled MP_0, MP_1, MP_2 and MP_3 macrophage subsets in FNAs between early and advanced MASH (left) and the corresponding enriched pathways in early (middle) and advanced (right) MASH. One-sided hypergeometric test with Bonferroni’s correction was used. TriC/CCT, T-complex protein Ring Complex also known as Chaperonin Containing TCP-1 (CCT). h, IF showing staining for S100A9 and CD68 in livers from patients with SS (n = 1), MASH F1 (early MASH, n = 1) and MASH F2 (advanced MASH, n = 1) (top and bottom left) and percentage of S100A9⁺ cells among CD68⁺ cells in patients with NS (n = 1), SS (n = 1), MASH F1 (n = 1) and MASH F2 (n = 1) (bottom right). i, Representative immunohistochemistry images of S100A9 staining in male wild-type C57BL/6J mice fed with HF-CDAA every day for 4 (n = 4), 8 (n = 4) and 12 (n = 4) weeks, starting at 8 weeks of age. ER, endoplasmic reticulum.

Fig. 4 |. DC subsets and association with MASLD/MASH progression.

a,b, UMAPs of DCs color coded by subcluster (a) and compartment (liver FNAs or PBMCs) (b) in cells as in Fig. 1b (clustering settings: 15 PCs; Louvain resolution, 0.13). c,d, Marker gene expression in DCs as in a (c) and heatmaps of differentially expressed genes in cDC1, HLA^hiS100^hi cDC2, S100^hiHLA^lo cDC2 and pDC subclusters (d) as in a. e, Expression dot-plot showing marker genes, including SREBF1/2, LDLR, HMGCSR, HMGCR, MVD, FDFT1, IDH1 and FAS in cDC1, HLA^hiS100^hi cDC2, S100^hiHLA^lo cDC2 and pDC subclusters. The color represents the average scaled expression of a gene in a cluster and size the percentage of cells that express a gene in a cluster. f, Enrichment of cDC1s (n = 86 cells), HLA^hiS100^lo cDC2s (n = 389 cells), S100^hiHLA^lo cDC2s (n = 406 cells) and pDCs (n = 363 cells) as in a in liver FNAs versus PBMCs. Data are presented as observed log₂(FC) ± 95% bootstrap CI; two-sided permutation test: significance determined as FDR < 0.05 and log₂(FD) > 0.58. g, Frequency changes of cDC1, HLA^hiS100^hi cDC2, S100^hiHLA^lo cDC2 and pDC subclusters as in a among total PBMCs or liver FNA cells in relation to MASLD/MASH progression measurements. The two-sided Spearman’s test was used.

Fig. 5 |. T cell subsets and association with MASLD/MASH progression.

a,b, UMAPs of T cells color coded by subcluster (a) and compartment (FNA or PBMC) (b) in cells as in Fig. 1b (clustering settings: 13 PCs; Louvain resolution, 0.6). c,d, Marker gene expression in T cells as in a (c) and heatmaps of differentially expressed genes for CD4⁺ T cell subclusters (left) and CD8⁺ T cell subclusters (right) (d) as in a. e, Expression dot-plot for marker genes in CD8⁺ T_N cell, CD8⁺ T_CM cell, CD8⁺ T_EM cell, CD8⁺ T_EMRA 1 cell and CD8⁺ T_EMRA 2 cell subsets as in a. The color represents the average scaled expression of a gene in a cluster and size the percentage of cells that express a gene in a cluster. f, Enrichment of T cell subsets (MAIT cells (n = 5,891 cells), CD8⁺ T_EM cells (n = 8,186 cells), γδ T cells (n = 3,744 cells), CD4⁺ CTLs (n = 4,227 cells), CD8⁺ T_EMRA 2 cells (n = 6,270 cells), CD8⁺ T_EMRA 1 cells (n = 6,403 cells), CD8⁺ T_CM cells (n = 2,070 cells), T_H17 cells (n = 7,037 cells), naive T_reg cells (n = 1,427 cells), CD8⁺ T_N cells (n = 7,468 cells), T_H1 cells (n = 2,679 cells), CD4⁺ T_CM cells (n = 21,366 cells), CD4⁺ T_N cells (n = 27,813 cells) and memory T_reg cells (n = 4,653 cells)) as in a in liver FNAs versus PBMCs. Data are presented as observed log₂(FC) ± 95% bootstrap CI; two-sided permutation test: significance determined as FDR < 0.05 and log₂(FD) > 0.58. g, Frequency changes of T cell subsets as in a among total PBMC T cells or FNA T cells in relation to MASLD/MASH progression measurements. The two-sided Spearman’s test was used. h, NMF-derived functional programs expression by T cell UMAP plots (top), top genes for each NMF program (middle) with literature-based delineated functions (bottom). GZMH, granzyme H; GZMB, granzyme B. i, Matrix plot of pairwise correlation coefficients among NMF programs. The two-sided Spearman’s test was used. j, Matrix plot showing the median NMF expression of each T cell subset as in a. k, NMF program expression changes of T cells in PBMCs or FNAs in relation to MASLD/MASH parameters. The two-sided Spearman’s test was used. l, Enrichment of each NMF program in FNAs compared with PBMCs. The two-sided Wilcoxon’s test was used.

Fig. 6 |. B cell subsets and association with MASLD/MASH progression.

a,b, UMAPs of B cells and plasma cells color coded by subcluster (a) and compartment (FNAs or PBMCs) (b) in cells as in Fig. 1b (clustering settings: 11 PCs; Louvain resolution, 0.5). c,d, Marker gene expression in B cells and plasma cells as in a (c) and heatmaps of differentially expressed genes per B cell subcluster as in a (d). e, Expression of ICOS ligand per B cell subcluster as in a. f, Frequency changes of B cell subclusters as in a among total PBMCs or liver FNA cells in relation to MASLD/MASH progression measurements. The two-sided Spearman’s test was used. g, Enrichment of plasma cells (n = 296 cells), memory B cells (n = 3,226 cells), naive B cells (n = 9,947 cells), FCRL5⁺ memory B cells (n = 540 cells), ITGB1⁺ memory B cells (n = 1,127 cells) and PLCG2⁺ naive B cells (n = 680 cells) as in a in FNAs versus PBMCs. Data are presented as observed log₂(FC) ± 95% bootstrap CIs; two-sided permutation test: significance determined as FDR < 0.05 and log₂(FD) > 0.58.

Fig. 7 |. NK cell subsets and association with MASLD/MASH progression.

a,b, UMAPs of NK cells color coded by subcluster (a) and compartment (FNAs or PBMCs) (b) in cells as in Fig. 1b (clustering settings: 16 PCs; Louvain resolution, 0.25). c,d, Marker gene expression in NK cells as in a (c) and heatmap of differentially expressed genes in NK cell subclusters as in a (d). e, Expression dot-plot for marker genes in NK cell subclusters as in a. The color represents the average scaled expression of a gene in a cluster and size the percentage of cells that express a gene in a cluster. f, Enrichment of tissue-resident NK cells (n = 2,237 cells), CD16⁻IL7R⁻ NK cells (n = 483), CD16⁺CCL4^hi NK cells (n = 596), CD16⁻IL7R⁺ NK cells (n = 1,549), CD16⁺ NK cells (n = 9,671), CD16⁺PTGDS⁺ NK cells (n = 1,251) and CD16⁺S100B⁺ NK cells (n = 290) as in a in liver FNA cells or PBMCs. Data are presented as observed log₂(FC) ± 95% bootstrap CIs; two-sided permutation test: significance determined as FDR < 0.05 and log₂(FD) > 0.58. g, Frequency changes of NK cell subsets as in a among total PBMCs or liver FNA cells in relation to MASLD/MASH parameters. The two-sided Spearman’s test was used. h, NMF-derived functional program expression by NK cell UMAP plots (top), top genes for each NMF program (middle) and with literature-based delineated functions (bottom). i, Matrix plot of pairwise correlation coefficients among NMF programs. The two-sided Spearman’s test was used. j, Matrix plot showing the median NMF expression of each NK cell subset as in a. k, NMF program expression changes of NK cells in PBMCs and liver FNA cells in relation to MASLD/MASH parameters. The two-sided Spearman’s test was used. l, Enrichment of each NMF program in liver FNA cells compared with PBMCs. The two-sided Wilcoxon’s test was used.

Fig. 8 |. Hepatic immune cell crosstalk and changes with MASLD/MASH progression.

a, Predicted ligand–receptor interactions for each pair of cell types that were captured by liver FNA cells, for NS, SS, early MASH (F0–F1) and advanced MASH (F2–F4). Values of zero (black) reflect insufficient cell numbers for interaction analysis. b, Radar plot of the total number of cell–cell interactions (both incoming and outgoing) between HSCs and other cell types in SS, early MASH (F0–F1) and advanced MASH (F2–F4). c, Ligand–receptor interaction probability per stage between HSCs and CD8⁺ T_EM cells (left) and macrophages (right).