Single Cell Analysis Reveals the Presence of Novel Intermediate Cells in Both Mice and Patients With Severe MASLD.

The authors' responsibilities were as follows: M.M., study design, data analysis and interpretation, manuscript drafting; P.P., A.S., A.C., A.M., E.M. data analysis and interpretation; E.P., M.L., M.B., D.D., S.G., and G.M.T. data generation; D.M., M.I., G.C. acquisition of spatial proteomics images; E.T. cell sorting and data analysis; A.L.F. and F.P. funding acquisition and manuscript revision; P.D. study design, funding acquisition, manuscript revision and has primary responsibility for final content. All authors read and approved the final manuscript.

This work was supported by Ministero della Salute (Italy Ministry of Health), GR‐2019‐12370172, RF‐2021‐12374481, PNRR‐MCNT2‐2023‐12378295, PNC‐E3‐2022‐23683266. Ministry of Education 80185250588. Dipartimenti di Eccellenza.

C57Bl/6 male mice (Taconic Biosciences) were housed at constant room temperature (23°C) under 12‐h light/dark cycles with ad libitum access to water in compliance with the Principles of Laboratory Animal Care (NIH publication 86–23). C57BL/6N mice were fed AMLN diet (40 kcal% fat, 20 kcal% fructose and 2% cholesterol) or standard diet (SD) for 14–22–28 weeks, starting from 6 weeks of age ( n  = 10 mice/group, referred to as steatosis, MASH, MASH‐fibrosis or control group, respectively). Food intake and body weight were recorded weekly. Before sacrifice, mice were fasted for 16 h, and the interventions were done during the light cycle. Blood and liver samples were collected at sacrifice. The complete biochemical and histological characterization of the mouse model has been described in Meroni et al. [ 20 ]. The experimental protocol was approved by Fondazione IRCCS Cà Granda, Ospedale Maggiore Policlinico Milano, and the Italian Ministry of Health Review Boards (protocol 10/2017‐UT). In n  = 2 mice/group, primary mouse hepatocytes and non‐parenchymal cells (NPCs) were isolated by a multi‐step ethylene glycol tetra‐acetic acid (EGTA)/collagenase perfusion technique, through microcannulation of portal vein to ensure efficient perfusion of tissues [ 34 ] (both HEPs and NPCs were isolated from the same mouse). Cells were isolated from the entire liver and centrifuged at 50 g x 3 min three times to separate hepatocytes (on bottom) and NPCs (on top). Next, cell suspensions were FACS sorted by using FACSAria Fusion (FACSAria III) to remove dead cells and cellular debris, exploiting 7‐Aminoactinomycin D (7‐AAD) staining, thus allowing viable cell enrichment. CD45 + positive cells were depleted by hepatocyte fractions. For cell sorting, we applied the following gating strategies: (1) exclusion of doublets and cellular debris (singlets gate); (2) selection of cells based on morphological parameters (morpho) of forward scatter and side scatter (FSC/SSC dot plot) to accurately separate hepatocytes and NPCs; (3) exclusion of cells of leukocyte nature (CD45‐ only for hepatocytes fraction); (4) selection of viable cells (viability). Data was analyzed by FACSDiva (version 8.0.3). Hepatocytes and NPCs live cells were separately loaded on a Chromium Single Cell Instrument (10× Genomics) to generate single‐cell gel beads in emulsion. RNA‐seq libraries were prepared with the Chromium Single Cell Library protocol (10× Genomics) according to manufacturer's instructions, and sequencing was performed on the Novaseq machine using paired‐end sequencing runs for gene expression analysis. Base calling, demultiplexing and adapter masking were performed with Illumina BCL Convert v3.9.3 by IGATech [ https://igatechnology.com ]. Alignment, filtering, barcode counting, and UMI counting were performed using 10× Genomics Cell Ranger v4.0.6 and reference transcriptome GRCm38 (“gex‐mm10‐2020‐A”). The resulting genes‐by‐cell count matrices were analyzed for filtering, normalization, clustering, data integration, UMAP visualization using Seurat V4 [ 10 , 35 ]. These data have been deposited in NCBI Gene Expression Omnibus and are accessible through GEO Series accession number GSE284186 ( https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc : GSE284186 ). Only cells with a number of expressed genes between 3000 (L1) and 7500 (L8), a number of UMIs between 10 000 (L5) and 40 000 (L8), and a fraction of mitochondrial gene less than 20% of the total numebr of expressed genes, were considered (Figure  S2A ). Low‐quality cells were removed. Filtered count matrices were normalized across cells with the centered log ratio (CLR) transformation method and then they were merged using a canonical correlation analysis method to perform pooled analyses. Firstly, cell clustering was performed using Seurat functions FindNeighbours() over the first 20 principal components (PCs) and FindCluster() with resolution equal to 0.8, to optimize cluster stability and number of clusters. Uniform Manifold Approximation and Projection (UMAP) visualization was obtained using the function RunUMAP() on the first 20 PCs. The existence of doublets was assessed using the R package DoubletFinder [ 36 ]. The relative number of cells of each cluster C i and sample S j was quantified by means of the Jaccard index J ij = C i ∩ S j C i ∪ S j . Gene set scores were calculated using the R package scMuffin [ 37 ], which compared to Seurat function AddModuleScore() , provides a more customizable implementation of the underlying algorithm. Cluster enrichment was assessed by the Cell Set Enrichment Analysis procedure. A similar approach has been applied to identify periportal, pericentral, or mid‐lobular hepatocytes by taking advantage of their transcriptional profile and previously described zonation markers. Gene expression changes among clusters were tested using “MAST” Seurat function FindAllMarkers() , setting “min.pct” to 0.1 and “logfc.threshold” to 0.15. Pathways were collected from the KEGG database. Only pathways with a number of genes between 10 and 250 were considered. Over representation analysis was performed using the hypergeometric test and adjusted p ‐values (FDR) were obtained by the Benjamini–Hochberg procedure. To perform gene functional clustering, a “significant pathway similarity” graph was defined considering pathways enriched in DEGs (fdr < 0.1 and at least 3 DEGs) and the overlap coefficient as a similarity measure between two pathways (based on their composition in DEGs). Secondly, clusters of pathways were identified by means of the fastgreedy algorithm. Relationships among cell differentiation states were inferred using RNA velocity (scVelo package, dynamical model) and diffusion maps (scMuffin package). Finally, to clarify the direction of cell differentiation, evolutionary trajectory inference was defined by the dynamical model scVelo. In detail, RNA velocity describes the rate of expression change for an individual gene at a given time point based on the ratio between its spliced and unspliced mRNA [ 16 ]. The MACSima Imaging System is a fully automated platform that combines liquid handling with widefield microscopy for cyclic immunofluorescence imaging. Liver frozen embedded tissues, 4 μm thick, were mounted on SuperFrost Plus slides and fixed with a 4% paraformaldehyde solution for 10 min at RT, rinsed three times with MACSima Running Buffer (130‐121‐565, Miltenyi Biotec, Bergisch Gladbach, Germany), and placed onto MACSwell Four Imaging Frames (130‐124‐676, Miltenyi Biotec). Sections were then pre‐treated with DAPI at a 1:5 ratio (130‐127‐574, Miltenyi Biotec) in the Running Buffer for 20 min, washed three times with Running Buffer, and finalized with a sample volume of 250 μL. The primary fluorochrome‐labeled antibodies were prepared in a MACSwell Deepwell Plate (130‐126‐865, Miltenyi Biotec) using MACSima Running Buffer (Table  S4 ), and then sealed with MACSwell Sealing Foil (130‐126‐866, Miltenyi Biotec) to prevent evaporation. DAPI, at a 1:50 dilution, was then incorporated into the antibody‐filled Deepwell Plate every eighth cycle. Detailed information on image processing was previously described [ 38 ]. In brief, individual raw Field of View (FoV) images were processed using the automated pipeline in the MACS iQ View software (Version 1.3.2). Data are represented as mean ± SD. Statistical analyses were performed using JMP 16.0 Pro (SAS, Cary, NC), R statistical analysis version 3.3.2 ( http://www.R‐project.org/ ) and Prism (version 9, GraphPad Software Inc., San Diego, CA, USA), by using one or two‐way analysis of variance (ANOVA) or chi‐square test, where appropriate, followed by post hoc t ‐tests adjusted for the number of comparisons when multiple groups were involved (Bonferroni correction or Benjamini‐Hochberg false discovery rate (FDR) correction, where indicated). p ‐values < 0.05 (two‐tailed) were considered statistically significant. For sc‐RNAseq analyses, p ‐values were corrected for multiplicity by the Benjamini–Hochberg method, and adjusted p  < 0.1 were considered statistically significant.

To then delve into the fine‐tuned mechanisms underlying steatosis onset/progression and to disentangle cell‐specific complexity and heterogeneity, we sought to perform a sc‐RNAseq on hepatocytes and non‐parenchymal cells (NPCs) separately and isolated from the entire liver of mice fed SD or AMLN for 14–22–28 weeks ( n  = 2/each group). The complete histological and biochemical characterization of this mice model has been recently published [ 8 ], it is represented in Figure  S1A,B and described in Supplementary results. The experimental plan is described in Figure  1A and Figure  S2A . Overall, we obtained about 30 000 single‐cell transcriptomes, including 7424 from control mice, 3300 from AMLN for 14 weeks, 8081 from AMLN for 22 weeks and 10 948 from AMLN for 28 weeks. An integrated Uniform Manifold Approximation and Projection (UMAP) plot has been created, and Seurat was used for clustering analysis according to similarity of expression across all samples [ 9 , 10 ]. We integrated the samples by means of canonical correlation analysis and identified 32 clusters (named from 0 to 31, Figure  1B and Figure  S2B ). Study design, clustering and annotation of cell types based on sc‐RNAseq. Schematic representation of the study design of sc‐RNAseq experiment, from cell isolation, FACS‐sorting, capturing and sc‐RNA sequencing. Hepatic cells were isolated from C57BL/6 mice fed SD or AMLN diet for 14–22–28 weeks ( n  = 2 mice/group) by in vivo perfusion. Hepatocytes (HEPs) and nonparenchymal cells (NPCs) were separated by low‐speed 50 g centrifuge; live cells were enriched by FACS‐sorting and then HEPs (CD45 depleted) and NPCs were separately loaded on Chromium 10× Genomics instrument (A). The Integrated Uniform Manifold Approximation and Projection (UMAP) plot has been created, based on 30.000 single‐cell transcriptomes. 32 clusters have been identified in an unsupervised manner (named from 0 to 31). Each cluster is represented by a unique color (B). Cell type annotation has been performed by considering the gene expression profile similarity, using a gene set scoring system and previously published canonical markers. Heatmap shows differential gene analysis and correlation analysis performed for each cell type (C). The 32 clusters have been subdivided into HEPs (15), CHOLs (1), HSCs (2), ENDO (5), KCs/MoMF (5), N (1), DCs (1), T cells and NKs (1), and B cells (1). Each cell type is represented by a unique color (D). Specifically, by assessing the enrichment of expression markers we found the following liver cell types: HEPs, CHOLs, HSCs, ENDO, KCs, MoMf, N, DCs, T cells, NKs and B cells (Figure  1C,D ). A complete list of canonical markers exploited to identify cell types is presented in Figure  1C and Tables  S1 and S2 . Although the unsupervised analysis has included clusters 10, 17, and 26 in the KCs/MoMf population, cluster 15 in ENDO cells, and cluster 21 in HSCs, they also expressed markers shared by other cell types (lighter blue, turquoise, orange colors in Figure  1C,D ). HEPs represented the largest cell type (16 600 cells, about 50% of the total) and included 15 clusters (Figure  1C,D ; Tables  S2 and S3 ). Clusters 2 and 3, followed by 9 and 6 showed the highest percentage of cells in the control group, and their expression profile could feature hepatocytes in the physiological state. Conversely, they decreased in the steatosis condition in which HEPs were mainly distributed into clusters 4, 11, 12, 14, 18 and 19, poorly represented in the control group (Figure  2A,B ). As expected, these clusters, showing a specific expression of genes involved in fatty acid metabolism and biosynthesis, cholesterol and triglyceride metabolism, bile secretion and TCA cycle, may specifically define the presence of steatosis (Figure  2C,D and Figure  S2D ). Changes in HEP clusters upon AMLN diet exposure. Uniform Manifold Approximation and Projection (UMAP) plots are representative of alterations of cell clusters according to the time of AMLN exposure. Numbers on the UMAP indicate clusters (A, upper panel). HEPs are highlighted in violet in the integrated UMAP. Bar plots show the number of cells (Jaccard index) in HEP clusters across disease severity (A, lower panel). Heatmap illustrates the average expression () of the top markers across the clusters of HEPs (2, 3, 4, 6, 7, 8, 9, 11, 12, 14, 18, 20, 28, 30, 19). The intensity of blue shading is proportional to the induction (B). The average expression ( z ‐score) of the top differentially expressed genes (DEGs) for each experimental condition vs control group was functionally clustered using KEGG pathways for HEPs. Red shading indicates induction, and blue shading indicates repression. Related pathways are clustered together in meta‐pathways by similarity (C). In MASH, clusters 4 and 12 further increased compared to steatosis, cluster 19 disappeared whereas clusters 7, 28, and 30 seemed to be specific only to this condition (Figure  2A ). Pathway analysis revealed that lipid metabolism, TCA cycle, and oxidative phosphorylation were strongly downregulated in MASH compared to steatosis (Figure  2C and Figure  S2D ). To note, in MASH condition, few cells belonging to the cluster 11 of HEPs appeared near cluster 27 of CHOLs in UMAP, probably due to dedifferentiation processes of HEPs into hepatic progenitor cells (HPCs) or oval cells, thus contributing to enforce the CHOLs pool as previously described [ 11 ]. Finally, in MASH‐fibrosis we observed clusters 4 and 7 which were present also in steatosis and MASH, whereas cluster 12, 18, 28 and 30 were almost absent (Figure  2A ). Both the latter two clusters (28 and 30) mainly shut down the overall gene expression compared to the other HEPs, maybe as a result of severe cell injury during MASH, leading to cell death and explaining their disappearance in MASH‐fibrosis (Figure  2A–C ; Figure  S2D ; Table  S3 ). In addition, we could hypothesize a progressive dedifferentiation of HEPs from cluster 9 to cluster 4 passing by 6, since we observed a downregulation of Serpina3k , Apoc3 and Hamp (Figure  2B and Figure  S2C ). The progression in the disease severity was corroborated by the enrichment in pathways such as vascular muscle contraction, retinol metabolism, antigen processing and presentation, and oxidative phosphorylation (Figure  2C ). By assessing zonation markers, we revealed an altered HEPs subset distribution in mouse livers across disease severity (Figure  3A and Figure  S3A ). Indeed, in normal livers and steatosis, some clusters have a pericentral distribution whereas others have a periportal one, reflecting a more physiological liver zonation. Conversely, there was a marked increase in clusters with periportal localization in MASH and MASH‐fibrosis, possibly reflecting the advanced liver damage observed in this area. For instance, in the control group, we found the prevalence of clusters 2, 3, 6, and 9 which showed a mixed pericentral/periportal pattern of distribution, suggested by the co‐expression of Cyp2e1 and Cps1 marker 5 . Changes in ENDO clusters upon AMLN diet exposure. The heatmap represents the average expression ( z ‐score) of zonation markers (periportal vs. pericentral) across different clusters of HEPs. Red shading indicates induction, and blue shading indicates repression. (A) ENDOs are highlighted in blue in the integrated UMAP. Bar plots show the number of cells (Jaccard index) in ENDO clusters across the disease severity. (B) Heatmap illustrates the average expression () of the top markers across the clusters of ENDOs (0, 1,16, 22, 15). The intensity of blue shading is proportional to the induction (C) The average expression (z‐score) of the top differentially expressed genes (DEGs) for each experimental condition vs control group was functionally clustered using KEGG pathways for ENDOs. Red shading indicates induction, and blue shading indicates repression. Related pathways are clustered together in meta‐pathways by similarity (D). During steatosis, we observed the rising of cluster 4, mainly located in the midlobular zone, although the expression of zonation markers was lower compared to that in other HEP clusters in line with Serpina3k , Apoc3 , and Hamp downregulation. Conversely, clusters 11, 12, and 18 had a periportal location (high Cps1 expression), and clusters 14 emerged to be closer to the centrilobular vein (high Cyp2e1 expression). Notably, in MASH condition, clusters 4 and 28 were even more present and showed a midlobular localization, together with the periportal 7, 11, 12, 18, and 30 which became more prominent (high Ass1 and Cps1 expression). Finally, in MASH‐fibrosis clusters 28 and 30 almost disappeared whereas the de‐differentiated cluster 4 and the periportal 7 were still present. We first revealed that during MASLD course novel subtypes of HEPs are generated in different areas across the hepatic lobule, particularly in periportal regions during severe disease. Specifically, the novel type of HEPs which feature steatosis are mainly involved in lipid metabolism, whereas those enriched in MASH and MASH‐fibrosis with a periportal localization may possibly affect oxidative phosphorylation (Figure  S2D and Figure  S3A ). Subpopulation analysis of HEPs belonging to clusters 14 and 11 is represented in Figure  S3B and described in Supplementary results. As previously reported by Xiong et al. 2019 11 , ENDOs represent the largest cell population among the NPCs (Figure  3B ). Nearly 6700 ENDO single‐cell transcriptomes from controls and MASLD livers comprised five clusters (0, 1, 15, 16 and 22), and the total number of ENDOs is remarkably enriched in the MASH‐fibrosis stage, supporting their driving role in the fibrotic damage [ 12 ] (Tables  S2 and S3 ). Cluster 1 was the most representative of the control group, while cluster 16 was enriched in steatosis, and it was maintained also in MASH and MASH‐fibrosis (Figure  2A and Figure  3B ). To note, cluster 1 showed the expression of Fcgr2b and Gpr182 , which usually marked classic Liver Sinusoidal Endothelial Cells (LSECs). Clusters 16 and more so cluster 0 were mainly composed of periportal ENDOs, both positive for Ly6a, Efnb1, Ltbp4 and Ednrb , in addition to genes involved in vascular angiogenesis, including Notch1 and the ligand Dll4 (Figure  3C ). Pathway‐enrichment analysis demonstrated an upregulation of complement and coagulation cascades, PPAR signaling, and biosynthesis of amino acids during steatosis (Figure  3D and Figure  S2E ). Cluster 16 was also present during MASH and expressed genes involved in leucocyte trans‐endothelial migration, as various integrins. Conversely, cluster 15 disappeared, possibly suggesting a transcriptome reprogramming of these cells during this condition (Figure  2A and Figure  3B ). The enrichment of pathways observed in steatosis was confirmed in MASH (Figure  3D and Figure  S2E ). Finally, MASH‐fibrosis was characterized by an increase of the LSEC cluster 0 and compared to cluster 1 and 16, it further induced the expression of genes involved in neoangiogenesis ( Notch1, Dll4 ), VEGF receptors ( Kdr, Flt1, Flt4, Nrp1, Nrp2 ), TGFβ receptors ( Tgfβr3, Bmpr2, Eng, Col4a2, Col4a1 ), Ephrin B receptor ( Ephb4 ), and receptor tyrosine kinase ( Tek, Tie1 ) and phosphatase ( Ptprb ) signaling pathways, in addition to retain classical ENDOs markers (Figure  3C ). Moreover, cluster 0 strongly upregulated genes codifying cytokines and their receptors (i.e., Tnfaip1, Bmp2/6, Il1r1, Il6st, Lifr, Tgfβr3, Bmpr2 ), integrins and those related to cytoskeleton remodeling, thus corroborating its central role during fibrogenesis (Figure  3C ). This condition was featured by the strong upregulation of pathways involved in cytoskeleton remodeling, phagocytosis, RAP1, TGFβ, and chemokine signaling (Figure  3D and Figure  S2E ). Cluster 22, which was composed of a few cells, slightly increased during MASH‐fibrosis, and its expression pattern may indicate a similarity with both pericentral ENDOs ( Rspo3+, Wnt9b+, Cdh13+, Wnt2+ ) and fibroblasts. Indeed, they showed an upregulation of typical markers of fibrocytes, including S100a6 , Smad3/6 , and Ecm1 compared to other ENDOs. Furthermore, cluster 22 was characterized by a higher expression of Klf4, Vwf, Edn1, Nos1ap , which are induced in hypoxic conditions (Figure  3C ). It has been described that hypoxia promotes mesenchymal transition in ENDOs [ 13 ]. According to this possible mechanism named endothelial mesenchymal transition (EndMT), we found that these cells gradually lose the expression of the adhesion molecule PECAM1/CD31 from the control group and steatosis to MASH and MASH‐fibrosis. Conversely, they acquired mesenchymal cell‐specific genes, like extracellular matrix (ECM) macromolecules such as fibulin2/5 and bone morphogenetic protein 4 ( BMP4 ), paralleled by the activation of Wnt signaling and of genes involved in the passaging of immune cells throughout endothelia and in the regulation of vascular permeability, such as Plvap (Figure  3C ). These pronounced alterations may be the result of a severe disruption of the sinusoidal capillaries during MASH and fibrosis. Nearly 4200 KC/MoMf single‐cell transcriptomes from controls and MASLD livers encompass five clusters (5, 10, 17, 25 and 26) (Figure  4A ; Tables  S2 and S3 ). Cluster 5 was present in control; it was enriched in Vsig4  + ve, Clec4f  + ve, Cd5l  + ve, and Ccr2 ‐ve cells, thus specifically defining resident KCs, and it sharply boosted in MASH‐fibrosis, suggesting an enhanced infiltration of macrophages during the course of the disease (Figure  4A ). In particular, cells belonging to cluster 5 were positive for CD86 and HMOX1/2 , which define a M1 phenotype, and they modulated cytoskeleton‐codifying genes mainly in the context of MASH‐fibrosis, probably due to the acquisition of the ability to migrate along hepatic sinusoids (Figure  4B ). Changes in KCs/MoMF and HSCs clusters upon AMLN diet exposure. KCs/MoMF and HSCs are highlighted in light blue and orange in the integrated UMAP, respectively. Bar plots show the number of cells (Jaccard index) in KCs/MoMF and HSCs clusters across the disease severity (A). Heatmap illustrates the average expression () of the top markers across the clusters of KCs/MoMF (5, 25, 10, 17, 26). The intensity of blue shading is proportional to the induction (B). The average expression ( z ‐score) of the top differentially expressed genes (DEGs) for each experimental condition vs control group was functionally clustered using KEGG pathways for KCs/MoMF. Red shading indicates induction, and blue shading indicates repression. Related pathways are clustered together in meta‐pathways by similarity (B). Heatmap illustrates the average expression () of the top markers across the clusters of HSCs (13, 21). The intensity of blue shading is proportional to the induction (C). The average expression (z‐score) of the top differentially expressed genes (DEGs) for each experimental condition vs control group was functionally clustered using KEGG pathways for HSCs. Red shading indicates induction, and blue shading indicates repression. Related pathways are clustered together in meta‐pathways by similarity (C). Cluster 10 was highly present in steatosis and more so in MASH, and it was positive for both immune and hepatocyte markers (Figure  4B ). Conversely, in MASH‐fibrosis, we evidenced a spike in the prevalence of clusters 25 and 26 compared to the other conditions (Figure  4A ). Cluster 25 may possibly represent immune cells recruited from the circulation since they include Adgre1  + ve, Ccr2  + ve, and Cx3cr1  + ve cells, but Vsig4‐ ve and Clec4f‐ ve (Figure  4B and Figure  S4A ). In addition, cluster 25 seemed to have a high replicative potential, evidenced by a great number of Cyclins and Cdks regulating the cell cycle (Figure  4B ). Cluster 26, instead, possibly represented an ENDOs/KCs intermediate since it was positive for markers of both cell types and appeared polarized towards the M2 phenotype, which has been correlated with hepatic injury in MASLD, thus favoring liver remodeling and tissue repair (Figure  S4B ). Lastly, the total number of cells belonging to Cluster 17 was enhanced during MASH‐fibrosis and similarly to Cluster 10, it was positive for both specific markers thus representing a KCs/HEPs intermediate (Figure  4A ). The greatest enrichment of pathways related to KCs/MoMf has been observed in the MASH‐fibrosis group, highlighting those involved in leukocyte trans‐endothelial migration, oxidative phosphorylation, glycolysis/gluconeogenesis, cholesterol metabolism, RAP1, TNFα signaling, and antigen processing and presentation (Figure  4B and Figure  S4C ). A total of 1182 HSCs single‐cell transcriptomes from controls and MASLD livers comprised two clusters (21 and 13) which expressed classical markers ( Reln, Rgs5, Lrat, Lama1, Bmp5, Pdgfrb positives) (Tables  S2 and S3 ). However, cluster 21 is an HSCs/ENDOs intermediate mainly present in healthy livers (Figure  4A ). During the disease worsening, these cells not only overexpressed several collagens, ECM and Pdgfrb signaling, as occurred in the cluster 13 subset, but also LSECs markers and molecules implicated in neoangiogenesis (Figure  4C and Figure  S5A ). Next, to further investigate whether cluster 13 was split into two areas with the smallest one appearing as an “island” in the UMAP, we subdivided it into 6 sub‐clusters, named 13_0, 13_1, 13_2, 13_3, 13_4, and 13_5 (Figure  S5B ). Specifically, we determined that the 13_5 islet defined myofibroblasts, positive for the expression of Igfbp5 , Igfbp6 , S100a6 , Krt19 , and it displayed a higher number of cells in MASH‐fibrosis (Figure  S5C ). Compared to the other HSCs, this subgroup enhanced the expression of collagens and pro‐proliferative genes such as Cdks , Ccns , FOS , and Vim , which may explain its increase alongside disease severity (Figure  4C ). As expected, the larger modulation of pathways featured MASH‐fibrosis, in which HSCs strongly upregulated genes involved in smooth muscle contraction, cytoskeleton remodeling, retinol metabolism, and cytokine production (Figure  4C ). Different intermediate populations have been previously identified by using sc‐RNAseq, although their clear existence and their functions are a matter of debate [ 14 ]. This phenomenon might be driven by the switch on/off of specific signaling pathways as a consequence of noxious stimuli during disease progression. Here, we identified 5 clusters of intermediate populations (10, 15, 17, 21 and 26; Figure  5A ) and to ensure the reliability of them, we evaluated the distribution of doublets and proliferating cells in the integrated UMAP plot (Figure  S6A–D ; Supplementary results). Identification of intermediate cell populations. Distribution of gene set scores for HEPs, ENDO, KCs/MoMF and HSCs marker genes in intermediate clusters (A). The average expression (z‐score) of the top differentially expressed genes (DEGs) for each experimental condition vs. control group was functionally clustered using KEGG pathways for intermediate clusters. Red shading indicates induction, and blue shading indicates repression. Related pathways are clustered together in meta‐pathways by similarity (B). Cluster 10 was composed of KCs ( Vsig4+, Adgre1+, Clec4f+, Mrc1+, Hmox1+, Cd5l +  vs. HEPs) which also expressed HEPs markers ( Cps1+, Cyp2e1+, Alb+, Hamp+, ApoA1+, ApoA2+, ApoB+ vs. KCs). It is peculiar to steatosis and MASH and progressively dampen the expression of classic KCs markers as Vsig4 , Cd5l and Clec4f according to disease severity. Another cluster intermediate between KCs ( Vsig4+, Adgre1+, Clec4f+, Mrc1+, Hmox1+, Cd5l +  vs. HEPs) and HEPs ( Cps1+ Alb + Apoa1+ vs. KCs), was cluster 17 which differently from cluster 10, is prominent during MASH‐fibrosis acquiring the expression of Vsig4 , Cd5l and Clec4f (Figure  4A and Figure  S6C ). This evidence raises the possibility that in certain conditions cells may acquire markers typically associated with a different cell type or change their phenotype in response to environmental cues. Next, we identified cluster 26, which was featured by the expression of both KCs ( Adgre1+, Vsig4+, Clec4f+, Hmox1+, Cd5l+, Esam+, Lyve1+, but Ccr2‐ vs. ENDO) and ENDOs ( Ptprb+, Clec4g+, Pecam1 (CD31)+, Dnase1L3+ vs. KCs) markers. These cells are greatly heightened in MASH‐fibrosis, and they could be useful not only in immune response but also in the preservation of the vascular system by interacting with endothelia during fibrogenesis. Concerning cluster 15, it represented an intermediate between ENDOs and HEPs as it was positive for Ptprb, Clec4g, Pecam1 (CD31), Dnase1L3 (vs. HEPs) and for Cps1, Pck1, Slc7a2, ApoA1/A2, Alb and Bhmt (vs. ENDOs). The latter genes are classically expressed by HEPs with periportal localization. These cells were highly present in control and steatosis groups, whereas they disappeared during MASH (Figure  3B ). Finally, cluster 21 represented an ENDOs/HSCs population as it expressed Ptprb , Clec4g , Kdr +, Ftl1 +, Adam23 +, Pecam1 (CD31), Dnase1L3 (vs. HSCs) and Reln, Lrat, Lum, Rgs5, Ecm1 (vs. ENDOs), which are markers of endothelial and stellate cells, respectively (Figure  4A and Figure  S6C ). We performed pathway enrichment analysis, focusing on intermediate clusters. The biphenotypic behavior of cluster 15 (ENDOs/HEPs) was confirmed by the over‐representation of pathways specific to HEPs, such as cholesterol and amino acid metabolism, and PPAR signaling during steatosis. These processes were downregulated in MASH and MASH‐fibrosis. In the latter condition, the ENDO component predominated, as testified by the enrichment of cytoskeleton remodeling, vascular muscle contraction, and platelet activation, cytokine production through RAP1, and antigen processing and presentation (Figure  5B ). Regarding the double positive KCs/HEPs cluster 10, antigen processing and presentation, leucocyte trans‐endothelial migration, chemokine/TGFβ signaling, phagocytosis, autophagy, and hypoxia, which are biological processes typical of immune cells, were primarily fostered in MASH‐fibrosis. Conversely, the PPAR pathway and cholesterol metabolism predominated in the control group (Figure  5B ). Similarly to cluster 10, the intermediate population 17 (KCs/HEPs) showed a superimposable trend in pathways representation. In addition, we found an enforcement of the lysosome pathway and ECM‐receptor interaction in MASH‐fibrosis (Figure  5B ). Another intermediate between immune cells is cluster 26, which mixes KCs and ENDO characteristics. During fibrogenesis, it boosted pathways distinctive of KCs, like antigen processing and presentation, phagocytosis, autophagy, chemokine signaling, and leukocyte migration. Moreover, the ENDO trait was evidenced by the activation of neoangiogenesis and vasculature remodeling (Figure  5B ). In keeping with the HSCs component of the intermediate cluster 21 (HSCs/ENDOs), we revealed an activation of smooth muscle contraction, cytoskeleton remodeling, and ECM‐receptor interaction during MASH‐fibrosis. On the other hand, the ENDO identity emerged from the induction of neoangiogenesis, regulation of vasculature, and platelet activation (Figure  5B ). To uncover the transcriptional changes that affected intermediate cells and to place them along a temporal trajectory, we applied RNA velocity analysis [ 15 , 16 ]. HEPs (in violet) underwent significant dynamic changes as demonstrated by their velocity lengths and directions, differentiating into various subtypes after AMLN exposure (Figure  6A and Figure  S7A ). Among others, clusters 11, 14, and 12 showed robust directional transitions, thus suggesting an activation of novel gene signatures. Moreover, we also observed a transition from cluster 9 to 6 and then to 4, in line with the progressive loss of hepatocytic marker Serpina3k (Figure  S2C ), which also represented one of the fastest genes in these clusters. Tracking of intermediate cell populations. Evolutionary trajectory inference analysis has been used to explore differentiation patterns of intermediate cell subsets. Arrows track the direction of the inferred directionality of cellular transition and potential future state of cells based on RNA velocity. The transitional velocities are reflected by the size of vectors. The integrated UMAP is split into different cell types (in different colors). Insets in the circles represent the different clusters (A). UMAP plot of clusters 21 and 13, with indication of the 6 sub‐clusters named 13_0, 13_1, 13_2, 13_3, 13_4 and 13_5 in different colors. Diffusion map indicates a specific lineage, showing progression from a progenitor state of HSCs (13_4 in yellow) to a more differentiated/activated one (13_5 and 21, in dark red and blue, respectively). Distribution of “diffusion pseudo‐time” score (dpt) by sub‐clusters of (13). The diffusion map indicates that the more differentiated/activated status of HSCs is observed during MASH‐fibrosis (B). Furthermore, this analysis outlined a trajectory from cluster 4 of HEPs to cluster 1 of ENDOs passing towards the intermediate population 15. The latter displayed high velocity length indicating a rapid change in gene expression, a high confidence which assures that cells in this region are going in the same direction and high latent time (in yellow), resulting in a more differentiated state (Figure  S7A ). Notably, this trend of differentiation was intensified during MASH‐fibrosis (Figure  S7B ). Among the fastest genes supporting this trajectory, there were Ptprb , Adam23 , Ehd3 , Psd3 , and Ushbp1 which are specific ENDO markers. Concerning the other intermediate populations, latent time analysis demonstrated that they had a more elevated differentiation potential, compared to cluster 15 which has reached a final differentiated state (Figure  S7A ). Accordingly, clusters 10 and 17 seemed to be two different stages of the same process of differentiation, and they acquired a more dynamic state in MASH and MASH‐fibrosis although in an opposite direction (Figure  S7B ). In this regard, during MASH the inferred transition was from cluster 5 of KCs to HEPs through cluster 10, corroborating the downregulation of KCs markers ( Vsig4 , Cd5l , and Clec4f ) (Figure  S7B ). Cluster 17, on the other hand, appeared mainly in MASH‐fibrosis, and the differentiation direction was from HEPs to KCs gaining the expression of Vsig4 , Cd5l , and Clec4f . As we expected, for both clusters 10 and 17, the genes driving the velocity entailed Cybb , Clec4f , Fyb , Ccl6 , Pltp , Cd5l , Ctss , Vsig4 , all implicated in immune response. In addition, cluster 10 was a starting point of another trajectory towards cluster 26 (Figure  S7B ). Cluster 21 was designed by RNA velocity analysis as an end point of a differentiation process starting from cluster 13 (Figure  6A ). 13_5 islet was not included in the trajectories due to its low number of cells. Therefore, to better define the role of cluster 21 in HSCs activation, we performed a new clustering on 13 and 21 cells together and then we exploited diffusion map analysis (Figure  6B ). The new clustering was tuned (resolution 0.3) to match the one already performed on cluster 13 (Figure  S5B ). The diffusion map suggested a possible differentiation process starting from sub‐cluster 13_4 and ending in 21 and 13_5. This process was in line with the velocity trajectory from cluster 13 to 21, indicating these cells as a terminal point of differentiation. Indeed, the UMAP plot revealed that sub‐cluster 13_5 might be an alternative endpoint compared to 21 (Figure  6B ). Fastest genes were mainly affected in epithelial to mesenchymal transition and extracellular matrix organization in cluster 21 (ENDOs/HSCs, including Dcn , Bmp5 , Ehd3 , Vipr1 , Pde3a ), whereas antigen processing and presentation and lysosomal degradation were in cluster 26, 10, and 17. To further corroborate sc‐RNAseq data, we spatially mapped the position of the main cell type of the liver across the disease stages by leveraging a MACSIMA approach [ 17 ]. All cell types were identified in hepatic tissue in each condition. The hepatic zonation marked by the pericentral Glul and periportal Cps1 antibodies is more preserved in the control group compared to the pathological conditions (Figures  S8 and S9 ). In the latter, we observed a progressive increase of CD45 (violet) and F480 (white) positive cells. Specifically, we detected a higher prevalence of recruited inflammatory cells, including DCs (MHC‐class II, teal color), MoMf (CD192, light blue), and T cells (CD5) during MASH and MASH‐fibrosis. The latter showed an even stronger positivity for CD163 cells (dark violet) (Figure  S8 ). Next, to confirm the presence of intermediate cells, we investigated the co‐localization of specific markers in the liver microenvironment. We proved the co‐expression of CD31/CD73 (cluster 21), CD31/Cps1 (cluster 15), F480/Cps1 (clusters 10 and 17), and F480/CD31 (cluster 26) (Figure  7 ). Spatial validation of sc‐RNAseq and localization of intermediate cells in liver specimens. Validation of intermediate cell populations across disease severity by analyzing the co‐expression patterns of specific markers using MACSima Imaging System. DAPI staining has been employed to highlight nuclei (scale bar 20 μm; scale bar inset 4 μm). We revealed that CD31/CD73 double‐positive cells (cluster 21) were found mainly in controls, distributed along healthy vessels. Their presence was also detectable in steatosis and MASH in a widespread manner, whereas the reduction of cluster 21 during MASH‐Fibrosis may exacerbate liver damage. Cluster 15 (CD31/CPS1) was mainly present in early stages of the disease, likely located along sinusoids and blood vessels. Both types of CPS1/F480 positive cells (clusters 10 and 17) were prominent during MASH and more so in MASH‐Fibrosis along with the CD31/F480 positive intermediate (cluster 26) that was located close to blood vessels (Figure  7 ). To validate intermediate cells in patients with MASLD, we exploited the MACSIMA approach in hepatic biopsies (normal liver vs MASH‐fibrosis, n  = 2/group). As control, we selected patients who underwent hepatic biopsies for suspected microcirculation defects, negative for steatosis, lobular inflammation, ballooning, and fibrosis. Conversely, patients with MASH‐fibrosis had a histological diagnosis of MASH and significant fibrosis (F2). As expected, hepatic zonation was disrupted in MASH‐fibrosis, also characterized by a massive inflammatory infiltration and a marked positivity for CD163 which tags M2 polarized macrophages (Figures  S8 and S9 ). Intriguingly, we confirmed the presence of intermediate populations also in samples obtained from patients with a pattern of expression similar to that observed in mice. Specifically, clusters 21 (Reelin/CD31) and 15 (Cps1/CD31) mainly feature normal liver with a localization close to the blood vessels. Conversely, the intermediate 10/17 (CD45‐CD11b/Cps1) and 26 (CD45‐CD11b/CD31) were mostly present in MASH‐fibrosis where the former was placed near fibrotic areas and the latter close to the vessels, thus recapitulating the pattern described in mice with similar disease stage (Figure  8 and Figure  S10 ). Spatial visualization of intermediate cells in human liver specimens. Validation of intermediate cell populations in human hepatic samples of normal liver and MASH‐fibrosis by analyzing the co‐expression patterns of specific markers by using MACSima Imaging System. DAPI staining has been employed to highlight nuclei (scale bar 20 μm; scale bar inset 4 μm).

This study provides novel insights into the cellular heterogeneity and transcriptional profiles that underlie the progression of MASLD towards advanced disease. By exploiting sc‐RNAseq paired with a spatial proteomic approach, our model was enabled to reveal novel cell populations and dynamic changes occurring in HEPs, ENDOs, HSCs, and immune cells across the spectrum of the disease, thus further corroborating the notion that the evolution of liver injuries is driven by complex interactions between the cellular milieus. Here, we exploited the AMLN murine model of MASLD, which recapitulated the human MASLD, from simple steatosis to MASH and MASH‐fibrosis according to time of exposure [ 18 , 19 , 20 ]. To delve into the cellular and molecular heterogeneity underlying the disease course, we handled a sc‐RNAseq experiment on HEPs and NPCs separately to amplify the resolution of the nonparenchymal populations 12 , thus minimizing the possibility of cross‐contamination. Our strategy allowed us to recognize 32 distinct cell populations in the liver, classified in different types of HEPs, KCs, HSCs, ENDOs, CHOLs, B cells, T/NK cells, and neutrophils, based on canonical markers [ 21 ]. By integrating data obtained from mice at different time points, we detected a higher number of cell clusters compared to the previous studies [ 6 , 12 , 21 , 22 ], half of which were composed of hepatocytes. Indeed, we classified 15 clusters of HEPs, with a significant variation in their quantity and in the transcriptional profile across the different disease stages [ 11 ]. Notwithstanding, differently from Su Q. et al. [ 21 ] which obtained few hepatocytes from the latter stages of the disease by collecting them from the NPC suspension, our experimental approach allowed us to identify a greater number of HEPs in MASH‐fibrosis, thus indicating the efficacy of the perfusion/isolation process. In this regard, our findings underscored even the impact of advanced disease on liver zonation. In detail, we showed a subversion in hepatocyte zonation upon the AMLN diet, with a shift towards periportal clusters, especially in MASH and MASH‐fibrosis, thus suggesting a spatial reorganization due to severe damage in these regions. Conversely, in normal liver some clusters have a pericentral distribution whereas others have a periportal one, reflecting a more physiological liver zonation. Accordingly, a recent study emphasized that the importance of hepatic zonation relies on the ability of HEPs to adapt their metabolic activities according to the position throughout the lobule, thus defining hepatocyte heterogeneity [ 5 , 20 ]. The latter adaptability acquired even more relevance under stress conditions such as severe MASLD. Furthermore, the disappearance of specific clusters and the appearance of others during disease progression could reflect HEPs dedifferentiation or cell death [ 23 , 24 ]. Trajectory analysis demonstrated that HEPs were the main cell population who underwent a substantial change in transcriptomic profiles. Indeed, it showed that some clusters (i.e., 11 and 12) acquired novel gene signatures including those involved in lipid and cholesterol metabolism. In addition, it indicated the transition of HEPs towards different cell types and finally confirmed the de‐differentiation process involving clusters 9 to 4. The latter, from a dedifferentiated status, could possibly move towards cluster 1 of ENDO, passing by cluster 15. During this transition, the expression of specific endothelial genes ( Ptprb , Adam23 , Ehd3 , Psd3 and Ushbp1 ) was acquired, thus representing periportal hybrid hepatocytes or progenitor‐like cells located at the endothelial‐hepatocyte interface. The ENDOs/HEPs intermediate cells may likely arise from physiological hepatocyte–LSEC zonation in the healthy liver, where parenchymal and sinusoidal endothelial cells engage in tightly coordinated angiocrine and metabolic programs along the porto‐central axis. This zonal program is progressively disrupted during MASLD progression, leading to the loss of this intermediate transcriptional state. The existence of the intermediate ENDOs/HEPs population was testified also by Su Q. and Xiong et al. who did not report any specific distribution of these cells at histological level [ 21 , 25 ]. The dedifferentiation process has also been described in primary hepatocytes which lose their specific functions, changed morphology, and downregulated metabolic enzymes [ 26 , 27 ]. Similarly, the reduction of synthetic capacity has been observed also during chronic liver injury and liver tumors [ 28 , 29 ]. Among NPCs, the largest cell population is constituted by ENDOs, which were outlined in 5 clusters, characterized by a different prevalence during the disease. The enrichment of ENDOs in the latter stages highlighted the importance of these cells in the fibrosis onset and pinpointed the pivotal role of EndMT as a contributor to liver fibrogenesis. Indeed, during the disease worsening, ENDO cells seemed to lose their classical markers such as PECAM1 (CD31), gaining in turn mesenchymal ones and upregulating genes involved in collagen formation, ECM proteoglycans, neoangiogenesis, fibrovasculogenesis, and TGFbeta signaling. Alongside, KCs and MoMf were represented in 5 clusters; cluster 5 was enriched in classical resident KCs with M1 phenotype, and it was progressively enlarged with disease severity, even acquiring migration properties. In contrast, cluster 10 was defined by a mixed immune‐hepatocyte phenotype and was dominant in steatosis and MASH, probably to aid tissue repair. During MASH, increased hepatocytes death and steatotic debris raise the burden of apoptotic bodies and lipid droplets that are cleared by lipid associated macrophages (TREM2+ LAM) [ 30 ]. We found that cluster 10 was TREM2+ possibly representing LAM that have engulfed debris of apoptotic or necrotic hepatocyte fragments thus acquiring a combined transcriptional signature. Another type of intermediates between macrophages and hepatocytes was cluster 17, which was enhanced in MASH‐fibrosis. In the latter stage, it gained the expression of Vsig4, Cd5l and Clec4f , to promote tissue repair and scar formation [ 22 , 31 ]. Thus, it may be associated with the acquisition of a more inflammatory phenotype of hepatocytes (by expressing MHC, TLRs, cytokines), related to tissue remodeling. Conversely, clusters 25 and 26 emerged in MASH‐fibrosis, potentially representing circulating immune cells and ENDOs/macrophage intermediates, respectively. In particular, cluster 25 appeared to be positive for Ccr2 and Adgre1 , supporting their recruitment from blood circulation. Instead, cluster 26 may exert a dual role in immune response and vascular maintenance, since it expresses classical genes of both cell types. The occurrence of the intermediate cluster 26 is consistent with recent data [ 6 , 21 , 25 ]. These cells have been discovered in severe injuries and in liver tumors, acquiring a crown‐like disposition around injured HEPs. They seem to be phenotypically closer to tumor‐associated macrophages and, in keeping with this notion, we found them strongly enhanced in severe disease. Existing studies demonstrated that bone‐marrow (BM)‐derived progenitor cells (earliest precursors of erythrocytes, megakaryocytes, and macrophages) may provide a complementary source of ENDOs, thus contributing to the formation of liver sinusoidal vessels after endothelial damage [ 32 ]. Accordingly, we found them strongly enhanced in severe disease, possibly representing a particular type of macrophages adhering to damaged endothelial cells in response to advanced disease. Then, we identified 2 clusters of HSCs, both expressing classical markers but with peculiar features. Cluster 21 had an ENDO‐HSC bi‐phenotype, indicating a novel HSC subtype with vasoactive properties affecting both vascular and extracellular matrix remodeling. The existence of these cells has been already supported by Su Q. et al. who demonstrated that although their number was higher in physiological conditions, they participated in liver damage modifying their localization similarly to activated HSCs [ 21 ]. Notably, it should be considered that HSCs are more difficult to isolate from fibrotic livers, and we cannot rule out that we have partially lost a great number of cells in the latter stages during the sc‐RNAseq procedure. Concerning the other HSCs cluster, we subdivided the 13 into six distinct populations, and in particular the 13_5 was classified as myofibroblasts. These cells paralleled the disease severity by heightening their number during MASH‐fibrosis and reflecting their role in collagen production. Diffusion map indicated that clusters 13_5 and 21 may be alternative endpoints of differentiation during HSCs activation. Interestingly, the identification of novel cell populations bearing intermediate phenotypes, exhibiting features of both ENDOs and HSCs or HEPs and KCs, is noteworthy. Our findings overcome the previously reported, since we highlighted a greater number of intermediate populations compared to previous sc‐RNAseq experiments, including two novel intermediates between HEPs and KCs. Applying RNA velocity, which leverages the ratios of spliced versus unspliced RNA molecules, we attempted to predict the fates of these cells. We determined that while cluster 15 had reached a defined differentiated state from HEPs to ENDO cells, clusters 10 and 17 had a greater differentiation potential, shifting from KCs to HEPs and vice versa . Under pathological conditions, ENDOs and HEPs adopt a more pro‐inflammatory and pro‐fibrotic phenotype contributing to fibrosis and potentially interacting with immune cells and HSCs [ 33 ]. The discovery of these intermediate cells raises the idea that they may represent transitional states occurring as a maladaptive response to chronic injuries and stress. We can hypothesize that ENDOs‐HSCs hybrids may be involved in the initiation and expansion of fibrosis, acquiring both pro‐fibrotic and pro‐inflammatory properties characteristic of activated HSCs and ENDOs, respectively. Likewise, intermediates between HEPs and KCs or HEPs and ENDOs may perpetuate inflammation combining metabolic and immune capabilities. This double identity could supply cells with the ability to respond to microenvironmental cues, making them key players in the disease course. Spatial proteomics confirmed the co‐localization of specific markers in the liver microenvironment, thus validating intermediate cells. Cluster 21 was found mainly in controls, distributed along healthy vessels, and may be exerting anti‐inflammatory effects. Likewise, cluster 15 featured the early stages of the disease, and it was located along sinusoids and blood vessels, in the attempt to compensate for endothelial dysfunction. As suggested by trajectory analysis, in the latter stages, this cluster is poorly detectable since it fully differentiated into ENDOs, thus probably participating in vascular remodeling. Both clusters 10 and 17 portrayed MASH and MASH‐Fibrosis respectively, thus reflecting an active immune response to hepatocyte damage and indicating a hepatocyte‐immune crosstalk. Finally, we found cluster 26 in MASH‐fibrosis close to blood vessels where immune responses and vascular damage converge. Notably, we confirmed the presence of intermediate populations also in samples obtained from MASLD patients with a pattern of expression and localization similar to that observed in mice. The main limitation of this study consists in the lack of lineage tracing and functional assays to fully understand the origin and function of intermediate populations, although we have firstly confirmed their presence also in human livers with a spatial localization comparable to that observed in mice. In conclusion, this study ameliorates our understanding of MASLD pathophysiology by highlighting novel complex cellular scenarios and revealing intermediate cell populations possibly driving the disease worsening. The integration of single cell technology implemented with RNA velocity analysis and spatial proteomics was crucial to unraveling the cellular network which participates in MASLD, to clarify the dynamic changes of hepatocytes during the disease course and to identify novel cell subtypes that may contribute to MASLD evolution.

Metabolic dysfunction‐associated steatotic liver disease (MASLD) is the most common chronic liver disease worldwide affecting around 25% of the adult population [ 1 , 2 ]. Over the past decades, it has become a global health challenge due to its increasing diagnostic incidence, even in pediatrics [ 3 ]. MASLD encompasses a wide spectrum of clinical phenotypes spanning from steatotic liver disease (SLD) to metabolic dysfunction‐associated steatohepatitis (MASH), which is characterized by inflammation and hepatocellular degeneration, fibrosis, cirrhosis, and hepatocellular carcinoma (HCC) [ 4 ]. The pathogenesis of MASLD consists of an intricate series of changes which occur in the different cell types of the liver, thus precipitating simple steatosis into MASH and fibrosis. However, the fine‐tuned mechanisms regulating the transition towards more severe forms of MASLD are not completely clarified. Specifically, a deep understanding of the heterogeneity of cell populations involved in the processes succeeding fat deposition into the hepatic cells, of their peculiar transcriptomic profiles and cross‐talks remains so far [ 5 ]. In the era of ‘omics’, massive RNA‐sequencing data is no longer adequate to meet the need to disentangle this mess. Indeed, bulk RNA‐sequencing analysis does not provide an indication of the different sub‐types of resident or infiltrating cells driving the onset of inflammation and fibrogenesis and their pattern of distribution, thus limiting our insights. For instance, it is well established that hepatocytes (HEPs) can exhibit various metabolic and functional attributes based on their transcriptional profiles and localization in the lobule architecture, thus defining their cellular identities. Similar properties have also been ascribed to the other hepatic cell types. For example, endothelial cells (ENDOs) may be featured by heterogeneity thus expressing markers of Kupffer cells (KCs) and playing a role in tumor angiogenesis [ 6 ]. Overall, the cellular composition defines the cell atlas [ 7 ], which may outline the healthy liver. Hence, it can consequently be assumed that variations in the cellular landscape, in terms of quality, quantity, transcriptional activity, and localization may delineate the different stages of the disease. Therefore, we aimed to investigate at single cell resolution these differences in cell number, type and transcriptomic signature underlying the switching from simple steatosis to MASH and fibrosis by exploiting a murine model of high fat high fructose high cholesterol (HFHFHC) diet, also referred to as Amylin liver NASH (AMLN) diet, to resemble human MASLD [ 8 ]. Here, we leveraged a single cell RNA sequencing (sc‐RNAseq) approach coupled with a spatial proteomic validation to identify and map novel cell clusters with an aberrant phenotype which populate livers in the context of progressive disease. In particular, by considering their transcriptional profiles, we clustered 32 cell populations, including distinct types of HEPs, ENDOs, KCs, monocytes and monocytes‐derived macrophages (MoMf), hepatic stellate cells (HSCs), cholangiocytes (CHOLs), B cells, T cells/natural killer (NK) cells, dendritic cells (DCs) and neutrophils (N). Next, we identified novel disease‐related cells with an intermediate phenotype between endothelial and stellate cells, hepatocytes and KCs, and we tracked their localization in the parenchyma across disease stages.

The authors declare no conflicts of interest.

Figure S1: fsb271673‐sup‐0001‐FigureS1.tif. Figure S2: fsb271673‐sup‐0002‐FigureS2.tif. Figure S3: fsb271673‐sup‐0003‐FigureS3.tif. Figure S4: fsb271673‐sup‐0004‐FigureS4.tif. Figure S5: fsb271673‐sup‐0005‐FigureS5.tif. Figure S6: fsb271673‐sup‐0006‐FigureS6.tif. Figure S7: fsb271673‐sup‐0007‐FigureS7.tif. Figure S8: fsb271673‐sup‐0008‐FigureS8.tif. Figure S9: fsb271673‐sup‐0009‐FigureS9.tif. Figure S10: fsb271673‐sup‐0010‐FigureS10.tif. Figure S11: fsb271673‐sup‐0011‐FigureS11.tif. TableS1‐S4: fsb271673‐sup‐0011‐TableS1‐S4.docx.