Spatial proteogenomics reveals distinct and evolutionarily conserved hepatic macrophage niches

1The immense advances in single-cell transcriptomics have enabled a better understanding of the cellular composition of different organs across species. However, we still lack information regarding how these cells are organized in their distinct microenvironmental niches. Moreover, the specific cell-cell interactions determining the identity of individual cells within tissues remain to be defined. While the spatial organization of hepatocytes within the liver is understood, that of non-parenchymal liver cells remains unclear. This is the case for the mouse liver, but even more so for the human liver, where the identity and the precise localization of most hepatic cells is unknown. Moreover, the link between the transcriptome and the proteome has not been studied, resulting in a lack of reliable surface markers to identify these cells by flow cytometry and confocal microscopy. Here, we used proteogenomic techniques including cellular indexing of transcriptomes and epitomes by sequencing (CITE-seq) and spatial approaches to identify all cells and their specific locations within the healthy and obese livers of mice and humans. By doing so, we have developed strategies for the identification and further study of hepatic cells. Demonstrating the usefulness of this approach, we identify the evolutionary conserved and spatially restricted signals driving the distinct hepatic macrophage phenotypes.

2To generate a proteogenomic atlas of the liver, we first examined the optimal method for retrieving all hepatic cells. Using the murine liver, we compared single-cell RNA sequencing (scRNA-seq) using cells isolated via ex vivo or in vivo enzymatic digestion with single-nuclei RNA sequencing (snRNA-seq) (Figures S1A–S1C). We did not observe any differences in the number of genes/cell between the two digestion methods (Figures S1D and S1E), but snRNA-seq typically yielded a lower number of genes/cell (Figures S1D and S1F). This did not prevent distinct cell types from being identified in the snRNA-seq dataset as both scRNA-seq and snRNA-seq identified highly expressed genes in each population. However, expression was often higher in the scRNA-seq (Figure S1F). Additionally, we observed a signature of digestion-associated genes and snRNA-seq-associated genes across cell types (Figure S1F). While in terms of genes/cell snRNA-seq is inferior to scRNA-seq, this method best recapitulated the cell frequencies observed in vivo (Figures S1D–S1N). As each method has advantages and disadvantages, the optimal method to use depends on the biological question being addressed. Here, as we sought to generate a proteogenomic atlas of all hepatic cells, we thus opted to use a combination of all protocols. To investigate mRNA and protein expression at single-cell resolution, we used CITE-seq, staining a selection of the scRNA-seq samples with 107–161 oligo-conjugated antibodies (Figure 1A). Data were pooled together for a single analysis where, with TotalVI, both the protein and mRNA profiles were considered for clustering (Figure 1B). Analysis of the differentially expressed genes and proteins (DEGs/DEPs; Figures S2A and S2B; Table S1) identified 17 cell types (Figure 1B), which were differentially represented with each isolation method (Figure S2C). Addition of antibodies in the CITE-seq analysis enabled surface markers for all cells to be identified, including VSIG4 and FOLR2 for Kupffer cells (KCs) (Figure S2B), without affecting the quality of the transcriptomic data (Figure 2D). This analysis identified one subset of KCs. However, 2 subpopulations termed KC1 and KC2 have recently been described. As our CITE-seq analysis identified that the markers used to identify KC2s, namely CD206 and ESAM, are largely expressed by liver sinusoidal endothelial cells (LSECs) (Figure S2E), we next sought to determine if we had previously removed any potential KC2s in our initial QC steps as LSEC-KC doublets. To examine this, we generated a UMAP of samples containing CD206 and ESAM in the CITE-seq panel and performed a minimal QC filtering on number of genes and percentage of mitochondrial genes. In this UMAP, we found multiple subpopulations including 3 populations of KCs and 2 populations of B cells (Figure S2F). To determine if any of these populations could be KC2s, we harnessed the power of CITE-seq to recreate the gating strategy recently proposed to identify KC2s. Converting the CITE-seq data into an FCS file and analyzing this in FlowJo identified 2 of the KC populations to be KC2s expressing CD206 and ESAM, and these cells were present in a similar ratio among KCs as reported by Blériot et al. (Figure S2G). However, employing the same gating strategy with the B cells also enabled a B cell2 population to be identified (Figure S2G). The KC2 and B cell2 populations also expressed other protein markers associated with LSECs including CD26, CD31, and CD38, suggesting that these may be doublets (Figure S2H). Consistent with this, we did not uncover any DEGs specifically expressed in the KC2s or B cell2s, rather these cells had an intermediate profile between either KC1s or B cell1s and LSECs (Figure S2I), as would be expected for doublets. Finally, we perfused the livers with antigen fix to inflate the LSECs to be able to distinguish more readily between KCs and LSECs and performed confocal microscopy. 3D reconstruction of these images indicates that CD206 expression is observed primarily in the LSECs with which the KCs are intertwined. These microscopy images also show that KCs express some CD206, however, consistent with our CITE-seq analysis, this was observed across the KC population rather than in a KC subset further suggesting the existence of only one-KC population (Figure S2J). With this in mind, we decided to continue to apply our initial strict QC controls eliminating these potential doublets from further analyses.

3To locate the cells identified we performed spatial transcriptomics analysis using Visium. For this, we cut the liver in two distinct orientations to profile both the liver tissue and the capsule (Figures 1C and S2K). We ordered each Visium spot along a spatial trajectory, and annotated portal, periportal, mid, and central zones based on known hepatocyte zonation markers (; Figures 1D, 1E, and S2L) and confirmed this annotation using confocal microscopy (Figure S2M). By using the reference sc/snRNA-seq data, we then deconvolved each spot into its constituent cell types and investigated how cell abundance changed with zonation (Figures 1F, 1G, and S2N). Validating this approach, cholangiocytes mapped specifically to the portal zones, while KCs were preferentially located in periportal and mid zones. While KC location is zonated, we did not identify a strong zonation pattern in the gene expression profiles of KCs (data not shown). We further identified B cells, T cells, endothelial cells (ECs), and stromal cells (SCs) across all zones, while conventional dendritic cells (cDCs) were found at the portal vein (PV), with a minor presence at the central vein (CV) (Figures 1F, 1G, and S2N).

4To validate these locations at single-cell resolution, we next sought to identify the best cell-specific surface markers that would also work by confocal microscopy. As the fixation step utilized for confocal microscopy often affects the integrity of protein epitopes, it is not possible to predict which antibodies will work spatially on fixed tissue slices. Therefore, to simultaneously screen multiple antibodies to identify those working by microscopy, we performed a second Visium analysis which we complemented with 100 oligo-conjugated antibodies, chosen based on the CITE-seq results (Figure 1H). The antibodies identified to work spatially were then validated at single-cell resolution, using MACSima™ Imaging Cyclic Staining (MICS) technology and a 60-plex antibody panel (Figure 1I). Unfortunately, we could not identify useful surface markers for all populations. For example, we did not identify enough discriminatory surface markers that worked by confocal microscopy to distinguish the cDC subsets. Notably, while CD103 has previously been proposed to distinguish hepatic cDC1s from cDC2s, our data demonstrate that this is only expressed by a fraction of these cells (Figure S2O). To confirm the locations of the cDC subsets, we thus turned to Molecular Cartography™ (Resolve BioSciences) that allows for 100-plex spatial mRNA analysis. Genes were selected based on the DEGs from the sc/snRNA-seq data that were also spatially resolved according to Visium. We also identified the portal-central trajectory in this dataset using cholangiocytes genes (Epcam and Spp1) and known zonated hepatocyte genes (Glul, Cyp2e1, Hal, and Sds; Figure S2P). Using expression of Xcr1, Clec9a (cDC1s) and Cd209a, Mgl2, and Clec10a (cDC2s), we confirmed that both cDC1s and cDC2s were localized primarily at the PV (Figure 1J). As cDC2s shared a number of genes with monocyte-derived cells (Figure S2A), we also examined the expression of general monocyte/mac (Cd14, Adgre1, Axl, Mafb, Cx3cr1, and C5ar1) and KC-specific genes (Cd5l and Vsig4) to further validate their identification as bona fide cDC2s (Figure 1J). The punctate nature of mRNA expression in these analyses combined with the dendritic shape of myeloid cells renders it difficult to convincingly determine cell boundaries and to conclude these cDC2s and macs were distinct cells. To validate this, we therefore developed a protocol that combines mRNA detection (RNAScope) with surface protein detection. Examining expression of cDC- or mac-specific mRNAs combined with protein surface markers confirmed the presence of PV cDC1, cDC2s, and non-KC macs (Figure 1K). Taken together, by combining multiple spatial transcriptomic and proteomic approaches, we located all the cells within the murine liver and identified additional heterogeneity within the myeloid cells, not revealed when examining the sc/sn-RNA-seq dataset in isolation. This highlights the power of combining single-cell and spatial proteogenomic techniques to investigate cellular heterogeneity.

5To better understand the non-KC macs, we zoomed in on myeloid cells (cDCs, KCs, monocytes, and monocyte-derived cells) in our sc/snRNA-seq analysis defining 11 populations (Figures 2A–2C; Table S2). This included KCs, 3 populations of non-KC macs and cells that had a profile intermediate between monocytes and patrolling monocytes or macs, termed transitioning monocytes. Closer inspection of the non-KC macs identified cluster6 as peritoneal macs (Figure 2B). The DEGs between the remaining populations suggested that cluster7 likely resembles capsule macs, expressing Cd207 and Cx3cr1 while cluster8 resembles Gpnmb+Spp1+ lipid-associated macrophages (LAMs) we recently described in the fatty liver (; Figures 2A–2D). Conversion of the CITE-seq data into an FCS file allowed an in silico gating strategy to be defined (Figure S3A). Validating this, we utilized the strategy to FACS-purify the populations and assess gene expression (Figures S3B–S3D). Fitting with a recent report, washing the liver prior to digestion enriched the peritoneal macs in the wash fraction, demonstrating these were contaminants on the liver surface rather than being present in the liver tissue itself (Figure S3E). While the CITE-seq markers did not discriminate between cluster7 and cluster8, adding CD207 to the panel enabled the non-KCs to be divided into CD207+ and CD207− macs (Figure S3F). Fitting with their designation as capsule macs, the relative abundance of CD207+ macs was increased if we dissected and digested the capsule (Figure S3F). However, although Molecular Cartography confirmed the presence of Cd207+ macs in the capsule, it also revealed Cd207+ macs at the CV, which were rarely found at the PV (Figures 2E–2H and S3G–S3J). Thus, cluster7 consists of both capsule and CV CD207+ macs. This finding further demonstrates the need for spatial approaches to confirm cell identities. Molecular Cartography also identified macs at the PVs and CVs expressing Ccr2 and Chil3 (Figures 2G, 2H, S3H, and S3J), resembling transitioning monocytes (cluster11). Finally, a population of Gpnmb-expressing macs were found to be enriched around the bile ducts (Figures 2G, 2H, S3H, S3K, and S3L). As Gpnmb expression is cluster8 specific (Figure 2B), and these cells resemble LAMs, we termed these cells bile-duct LAMs.

6We next sought to investigate the differences between KCs and LAMs. Analysis of GO terms associated with biological processes for these cells suggested that KCs may play a role in regulating humoral responses, while LAMs were more broadly associated with immune responses (Figures 2I and 2J). Consistent with this, in the 100-plex protein microscopy data we noted that a significant proportion of the B cells present were interacting with KCs, which was not observed with T cells (Figure 2K). This suggests cross-talk between these two populations, potentially linked to the high expression of the B cell chemokine Cxcl13 by murine KCs (Table S2). To assess the inflammatory nature of LAMs compared with KCs, we FACS-purified the cells and performed qPCR analysis to examine expression of various cytokines. To enable LAM purification, we eliminated the capsule prior to digesting the tissue. Fitting with the GO analysis, LAMs expressed more Il1b at steady state compared with KCs (Figure 2L). However, despite this, upon in vivo TLR4 stimulation, they were less responsive than KCs, both in terms of pro- and anti-inflammatory cytokines (Figure 2L), possibly indicative of LPS tolerance. This may result from their location at the PV and hence exposure to blood from the intestine, although this remains to be tested. Taken together this highlights the distinct nature and functions of these cells; however, further research is required to determine the precise roles of these cells.

7As all the mac populations are in close contact with CD45− cells in their local environment (; Figure S3M), we further analyzed the CD45− cells, identifying multiple subsets of ECs and SCs and a gating strategy to distinguish them (Figures 3A–3C and S4A–S4C; Table S3). ECs could be further subdivided into 4 distinct clusters and analysis of their locations allowed them to be identified as CV ECs (cluster10), LSECs (cluster9), PV ECs (cluster11), and lymphatic ECs (LECs; cluster12) (Figures 3D, 3E, S4D, and S4E). As Visium found fibroblasts at both the PVs and CVs (Figure 3D), and as a previous report has suggested the presence of distinct subsets within these cells, we further zoomed in on the SCs to better assess their heterogeneity (Figures 3F, 3G, S4A–S4C, and S4F; Table S4). This revealed subsets of mesothelial cells and fibroblasts restricted to the capsule (Figures 3H and 3I). Myh11+ vascular smooth muscle cells (VSMCs) were localized around hepatic arteries, PVs and CVs (Figures 3F–3H, 3J, and S4G) and Mfap4+Svep1+Clic5−Reln− fibroblasts were found to be CV fibroblasts (Figures 3G, 3H, 3J, and S4G). Finally, we identified a subset of Clic5+Reln+ fibroblasts (cluster3), which were localized around the cholangiocytes, that we termed bile-duct fibroblasts (Figures 3J and S4G). Taken together, the presence of these spatially distinct subsets of ECs and SCs highlights the uniqueness of the specific microenvironments in which the distinct mac populations reside.

8Finally, in addition to providing gating strategies for myeloid and CD45− cells, we also wanted to investigate if the CITE-seq data would allow us to develop similar strategies for the lymphoid populations. To examine this, we re-clustered the T cells, NK cells, ILCs, B cells, and pDCs from Figure 1B identifying 12 distinct populations (Figure S5A). Analysis of the DEGs highlighted that while B cells, NK cells, ILC1s, and pDCs were distinct populations, there was considerable overlap between the transcriptomic profiles of the T cells (Figure S5B; Table S5). However, by analyzing the DEPs we were able to define distinct subsets including naive CD4 and CD8 T cells, TRegs, TH1s, CTLs, and TH17s (Figures S5A and S5C; Table S5). Moreover, we were able to design a gating strategy to isolate the distinct populations (Figure S5D).

9To determine the degree of conservation between the mac subsets and their different microenvironmental niches between the mouse and the human liver, we next generated a proteogenomic atlas of the human liver using sc/snRNA-seq and CITE-seq on 19 liver biopsies (Figures 4A, 4B, S6A, and S6B; Tables S1 and S5). Of these, most were histologically healthy with only 5 patients showing >10% hepatic steatosis in the absence of any significant fibrosis (Table S6). Cellular proportions varied according to the isolation technique used, and while there was some variability between patients, this was not linked to the surgery (Figures S6C–S6E). Further confirming the lack of fibrosis in the steatotic livers, we did not detect any increase in CTLs (Figures S5E–S5I; Table S5), which have been shown to correlate with non-alcoholic steatohepatitis (NASH). A significant increase in CTLs was detected in the setting of murine NASH induced by feeding a western diet (WD) for up to 36 weeks, demonstrating that this is not due a limitation in detecting these differences using CITE-seq (Figures S5J–S5M; Table S7). As Visium reliably located murine hepatic cells, we used this to locate the cells of the human liver in 4 biopsies (Figure 4C). However, the Visium spots from patients with >10% steatosis were found to cluster separately from the healthy samples (<10% steatosis; Figures 4D and 4E). We therefore used the healthy samples to calculate a baseline zonation and then transferred this trajectory onto the steatotic samples (Figures 4F and S6F). This identified the steatosis to be predominantly present in regions expressing peri-central zonation genes like CYP2E1. This zonation pattern was further validated using Molecular Cartography (Figure S6G). This fits with previous clinical studies demonstrating peri-central steatosis to be most common in non-alcoholic fatty liver disease (NAFLD) patients, especially in early disease. However, as these peri-central regions are larger than in the healthy controls, it could also imply that the presence of steatosis alters expression of the zonated hepatocyte genes, but this remains to be tested. Notably, the overall cellular distribution was not impacted by the presence of steatosis, although neutrophils and monocytes and monocyte-derived cells were preferentially localized peri-centrally in the steatotic patients correlating with the presence of steatosis (Figure 4G). MICS 100-plex protein analysis further validated the cellular distributions predicted by Visium and confirmed the increase in neutrophils in the steatotic livers (Figures S6H and S6I).

10To date, no validated markers of bona fide human KCs have been described. Explaining the difficulty to accurately define human KCs, we found monocytes and macs formed a single continuum in the human sc/snRNA-seq data, preventing a simple definition of human KCs (Figure 4B). Notably, a similar continuum from monocyte to KCs was also observed in the NASH murine liver (Figures S5J and S6J–S6N, Table S8). Consistent with our previous report, we observed both long-term resident Timd4-expressing KCs and recently recruited Timd4− monocyte-derived KCs (moKCs; Figures S6J–S6N). As the presence of such a continuum in the human liver suggests that there may also be monocyte contribution to the KC pool in the healthy human liver, we next zoomed in on myeloid cells to examine this, identifying 10 clusters (Figures 4H, S6O, and S6P; Table S2). To define the KCs, we examined expression of the top 25 murine KC genes by these clusters, which identified cluster10 to be the genuine human KCs (Figure S6Q). Unlike in mice, these were preferentially located in the mid zone (Figure S6R). Cluster9 also expressed many of these genes but lacked TIMD4 (Figure S6Q), suggesting that these cells may be recently recruited moKCs. The presence of moKCs in the liver is consistent with reports that host-derived macs are identified in transplanted donor livers and suggests that the KC population may be a mix of embryonic and monocyte-derived cells. Although not the case at mRNA level, VSIG4 was found to be the best human KC protein marker in the CITE-seq data, while FOLR2, CD163, and CD169 were also identified as useful markers of these cells for flow cytometry and confocal microscopy on frozen and paraffin sections (Figures 4I and S6S–S6X). Co-staining human livers for VSIG4 protein and KC-specific CD5L mRNA and MICS 100-plex protein analysis also confirmed the mid-zonal localization of KCs (Figures 4J and S6U). To assess if KC identity was further conserved in evolution, we profiled macaque, pig, hamster, chicken, and zebrafish livers (Figure 4K). We identified the KCs in an unbiased manner by mapping the conserved human-mouse KC signature onto the datasets (Figures 4K and S7A–S7C). We then examined the main features of each KC population identified (Figures S7D–S7H; Table S9). A strong overlap in transcriptomes across species was observed likely due to the conserved expression of core KC transcription factors (Figure S7I). However, each species also harbored a number of unique KC genes (Figure S7J; Table S9). VSIG4 protein expression was also conserved in pig and macaque KCs (Figures S7K–S7M). Similarly, we were also able to identify most of the other hepatic cells across species on the basis of conserved genes (Figures S7N and S7O). cDC2s were the main exception to this, as specific cDC2 marker genes were not conserved across all species (Figure S7O).

11Alongside KCs, we also identified distinct clusters of macs in the human myeloid cells (Figure 4H). To better understand the nature of these clusters we performed confocal microscopy to examine the specific locations of CD68+VSIG4− macs in the human liver (Figure S8A). This identified CD68+VSIG4− macs in the liver capsule, in close proximity to central and PVs as well as at bile ducts (BDs) (Figures 5A–5C and S8A). Similar populations were also observed at the PVs and CVs and at the BDs in the healthy macaque liver (Figure S7M). Examination of the scRNA-seq data and comparison with murine signatures identified immature and mature LAMs, with immature LAMs expressing some monocyte genes (Figures 4H, 5D, S8B, and S8C). Although recently suggested to be specific to fibrotic human livers, we identified LAMs in all patients profiled with scRNA-seq, but there was a trend toward increased proportions of LAMs in the livers with >10% steatosis (Figure S8D) consistent with the increased population of LAMs in murine NAFLD (Figures S6J–S6N). As in the healthy mouse, Visium identified human LAMs in portal zones of non-steatotic livers. However, in steatotic human livers, LAMs were primarily located peri-centrally, in zones with steatosis (Figure 5E), suggesting that monocytes are recruited to distinct locations in the healthy and obese liver where they then differentiate into LAMs. This altered location of LAMs was further validated by confocal microscopy and Molecular Cartography (Figures 5F–5H). However, this analysis did not identify any capsule macs, suggesting that these cells may be absent from our UMAP, likely as a result of the small amount of capsule tissue on a biopsy. The Mac1 population expressing IGSF21 was present in very low numbers throughput the tissue (Figures 5G and 5H). This coupled with their similar transcriptomic profile to moKCs could suggest that these are moKC precursors as observed in the mouse, but this requires further study. Focusing on the LAMs, the change in their location in the steatotic human liver was also observed in the murine NAFLD model. Here, LAMs were found across portal, periportal, and mid zones (Figures 5I, S8E, and S8F), fitting with the presence of steatosis in these regions and consistent with our previous report. Comparison of DEGs between LAMs in standard diet (SD) and WD-fed mice identified that LAMs had a more mature phenotype in WD-fed mice, downregulating their expression of some monocyte genes and increasing their expression of prototypical mac markers, consistent with the presence of both immature and mature LAMs in the human liver (Figure 5J). Fitting with a more mature phenotype, WD-derived LAMs also expressed lower levels of Il1b, Tnf, and Il10 compared with SD LAMs (Figure 5K). While further studies are required to assess the precise functions of these cells in NAFLD, this could further suggest a protective rather than a pathogenic role for LAMs.

12Given the altered localization of KCs in the healthy murine versus human liver and LAMs in the healthy versus steatotic liver of mice and humans, we next sought to investigate how the cells of the mac niches differed in these settings. Analysis of CD45− cells in the human liver identified similar (sub)populations of ECs, SCs, and hepatocytes as observed in the healthy mouse liver (Figures 6A and 6B; Table S3). However, while the 100-plex MICS analysis detected a second subset of portal fibroblasts expressing DESMIN (Figure S6H), we did not detect these cells in our UMAP, likely due to difficulties detecting DESMIN with snRNA-seq. No significant differences were found in terms of localization of the identified CD45− cells that could explain the altered location of KCs compared with the mouse (Figures 6C and 6D). With this in mind, we next utilized NicheNet to examine potential ligand-receptor pairs between KCs and CD45− cells present only in the human that could regulate KC location. This analysis identified CCL23 and CCL14 expression by LSECs and CCL16 by hepatocytes that binds to CCR1 expressed by human KCs but not murine KCs (Figure 6E). Crucially, these ligands were preferentially located peri-centrally (Figure 6F) and thus represent interesting targets for further study regarding the potential regulation of KC location in the human liver.

13We next aimed to investigate the role of CD45− cells in the regulation of LAM location in the healthy versus steatotic liver. Given the low number of steatotic human samples in our sc/snRNA-seq analysis and the variation between patients, we turned to the murine NASH model to investigate this. Zooming in on SCs (Figure 6G; Table S10), we identified that there was an increase in fibroblasts in the WD-fed mice. There was also a considerable overlap between fibroblasts and HSCs in terms of their gene expression profiles suggesting there may be a differentiation trajectory between the HSCs and the fibroblasts, although this remains to be validated in vivo. We also noted a sub-population of fibroblasts expressing Ccl2, a known ligand for CCR2 which is expressed on monocytes recruited to the liver during NAFLD, as well as Cd44 and Vcam1, two genes involved in monocyte recruitment and adhesion. Both fibroblast populations were enriched in the periportal steatotic regions (Figure 6I), in which we also observe enrichment of the LAMs (Figure 5I). We also found high CCL2 and Cd44 expression in human fibroblasts (Table S3), potentially recruiting bile-duct LAMs. Together, this suggests that fibroblasts may play an important role in recruiting LAMs.

14Having identified a potential role for the mac niche cells in regulating mac subset location, we next sought to determine the involvement of these cells in regulating mac phenotypes. To assess the roles of conserved cell-cell interactions in driving mac heterogeneity across species, we performed a differential NicheNet analysis between the distinct hepatic macs and the CD45− cells present in their respective niches focusing on ligands and receptors conserved in both human and mouse. This revealed very few specific ligand-receptor pairs for LAMs (Figure S8G; data not shown), hinting that local factors such as metabolites rather than unique cell-cell interactions may drive the LAM phenotype. Indeed, this would be consistent with the presence of LAM-like cells in multiple tissues including the obese adipose tissue, the brain, the lung, and the heart. In line with this, BM monocytes cultured with acetylated low-density lipoprotein expressed LAM-associated genes (Figure 7A), demonstrating a dominant role for lipids in inducing the LAM phenotype. Conversely, for KCs, we found multiple ligand-receptor pairs to be conserved between human and mouse (Figures 7B and S8G). One of these, an activin receptor-like kinase (ALK1)-bone morphogenic protein (BMP)9/10 circuit between KCs (ALK1; encoded by Acvrl1) and stellate cells (BMP)9/10 encoded by Gdf2/Bmp10 respectively) was found to be conserved in all 7 species and was predicted to control the expression of a number of the conserved KC genes (Figures 7B, 7C, and S8H). To validate a role for this axis in KCs, we generated Fcgr1-Cre × Acvrl1fl/fl mice, eliminating ALK1 specifically from CD64-expressing macs. This led to an almost complete loss of VSIG4+ KCs (Figures 8D–8F), demonstrating that evolutionarily conserved ALK1 signaling is crucial for KCs. To determine if ALK1 is required for KC maintenance we generated Clec4f-Cre × Acvrl1fl/fl mice, eliminating ALK1 only from differentiated KCs. This revealed a relatively similar phenotype than observed in the Fcgr1-Cre mice, suggesting that ALK1 is also required for KC maintenance (Figure S8I). To directly test the need for ALK1 in KC development, we generated BM chimeras whereby CD45.1 Clec4f-Dtr mice were irradiated with their livers shielded to avoid any radio damage. These mice were then reconstituted with either CD45.2 Acvrl1fl/fl or Fcgr1-CrexAcvrl1fl/fl BM. 4 weeks later, mice were given a single i.p. injection of DT to deplete the KCs and 7 or 13 days thereafter chimerism within the KC population was determined (Figure 7G). We observed that, already by day 7, KO BM-derived cells were almost completely outcompeted by WT BM-derived cells within the KC pool. This was not observed in monocytes but a similar pattern was observed in VSIG4− macs demonstrating that ALK1 is critically required for early KC development (Figure 7H). Finally, while NicheNet predicts that BMP9/10 from stellate cells would signal through ALK1 to induce KC development and maintenance, a prediction in line with our previous NicheNet analysis, another recent study has suggested that transforming growth factor (TGF)β signaling would be important for KCs. This prediction was made on the basis of SMAD4 signaling, which is a common downstream effector of both TGF-β Receptor and ALK1-induced signaling. As TGF-β was also found to participate in an ALK1-TGF-β receptor containing signaling complex, we thus examined whether TGF-β signaling is important for KCs. To this end we utilized an ALK1-Fc trap or TGF-β type II receptor (TGFβRII)-Fc trap alongside appropriate isotype controls. ALK1-Fc selectively sequesters BMP9/10, while TGFβRII-Fc selectively interferes with TGF-β1/TGF-β3 receptor binding. Clec4f-Dtr mice were depleted of KCs and simultaneously treated with either receptor-Fc traps or isotype controls and KC development was examined 7 days later (Figure 7I). In line with the chimera study, treatment with ALK1-Fc significantly abrogated KC development; however, blocking TGFβ signaling had only a minor effect on the proportion of VSIG4− macs (Figure 7J). Altogether this validates the NicheNet prediction and demonstrates that the evolutionarily conserved ALK1-BMP9/10 axis is crucial for the development and maintenance of KCs.

15To generate a practical cellular atlas of any human tissue and unravel the cell-cell circuits essential for the identities of cells inhabiting that tissue, four key pieces of information are required: (1) an inventory of all cells present, (2) the location of the different cells within the tissue to identify interactions between neighboring cells, (3) an alignment between the human and animal models allowing for any predicted cell-cell interactions to be perturbed, and (4) the identification of reliable antibody-based panels for the efficient screening of different patients and/or transgenic animals. Here, by integrating single-cell and spatial transcriptomic and proteomic data, we provide these 4 pieces of information for the liver and uncover evolutionarily conserved microenvironmental circuits controlling the development of hepatic macs.

16Unraveling the spatial localization of all hepatic cells, we identify LAMs around the bile ducts in the healthy mouse, human, and macaque liver. However, when steatosis is present, LAMs are preferentially recruited to the steatotic regions of the liver. This spatial information at least partially invalidates the hypothesis that LAM identity is specifically induced by fibrotic SCs. Rather, our data suggest that LAMs are induced by local lipid exposure. We also provide an alignment of the liver atlas across seven species. This reveals the conserved and unique transcriptomic programs of steady-state KCs and uncovers the spatially restricted and conserved ligand-receptor pairs between KCs and the cells constituting their niche. Underlining the need to first characterize the healthy tissue before attempting to understand how disease perturbs the cells, we identify the DLL-NOTCH interaction to be an evolutionarily conserved cross-talk between homeostatic LSECs and KCs and therefore not unique to hepatocellular carcinoma or fibrosis, as proposed. Similarly, we find that FOLR2 expression is not specific to tumor-associated hepatic macs but is expressed by KCs in the healthy mouse and human liver. Finally, we apply a proteogenomic pipeline starting from broad oligo-conjugated antibody panels for both single-cell and spatial profiling. This is crucial as transcriptomic profiling does not always correspond with the ability to detect proteins by flow cytometry or microscopy. By screening broadly, we identify the best surface markers for the isolation and localization of hepatic macs and their respective niche cells. This allows both the validation of the spatial location at the single-cell level, and the efficient screening of transgenic mouse models for the loss of KCs. Characterization of Fcgr1-CrexAcvrl1fl/fl mice using our defined panel readily demonstrates the cruciality of the ALK1-BMP9/10 axis in KC development emphasizing that mac-stromal-cell cross-talk goes much further than the exchange of growth factors. Moving forward, applying these relatively cheap antibody panels to large patient cohorts or multiple transgenic mouse models should enable any perturbations disturbing liver homeostasis to be efficiently identified.

17The current study has two main limitations. First, the analysis of the human liver remains restricted by the number of human patients included (19 patients for the sc/snRNA-seq, 4 patients for Visium, and 15 patients by microscopy). While this study provides markers to cheaply and efficiently screen large patient cohorts allowing this analysis to be extended, given the heterogeneity between patients, multiple studies will need to be integrated in a single analysis if we are to be able to interrogate transcriptomic differences in a particular cell subset related to age, sex, ethnicity, or pathological parameters. Second, this study highlights that more research will be needed to fully characterize human SCs. As we could only retrieve these cells through snRNA-seq this means we have not been able to identify good surface markers through CITE-seq. Additionally, for small populations of fibroblasts we did not recover enough nuclei to truly probe their heterogeneity. Better isolation protocols are required to retrieve and enrich these cells from the human liver in order to run a broad panel of CITE-seq markers to identify the different subsets and their corresponding surface markers.

18Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Martin Guilliams (martin.guilliams@irc.vib-ugent.be).

19This study did not generate new unique reagents.

20WT C57Bl/6J mice (Janvier) were used for this study. Male and Female mice were used for all experiments, with the exception of NAFLD experiments were only male mice were used. For NAFLD experiments, mice were put on the SD or WD at 5 weeks of age and sacrificed after either 24 or 36 weeks on the diet as indicated. Fcgr1-Cre mice (Scott et al., 2018) were obtained from Prof. Bernard Malissen, CIML, Marseille and Clec4f-Cre mice (Scott et al., 2018) were crossed with Acvrl1fl/fl mice (Park et al., 2008) obtained from Paul Oh, Barrow Neurological Institute, Florida, USA. Clec4f-Dtr mice (Scott et al., 2016) were crossed to CD45.1 mice (Janvier) for KC depletion and development experiments All mice were used between 6 and 12 weeks of age unless otherwise stated. All mice maintained at the VIB (Ghent University) under specific pathogen free conditions. All animals were randomly allocated to experimental groups. All experiments were performed in accordance with the ethical committee of the Faculty of Science, UGent and VIB.

21Piglets (female, 10 weeks old) were purchased at a local farm and transported to the animal facilities of the Faculty of Veterinary Medicine. The animals were housed in isolation units as blood donors and had access to feed and water ad libitum. At 30 weeks of age, the animals were euthanized by intravenous injection of sodium pentobarbital 20% (60mg/2.5kg) and livers were collected. The animal study was reviewed and approved by the Ethical Committee of the Faculty of Veterinary Medicine (EC2018/55).

22Study animals were clinically healthy Leghorn hens of approximately 58 weeks old collected from a commercial farm. The hens were housed at the Faculty of Veterinary Medicine according to acceptable welfare standards and were observed at least twice daily for health problems. Feed and water was offered ad libitum. The chickens were euthanized through intravenous injection (in the wing) with sodium pentobarbital. The EC approval number of this trial was EC2019/015.

23Male Cynomolgus macaques (≥2 years) were sourced from China and supplied by Guangzhou Xiangguan Biotech Co., Ltd and confirmed healthy before being assigned to the study. Animal handling, husbandry and euthanasia was performed by WuxiAppTec Co., Ltd., China according to local ethical guidelines (AAALAC accredited 2010). Study animals consists of 4 groups, orally dosed once with a Janssen proprietary immune modulator or vehicle control. Vehicle control animals only were used in this study. Liver tissue samples were snap-frozen immediately after euthanasia. Samples were thawed once before shipping to Ghent for snRNA-seq analysis.

24Female syrian hamsters (Janvier) were housed per one or two in ventilated isolator cages at a temperature of 21°C, humidity of 55% and 12:12 dark/light cycles, with access to food and water ad libitum and cage enrichment. All hamsters had SPF status at arrival and manipulations were performed in a laminar flow cabinet. Housing conditions and experimental procedures were approved by the ethical committee of KU Leuven (license P015-2020). Animals were euthanized at 6-8 weeks of age by intraperitoneal injection of 200 mg/mL sodium pentobarbital and livers were collected for analysis.

25Zebrafish were maintained under standard conditions, according to FELASA guidelines (Alestrom et al., 2019). All experimental procedures were approved by the ethical committee for animal welfare (CEBEA) from the ULB (Université Libre de Bruxelles) (Protocol #594N). The following transgenic lines at 6 months of age were used: Tg(mpeg1:EGFP)gl22; Tg(kdrl:Cre)s89; Tg(actb2:loxP-STOP-loxP-DsRedexpress)sd5 enabling macs to be sorted for sequencing as DsRed, GFP double positive cells.

26Patient studies were run in collaboration with Ghent University Hospital. Liver biopsies (1–2mm3) were isolated with informed consent from patients undergoing cholecystectomy or gastric bypass. In addition, liver biopsies were isolated from healthy adjacent tissue removed during liver resection due to colorectal cancer metastasis. In most cases, a second biopsy was also taken to evaluate liver histology. A full overview of all patient samples used in this study can be found in Table S6. Paraffin-embedded human liver samples were obtained through collaboration with Dr. Jan Lerut (Université Catholique de Louvain, UCL). All studies were performed in accordance with the ethical committee of the Ghent University Hospital (study numbers: 2015/1334 and 2017/0539).

27Liver cells were isolated by either ex vivo digestion (all species, except zebrafish) or in vivo liver perfusion (mice only) and digestion as described previously. Briefly, for ex vivo digestion, livers were isolated, cut into small pieces and incubated with 1mg/ml Collagenase A and 10U/ml DNAse at 37°C for 20 mins with shaking. For in vivo digestion, after retrograde cannulation, livers were perfused for 1-2mins with an EGTA-containing solution, followed by a 5min (6ml/min) perfusion with 0.2mg/ml collagenase A. Livers were then removed, minced and incubated for 20mins with 0.4mg/ml collagenase A and 10U/ml DNase at 37°C. All subsequent procedures were performed at 4°C. Samples were filtered over a 100μm mesh filter and red blood cells were lysed. Samples were again filtered over a 40μm mesh filter. At this point in vivo digestion samples only were subjected to two centrifugation steps of 1 min at 50g to isolate hepatocytes. Remaining liver cells (leukocytes, LSECs and HSCs; in vivo protocol) and total cells from the ex vivo digests were centrifuged at 400g for 5mins before proceeding to antibody staining for flow cytometry. A combination of Collagenase A and DNase were used to digest livers in both protocols to minimize cleavage of surface epitopes.

28Dissected livers from 6 months old transgenic zebrafish were triturated and treated with Liberase TM at 33°C for 20 min. Cells were then filtered through 40μm nylon mesh and washed with 2% FBS in PBS by centrifugation. Sytox Red was then added to the samples at a final concentration of 5nM to exclude nonviable cells before proceeding to flow cytometry. DsRed+GFP+ cells were then FACS-purified.

29Nuclei were isolated from snap frozen liver tissue with a sucrose gradient as previously described. Briefly, frozen liver tissue is homogenized using Kimble Douncer grinder set in 1ml homogenization buffer with RNAse inhibitors. Homogenised tissue is then is then subjected to density gradient (29% cushion – Optiprep) ultracentrifugation (7700rpm, 4°C, 30 mins). After resuspension, nuclei are stained with DAPI and intact nuclei were FACS-purified from remaining debris.

30To induce NAFLD and NASH, mice were fed a western diet (WD) high in fat, sugar and cholesterol for 24 or 36 weeks as described previously (Remmerie et al., 2020). This consisted of 58% fat, 1% cholesterol (Research Diets; D09061703i) and drinking water was supplemented with 23.1g/L fructose (MPBio) and 18.9g/L sucrose (VWR). Control mice were fed a standard diet with 11kcal% fat with corn starch (D12328i; Research Diets).

31Cells were pre-incubated with 2.4G2 antibody (Bioceros) to block Fc receptors and stained with appropriate antibodies at 4°C in the dark for 30-45 minutes. Cell viability was assessed using Fixable Viability dyes (eFluor780 or eFluor506; Thermo Fischer) and cell suspensions were analyzed with a BD FACSymphony or purified using a BD Symphony S6, BD FACSAria II or III. Nuclei were sorted on basis of DAPI positivity and size. Analysis was performed with FlowJo software (BD). Intracellular staining for CD207 was performed by fixing and permeabilizing extracellularly stained cells according to the manufacturer’s instructions using the FoxP3 Fixation/Permeabilization Kit (Thermo Fischer).

32Confocal staining was performed as described previously. Immediately after sacrificing mice with CO2, inferior vena cava were cannulated and livers were perfused (4 mL/min) with Antigenfix (Diapath) for 5 min at room temperature. After excision, 2-3 mm slices of livers were fixed further by immersion in Antigenfix for 1h at 4°C, washed in PBS, infused overnight in 34% sucrose and frozen in Tissue-Tek OCT compound (Sakura Finetek). 20 μm-thick slices cut on a cryostat (Microm HM 560, Thermo Scientific) were rehydrated in PBS for 5 min, permeabilized with 0,5% saponin and non-specific binding sites were blocked for 30 min with 2% bovine serum albumin, 1% fetal calf serum and 1% donkey or goat serum for 30 minutes. Tissue sections were labeled overnight at 4°C with primary antibodies followed by incubation for 1h at room temperature with secondary antibodies. When two rat antibodies were used on the same section, the directly conjugated rat antibody was incubated for 1h after staining with the unconjugated and anti-rat secondary and after an additional blocking step with 1% rat serum for 30 minutes. Slides were mounted in ProLong Diamond, imaged with a Zeiss LSM780 confocal microscope (Carl Zeiss, Oberkochen, Germany) with spectral detector and using spectral unmixing and analyzed using ImageJ and QuPath software.

33Experiments were performed using the RNAScope Multiplex Fluorescent V2 Assay kit (ACDBio 323100). Probes targeting intronic regions for Hs-Cd5l (ACDBio 850511), Mfa-Cd5l (ACDBio 873211), Mm-Cd5l (ACDBio 573271), Mm-Flt3 (ACDBio 487861), Mm-Xcr1 (ACDBio 562371), Mm-Mafb (ACDBio 438531) and Mm-Mgl2-O1 (ACDBio 822901) were custom-designed and synthesized. They were then labelled with TSA opal 520 (PerkinElmer FP1487001KT), TSA opal 540 (PerkinElmer FP1494001KT), TSA opal 570 (PerkinElmer FP1488001KT), TSA opal 620 (PerkinElmer FP1495001KT) or TSA opal 650 (PerkinElmer FP1496001KT). Tissues were fixed for 16 hours in AntigenFix (Diapath P0016), dehydrated and embedded in OCT as described above. Slices were pre-treated with hydrogen peroxide for 10 min and protease III for 20 min. The recommended Antigen retrieval step was not performed in order to preserve epitope integrity. Probes were hybridized and amplified according to the manufacturer’s instructions. Slides were then stained for protein markers as described above.

34Mice were euthanized by means of carbon dioxide (CO2) overdose. The liver was excised and consequently trimmed, on ice, to smaller tissue pieces fitting the 10X Visium capture area. Trimmed tissue pieces were embedded in Tissue-Tek® O.C.T.™ Compound (Sakura) and snap frozen in isopentane (Sigma) chilled by liquid nitrogen. Embedded tissue pieces where stored at -80°C until cryosectioning.

35A 10X Visium Spatial Gene expression slide was placed in the cryostat (Cryostar NX70 Thermo Fisher) 30 minutes prior to cutting. 10 μm sections where cut and placed within the capture area. Single 10X Visium Spatial Gene expression slides were stored in an airtight container at -80°C until further processing.

3610X Visium cDNA libraries were generated according the manufacturer’s instructions. In short: Tissue sections where fixed in chilled Methanol. A H&E staining was performed to assess tissue morphology and quality. Tissue was lysed and reverse transcription was performed followed by second strand synthesis and cDNA denaturation. cDNA was transferred to a PCR tube and concentration was determined by qPCR. Spatially barcoded, full length cDNA was amplified by PCR. Indexed sequencing libraries where generated via End Repair, A-tailing, adaptor ligation and sample index PCR. Full length cDNA and indexed sequencing libraries were analyzed using the Qubit 4 fluorometer (Thermo Fisher) and Agilent 2100 BioAnalyzer.

37Liver slices were prepared as described above for the classical Visium protocol. Slices were dried for 1 min at 37°C and subsequently fixed using 1% paraformaldehyde in PBS. Next, slices were blocked for 30 min (2% BSA, 0.1ug/ul Salmon Sperm, 0.5% Saponin, 1 U/μl protector RNase inhibitor (Roche) in 3X SSC) and incubated with the oligo-conjugated antibody staining mix (2% BSA, 0.1μg/μl Salmon Sperm, 0.5% Saponin, 1 U/μl protector RNase inhibitor, 10uM polyT-blocking oligo (TTTTTTTTTTTTTTTTTTTTTTT∗T∗T∗/3InvdT/), in 3X SSC) for 1h at 4°C. Slides were mounted (90% glycerol, 1 U/μl protector RNase inhibitor) and imaged on Zeiss Axioscan Z1 at 20X magnification. Samples were then processed for a transcriptomic experiment as per manufacturer’s instructions (Visium, 10X Genomics) with modifications to also capture antibody tags. In short, tissue was permeabilized using Tissue Removal Enzyme (Tissue Optimization kit, 10x Genomics) for 9 minutes, as determined by a tissue optimization experiment (10X Genomics, Visium Spatial Tissue Optimization). After reverse transcription, 2 μl of 100 μM FB additive primer (CCTTGGCACCCGAGAATT∗C∗C∗A) per sample was added to the second strand synthesis mix. During cDNA amplification 1 μl of 0,2 μM FB additive primer (CCTTGGCACCCGAGAATT∗C∗C∗A) was added. After cDNA amplification, antibody products and mRNA derived cDNA were separated by 0.6X SPRI select. The purified full-length cDNA fraction was quantified by qRT-PCR using KAPA SYBR FAST-qPCR kit on a PCR amplification and detection instrument. After enzymatic fragmentation indexed sequencing libraries were generated via End Repair, A-Tailing, adaptor ligation and sample index PCR. The supernatant containing antibody product was cleaned up by two rounds of 1.9X SPRI select. Next, 45 μl of the purified antibody fraction was amplified with a 96 deep well reaction module: 95°C for 3 min; cycled 8 times: 95°C for 20 s, 60°C for 30 s, and 72°C for 20 s; 72°C for 5 min; end at 4°C. ADT libraries were purified once more with 1.6X SPRI select. Full length cDNA, indexed cDNA libraries and antibody libraries were analyzed using the Qubit 4 fluorometer (Thermo Fisher) and Agilent 2100 Bioanalyzer. The separation of the cDNA and ADT libraries were performed according to the manufacturer’s instructions (10X genomics).

38The MACSima™ Imaging System is a fully automated instrument combining liquid handling with widefield microscopy for cyclic immunofluorescence imaging. In brief, staining cycles consisted of the following automated steps: immunofluorescent staining, sample washing, multi-field imaging, and signal erasure (photobleaching or REAlease).

39Cryosectioned slices on slides were taken out of the -80°C storage and the appropriate MACSWell™ imaging frame was mounted immediately on the slide. An appropriate volume of ice-cold 4% PFA solution was added (according to the MACSWell™ imaging frames datasheet) and incubated for 10 minutes at room temperature. The slide was washed three times with MACSima Running Buffer. After washing the appropriate initial sample volume of MACSima Running Buffer was added (according to the MACSWell™ imaging frames datasheet). Right before the start of the MACSima ™ instrument a DAPI pre-staining was performed: the MACSima Running Buffer was removed from the sample to be analysed and stained for 10 min with a 1:10 dilution of a DAPI staining solution (volume depends on working volume for the different MACSwell™ formats, see datasheet). The DAPI staining solution was removed and 3 washing steps were performed (MACSima Running Buffer). Finally, the initial sample volume of MACSima Running Buffer was added. Details of the antibodies used can be found in the key resources table.

40Liver was frozen and sectioned as described above for Visium analysis and liver slices were placed within capture areas on Resolve BioScience slides. Samples were then sent to Resolve BioSciences on dry ice for analysis. Upon arrival, tissue sections were thawed and fixed with 4% v/v Formaldehyde (Sigma-Aldrich F8775) in 1x PBS for 30 min at 4 °C. After fixation, sections were washed twice in 1x PBS for two min, followed by one min washes in 50% Ethanol and 70% Ethanol at room temperature. Fixed samples were used for Molecular Cartography™ (100-plex combinatorial single molecule fluorescence in-situ hybridization) according to the manufacturer’s instructions (protocol 3.0; available for download from Resolve’s website to registered users), starting with the aspiration of ethanol and the addition of buffer BST1 (step 6 and 7 of the tissue priming protocol). Briefly, tissues were primed followed by overnight hybridization of all probes specific for the target genes (see below for probe design details and target list). Samples were washed the next day to remove excess probes and fluorescently tagged in a two-step color development process. Regions of interest were imaged as described below and fluorescent signals removed during decolorization. Color development, imaging and decolorization were repeated for multiple cycles to build a unique combinatorial code for every target gene that was derived from raw images as described below.

41The probes for 100 genes were designed using Resolve’s proprietary design algorithm. Briefly, the probe-design was performed at the gene-level. For every targeted gene all full-length protein-coding transcript sequences from the ENSEMBL database were used as design targets if the isoform had the GENCODE annotation tag ‘basic’. To speed up the process, the calculation of computationally expensive parts, especially the off-target searches, the selection of probe sequences was not performed randomly, but limited to sequences with high success rates. To filter highly repetitive regions, the abundance of k-mers was obtained from the background transcriptome using Jellyfish. Every target sequence was scanned once for all k-mers, and those regions with rare k-mers were preferred as seeds for full probe design. A probe candidate was generated by extending a seed sequence until a certain target stability was reached. A set of simple rules was applied to discard sequences that were found experimentally to cause problems.

42After these fast screens, every kept probe candidate was mapped to the background transcriptome using ThermonucleotideBLAST and probes with stable off-target hits were discarded. Specific probes were then scored based on the number of on-target matches (isoforms), which were weighted by their associated APPRIS level, favoring principal isoforms over others. A bonus was added if the binding-site was inside the protein-coding region. From the pool of accepted probes, the final set was composed by greedily picking the highest scoring probes. Probe details are included in the key resources table.

43Samples were imaged on a Zeiss Celldiscoverer 7, using the 50x Plan Apochromat water immersion objective with an NA of 1.2 and the 0.5x magnification changer, resulting in a 25x final magnification. Standard CD7 LED excitation light source, filters, and dichroic mirrors were used together with customized emission filters optimized for detecting specific signals. Excitation time per image was 1000 ms for each channel (DAPI was 20 ms). A z-stack was taken at each region with a distance per z-slice according to the Nyquist-Shannon sampling theorem. The custom CD7 CMOS camera (Zeiss Axiocam Mono 712, 3.45 μm pixel size) was used.

44For each region, a z-stack per fluorescent color (two colors) was imaged per imaging round. A total of 8 imaging rounds were done for each position, resulting in 16 z-stacks per region. The completely automated imaging process per round (including water immersion generation and precise relocation of regions to image in all three dimensions) was realized by a custom python script using the scripting API of the Zeiss ZEN software (Open application development).

45The algorithms for spot segmentation were written in Java and are based on the ImageJ library functionalities. Only the iterative closest point algorithm is written in C++ based on the libpointmatcher library (https://github.com/ethz-asl/libpointmatcher).

46As a first step all images were corrected for background fluorescence. A target value for the allowed number of maxima was determined based upon the area of the slice in μm2 multiplied by the factor 0.5. This factor was empirically optimized. The brightest maxima per plane were determined, based upon an empirically optimized threshold. The number and location of the respective maxima was stored. This procedure was done for every image slice independently. Maxima that did not have a neighboring maximum in an adjacent slice (called z-group) were excluded. The resulting maxima list was further filtered in an iterative loop by adjusting the allowed thresholds for (Babs-Bback) and (Bperi-Bback) to reach a feature target value (Babs: absolute brightness, Bback: local background, Bperi: background of periphery within 1 pixel). This feature target values were based upon the volume of the 3D-image. Only maxima still in a z-group of at least 2 after filtering were passing the filter step. Each z-group was counted as one hit. The members of the z-groups with the highest absolute brightness were used as features and written to a file. They resemble a 3D-point cloud.

47To align the raw data images from different imaging rounds, images had to be corrected. To do so the extracted feature point clouds were used to find the transformation matrices. For this purpose, an iterative closest point cloud algorithm was used to minimize the error between two point-clouds. The point clouds of each round were aligned to the point cloud of round one (reference point cloud). The corresponding point clouds were stored for downstream processes. Based upon the transformation matrices the corresponding images were processed by a rigid transformation using trilinear interpolation.

48The aligned images were used to create a profile for each pixel consisting of 16 values (16 images from two color channels in 8 imaging rounds). The pixel profiles were filtered for variance from zero normalized by total brightness of all pixels in the profile. Matched pixel profiles with the highest score were assigned as an ID to the pixel.

49Pixels with neighbors having the same ID were grouped. The pixel groups were filtered by group size, number of direct adjacent pixels in group, number of dimensions with size of two pixels. The local 3D-maxima of the groups were determined as potential final transcript locations. Maxima were filtered by number of maxima in the raw data images where a maximum was expected. Remaining maxima were further evaluated by the fit to the corresponding code. The remaining maxima were written to the results file and considered to resemble transcripts of the corresponding gene. The ratio of signals matching to codes used in the experiment and signals matching to codes not used in the experiment were used as estimation for specificity (false positives).

50Final image analysis was performed in ImageJ using genexyz Polylux tool plugin from Resolve BioSciences to examine specific Molecular Cartography™ signals.

5140000-160000 cells of interest from livers of the different species were FACS-purified and pelleted by centrifugation at 400g for 5 mins. To ensure sufficient numbers of all cell types were present in our analyses, depending on the sample distinct populations of cells including Live CD45+, Live CD45-, Hepatocytes, Myeloid cells and Stromal cells were FACS-purified. When CITE-seq was to be performed, cells were then stained with 2.4G2 antibody to block Fc receptors and CITE-seq antibodies for 20mins at 4°C, before being washed in excess PBS with 2% FCS and 2mM EDTA. Antibody details are included in the key resources table. 40000-100000 nuclei were also FACS-purified based on DAPI expression. These were sorted into BSA coated tubes and pelleted by centrifuging for 3 mins at 400g and 5 mins at 600g sequentially. Cells/Nuclei were then resuspended in PBS with 0.04%BSA at ∼1000 cells/ml. Cell suspensions (target recovery of 8000-10000 cells) were loaded on a GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) to generate single-cell Gel Bead-in-Emulsions (GEMs). Single-cell RNA-Seq libraries were prepared using GemCode Single-Cell 3ʹGel Bead and Library Kit (10x Genomics, V2 and V3 technology) according to the manufacturer’s instructions. Briefly, GEM-RT was performed in a 96-Deep Well Reaction Module: 55°C for 45min, 85°C for 5 min; end at 4°C. After RT, GEMs were broken down and the cDNA was cleaned up with DynaBeads MyOne Silane Beads (Thermo Fisher Scientific, 37002D) and SPRIselect Reagent Kit (SPRI; Beckman Coulter; B23318). cDNA was amplified with 96-Deep Well Reaction Module: 98°C for 3 min; cycled 12 times: 98°C for 15s, 67°C for 20 s, and 72°C for 1 min; 72°C for 1 min; end at 4°C. Amplified cDNA product was cleaned up with SPRIselect Reagent Kit prior to enzymatic fragmentation. Indexed sequencing libraries were generated using the reagents in the GemCode Single-Cell 3ʹ Library Kit with the following intermediates: (1) end repair; (2) A-tailing; (3) adapter ligation; (4) post-ligation SPRIselect cleanup and (5) sample index PCR. Pre-fragmentation and post-sample index PCR samples were analyzed using the Agilent 2100 Bioanalyzer.

52RNA was extracted from 10000 sorted cells (gated using strategies shown) from livers of C57BL/6 mice using a RNeasy Plus micro kit (QIAGEN). Sensifast cDNA synthesis kit (Bioline) was used to transcribe total RNA to cDNA. Real-time RT-PCR using SensiFast SYBR No-Rox kit (Bioline) was performed to determine gene expression, therefore a PCR amplification and detection instrument LightCycler 480 (Roche) was used. Gene expression was normalized to β-actin gene expression. Primers used in the study can be found in the key resources table.

53sc/snRNA-seq libraries were loaded on an Illumina HiSeq or Illumina NovaSeq 6000 with sequencing settings recommended by 10X Genomics (26/8/0/98 – 2.1pM loading concentration, ADT and cDNA libraries were pooled in a 25:75 ratio). Visium sequencing libraries were loaded on an Illumina NovaSeq 6000 with sequencing settings recommended by 10X Genomics (28/10/10/75 – 2.1pM loading concentration). Sequencing was performed at the VIB Nucleomics Core (VIB, Leuven). The demultiplexing of the raw data was performed using CellRanger software (10x – version 3.1.0; cellranger mkfastq which wraps Illumina’s bcl2fastq). The reads obtained from the demultiplexing were used as the input for ‘cellranger count’ (CellRanger software), which aligns the reads to the mouse reference genome (mm10) using STAR and collapses to unique molecular identifier (UMI) counts. The result is a large digital expression matrix with cell barcodes as rows and gene identities as columns.

54To remove ambient RNA, the FastCAR R package (v0.1.0) with a contamination chance cutoff of 0.05 was run on the samples separately before merging them. The UMI cut off was determined individually for the different samples, using the CellRanger web_summary output plot (see GitHub). The Scater R package (v1.14.6) was used for the preprocessing of the data. The workflow to identify the outliers, based on 3 metrics (library size, number of expressed genes and mitochondrial proportion) described by the Marioni lab was followed. As a first step cells with a value x median absolute deviation (MADs) higher or lower than the median value for each metric were removed. This value was determined individually for the different datasets (see github). Secondly, the runPCA function (default parameters) of the Scater R package was used to generate a principal component analysis (PCA) plot. The outliers in this PCA plot were identified by the R package mvoutlier. By creating the Seurat object, genes that didn’t have an expression in at least 3 cells were removed. To normalize, scale and detecting the highly variable genes, the R package SCTransform (v0.2.1) was used. If batch correction (on sample level) was needed, the NormalizeData (log2 transformation), FindVariableFeatures and ScaleData functions of the Seurat R package (v3.1.2) were used in combination with the Harmony R package (v1.0). The Seurat pipeline was followed to find the clusters and create the UMAP plots. The number of principal components used for the clustering and the resolution were determined individually for the different datasets (see GitHub). On these initial UMAP plots we did multiple rounds of cleaning by removing proliferating and contaminating (e.g. doublets) cells. For non CITE-seq datasets the count data for the clean cells acquired by the previous steps were further processed with the scVI model (scvi Python package v0.6.7) (Lopez et al., 2018). Datasets including Cite-seq samples were further processed with the TotalVI model. The workflows described on scvi-tools.org were followed to generate new UMAPs, DEGs and DEPs. This information was further processed with the pheatmap R package (v1.0.12) to create heatmaps using the normalized values (denoised genes) calculated in the scVI/TotalVI workflow. The plots showing the expression of certain genes or proteins are created with the ggplot2 R package (v3.2.1) with a quantile cut off of 0.01.

55For mouse all the ABs from the whitelist (181 ABs) were loaded into TotalVI, while for the other species only the added ABs were loaded into TotalVI. For the ‘human liver-pool of techniques and patients’ we noticed that the batch correction (between samples) faced difficulties for the hepatocytes and stellate cells as the cells all originated from snRNA-Seq samples, while the other cell types originated from both snRNA-seq and scRNA-seq samples. To overcome this issue we randomly allocated 30% of the hepatocytes to scRNA-seq samples which were not CITE-seq samples. We did the same for 30% of the stellate cells.

56Heatmaps were made by scaling the normalized values (denoised values; calculated in the scVI/TotalVI workflow) using the scale_quantile function of the SCORPIUS R package (v1.0.7) and the pheatmap R package (v1.0.12). The plots showing the expression of certain genes or proteins were created based on the normalized values (denoised values) using a quantile cutoff of 0.99 and via either the ggplot2 R package (v3.2.1) or the scanpy.pl.umap function of the Scanpy Python package (v1.5.1).

57To find the conserved human and mouse KC markers we started by identifying the human KC markers. We mapped the annotation of the human myeloid UMAP on the human pool of techniques/patients UMAP to identify the real KCs in this last UMAP. The real KCs were identified as the top part of the mac cluster. Using this new annotation we then calculated the DE genes and DE proteins for each cluster. Some genes are listed as marker for multiple clusters, only for the cluster where the gene had the highest score (raw_normalized_mean1/raw_normalized_mean2∗lfc_mean), the gene was kept as marker. This way we found 110 potential human KC markers. We then created a heatmap of these 110 genes (using denoised gene values scaled between 0 and 1) and filtered this heatmap by removing the genes where the scaled normalized value was higher than 0.50 in more than 30% of the cells of a certain cell type other than KCs. Except for the macs, we only removed a gene when it had a scaled normalized value higher than 0.50 in more than 70% of the macs. After this filtering we ended up with 36 human KC markers. Next we converted these human gene symbols into MGI IDs via the BioMart tool on the HGNC website (https://biomart.genenames.org/martform/#!/default/HGNC?datasets=hgnc_gene_mart). We found a MGI ID for 30 genes. We then converted these MGI IDs into mouse gene symbols via the MGI webtool (http://www.informatics.jax.org/batch/).

58To identify the mouse KC markers we similarly mapped the annotation of the mouse myeloid UMAP on the mouse pool of techniques UMAP to identify the real KCs in this last UMAP. The real KCs matched with the mac cluster. Similarly as in human, the DE genes for each cluster was calculated and genes listed as marker for multiple clusters were dealt with in a similar way. This way we found 264 potential mouse KC markers. We then removed the genes that had a score (raw_normalized_mean1/raw_normalized_mean2∗lfc_mean) lower than 10 and ended up with 214 genes. We then created a heatmap of these 214 genes (using denoised gene values scaled between 0 and 1) and filtered this heatmap by removing the genes where the scaled normalized value was higher than 0.50 in more than 30% of the cells of a certain cell type other than KCs. After this filtering we ended up with 68 mouse KC markers. Next we converted these mouse gene symbols into MGI IDs via the MGI webtool (http://www.informatics.jax.org/batch/). We then converted these MGI IDs into human gene symbols via the BioMart tool on the HGNC website (https://biomart.genenames.org/martform/#!/default/HGNC?datasets=hgnc_gene_mart) and ended up with 60 genes.

59At this point we found 30 human KC markers and 60 mouse KC markers. In a next step, we only kept the human KC markers that we identified as a Highly Variable Gene (HVG) in the mouse pool of techniques UMAP (20 genes) and the mouse KC markers that were identified as HVGs in the human pool of the techniques UMAP (30 genes). We next put these 20 mouse KC markers in SingleCellSignatureExplorer to see where these genes are enriched in the mouse pool of techniques UMAP. In order to only get an enrichment in the KCs we decided to only use top 10 mouse KC markers (ordered on score), together with Slc40a1 and Hmox1. We then started to add the top human KC markers as long as we keep the enrichment solely in the KCs. This way we ended up with final list of 15 human-mouse conserved KC markers.

60We next converted these KC markers into the monkey, pig, chicken or zebrafish orthologs by looking up the human gene symbol on NCBI (https://www.ncbi.nlm.nih.gov/search/) and checking if there is an ortholog of the species of interest listed under the ‘Ortholog’ tab. The found orthologs were then used as input for the SingleCellSignatureExplorer tool.

61The protein normalized values (denoised values; calculated in the TotalVI workflow) were converted into an FCS file using the write.FCS function of the flowCore R package (v1.50.0).

62We first removed per sample all spots that were clear outliers compared to the location of the tissue. Each sample was then normalized individually using the SCTransform function of the Seurat R package (v3.2.3) with default parameters. All samples were then merged with the merge function of the Seurat R package (v3.2.3) with default parameters. Next, we determined the HVGs, created a PCA plot, performed clustering and created an UMAP plot as described in the spatial workflow available on the Seurat website (https://satijalab.org/seurat/articles/spatial_vignette.html). Clusters which showed high mitochondrial gene expression were removed. Spots located at the darker parts of the tissue were also removed as these parts are considered to be dead tissue or of bad quality.

63For modelling the cell type composition and zonation, spatial CITE-seq and transcriptomics data were analyzed using probabilistic graphical models, similar to what is used in tools such as cell2location and scVI. In brief, transcriptomics data was modelled as a NegativeBinomial distribution, parameterized with a mean and dispersion, the latter optimized as a free parameter for each gene. Visium Highly Multiplexed Protein data was modelled as a mixture of NegativeBinomials, with a and and a shared dispersion. The actual foreground/background signal within a modality was modelled as a that depends on the latent space, and which is multiplied with the empirical library size to get. For Visium Highly Multiplexed Protein, was modelled as a latent variable specific for each gene. for Visium Highly Multiplexed Protein and for RNA-seq were modelled as deterministic functions depending on the use case as described in the following paragraphs. The posterior of the probabilistic graphical model was inferred using black-box variational inference, in which the variational distribution was specified as a diagonal Normal distribution, transformed into the correct domain using transforms. Free parameters within this model were optimized using gradient descent, with the ELBO as loss function and Adam as optimizer as implemented in Pytorch (pytorch.org). We used a learning rate of 0.01 for variational parameters, and 0.001 for parameters of the amortization functions.

64To calculate the average expression of each gene within a cell type, we used a linear model in which both and were modelled as a latent variable specific for each gene and cell type. The for nuclei were multiplied with a gene-specific correction factor (optimized as a latent variable) that corrected for differences between scRNA-seq and snRNA-seq. Given that spatial transcriptomics data sequences the whole cell, the uncorrected values were used for spatial deconvolution.

65To infer the proportions of each cell type within a spot, we used a model in which the gene expression is modelled as a linear combination of cell type proportions and average expression in each cell type:

66For we adapted the values from the reference, but included- A capture bias per gene, which corrects for technical and biological differences between spatial and sc/sn-RNA-seq. The capture bias was modelled as a latent variable with prior - A red blood cell cell type, which was not included in the reference dataset but nonetheless had a dominant presence in the spatial data. The of this cell type was set to zero for all genes except Hbb-bt, Hbb-bs, Hba-a1, Hba-a2 for mouse and HBB, HBA1, HBA2 for human, which were modelled as free parameters.- Similarly, the expression of complement factors (C3, C2, C4B/C4b) within hepatocytes was modelled as free parameters.

67- A capture bias per gene, which corrects for technical and biological differences between spatial and sc/sn-RNA-seq. The capture bias was modelled as a latent variable with prior

68- A red blood cell cell type, which was not included in the reference dataset but nonetheless had a dominant presence in the spatial data. The of this cell type was set to zero for all genes except Hbb-bt, Hbb-bs, Hba-a1, Hba-a2 for mouse and HBB, HBA1, HBA2 for human, which were modelled as free parameters.

69- Similarly, the expression of complement factors (C3, C2, C4B/C4b) within hepatocytes was modelled as free parameters.

70A background signal shared for all spots was also modelled as follows:

71With a latent variable specific to each spot and a latent variable specific to each gene.

72A likelihood ratio test was used to assess whether a cell type was significantly present in a spot. Specifically, if is the gene expression of all genes at a particular spot, we used Monte Carlo samples from the posterior to estimate:

73A cell type was deemed significantly present if the log-likelihood was higher than 10.

74The zonation of spots was modelled as a univariate latent variable specific to each spot. This latent variable influenced the gene expression using a spline function by using a gaussian basis function with 10 knots at uniform fixed positions. The coefficients of this spline were modelled as a latent variable specific for each gene, with prior a Gaussian random walk distribution, and the step. was determined empirically as 2 times the standard deviation of the log1p transformed expression values in the whole dataset. The variational parameters of the zonation and were not optimized directly but were estimated using an amortization function. This amortization function used the count matrix as input, and estimated the variational parameters using the following layers: Linear (with 100 output dimensions), BatchNorm, ReLU, Linear (again with 100 output dimensions), ReLU, and a final Linear layer. This amortization function was used to transfer the zonation onto a different dataset, i.e., 1) to transfer the zonation trained on mouse spatial transcriptomics onto mouse Visium highly multiplexed protein and 2) to transfer the zonation trained on human low steatosis (<10%) onto human high steatosis (>30%).

75To determine the differential abundance of a cell type across zonation, the significant presence of a cell type within a spot was modelled using a spline function with the zonation of a cell type as input. The coefficients of this spline function were modelled as a latent variable with the step size. To determine differences in abundance between patients with high and low steatosis, we first modelled the zonation on human data on patients with steatosis < 10%. Potential interaction effects between zonation and steatosis status were then modelled using a spline function as before, but with a separate set of coefficients for both high and low steatosis. A likelihood ratio test was then used to determine whether this interaction was present significantly, by comparing the likelihood of this model with a model with shared coefficients.

76To analyze cell-cell communication in the hepatic mac niches, we applied Differential NicheNet, which is an extension of the default NicheNet pipeline to compare cell-cell interactions between different niches and better predict niche-specific ligand-receptor (L-R) pairs. It uses a flexible prioritization scheme that allows ranking L-R pairs according to several properties, such as niche- and region-specific expression of the L-R pair, ligand activity, and level of database curation. This in contrast to the default NicheNet pipeline which prioritizes expressed L-R pairs solely based on ligand activity predictions. All analyses were conducted according to the Differential NicheNet tutorial (https://github.com/saeyslab/nichenetr/blob/master/vignettes/differential_nichenet.md). As input to the Differential NicheNet pipeline, we used the data after normalization via SCTransform and integration of scRNA-seq and snRNA-seq according to the Seurat procedure for integration.

77For the mouse analyses, Differential NicheNet was first performed for the following 3 niche comparisons: 1) KCs versus central vein macs; 2) KCs versus capsule macs; 3) KCs versus LAMs. Following sender cell types were considered for these niches: KC niche: periportal hepatocytes, periportal LSECs, and periportal stellate cells; Central vein mac niche: central vein ECs and central vein fibroblasts; Capsule mac niche: mesothelial cells and capsule fibroblasts; LAM niche: cholangiocytes and bile duct fibroblasts.

78Because of the preferentially periportal localization of KCs in the mouse liver, we also included a ‘region specificity' factor in the Differential NicheNet prioritization framework. This was done to increase the ranking of ligands that are more strongly expressed in periportal than pericentral niche cells. Periportal sender cells were determined after subclustering based on the following markers: Hal and Sds for hepatocytes; Mecom, Msr1, and Efnb2 for LSECs; Ngfr, Igfbp3, and Dach1 for stellate cells.

79In the heatmap (Figure S8G), we show the prioritization scores of the top 40 ligands (and their highest scoring receptor) in the KC niche (score averaged over the 3 analyses), and of all the non-KC niche L-R pairs with a prioritization score ≥ the score of the lowest scoring KC L-R pair of this top 40. For each L-R pair/niche combination, we only displayed the score of the sender cell with the highest score (e.g. for the Csf1-Csf1r interaction in the KC niche, the score is shown for the LSEC-KC interaction because that score was higher than for Stellate–KC and Hepatocyte–KC; in the LAM niche, the score of Csf1-Csf1r is shown for the bile duct fibroblast – LAM interaction and not for the cholangiocyte–LAM interaction, etc.).

80Because of the strong concordance between the top-ranked L-R pairs in these 3 non-KC mac niches, it was decided to also conduct a subsequent analysis in which the KC niche is compared against all non-KC hepatic mac niches combined. For this final ‘KC versus all non-KC mac analysis’, KCs were compared to central vein macs, capsule macs, and LAMs together, with the same sender cell types as described here above (but now analyzed together).

81For the human analyses, Differential NicheNet was performed to compare the KC niche with the non-KC mac niches (similarly as the final analysis in mouse). For the KC niche, all hepatocytes, LSECs, and stellate cells were selected as sender cells; and KCs as receiver cells. For the non-KC mac niche, cholangiocytes, fibroblasts, and central vein ECs were considered as the sender cells; Mat. LAMs, Imm. LAMs, and Mac1s as the receiver cells (Figure 4H).

82To find KC-niche-specific L-R pairs that are conserved across mouse and human, the individual mouse and human prioritization scores were averaged to form a ‘conservation score’. The 40 ligands (and maximally 3 of their highest scoring receptors) with the highest conservation score were selected for further analysis (note: the L-R pair should be expressed by the same sender-receiver pair in both species). In the circos plot (Figure 6C; ), only a subset of these top L-R pairs is shown to keep the figure clearly interpretable. Following ligands were not shown: ITGA9, SEMA6D, JAM3, ITGB1 (stellate cells); ITGA9, F8, CD274, HSP90B1 (LSECs); C5, F9, F2, FGA, TF, TTR, COL18A1, COL5A3, SERPINA1, SERPINC1 (hepatocytes). The depicted target genes are KC-specific in both mouse and human, and a top-predicted target according to the NicheNet ligand-target regulatory potential scores. NR1H3 was manually added as a NOTCH2 target based on recent studies.

83Resected human liver was fixed in 4% formalin for 24-48h and subsequently embedded in paraffin. Samples were stored for 10-15 years at RT before analysis. Sections of 6 μm thick were cut using a Microm HM360 and mounted on a polarized glass slide. These sections were deparaffinized in xylene and rehydrated in a graded ethanol series. Antigen retrieval was performed by immersing the samples for 5 min in pH 8.3 TRIS-EDTA at 98°C. Slides were then cooled to RT and washed in PBS. Confocal staining was performed as described above.

84BM was isolated from the tibia and femur of mice by centrifugation. Red blood cells were lysed and single cell suspensions were stained with antibodies for flow cytometry. BM monocytes were sorted as live CD45+ CD11b+ Ly6G- Ly6C+ CD115+ cells using a BD FACSAria III. Monocytes were resuspended in DMEM/F12 media supplemented with 10% FCS, 30ng/ml CSF1, 2mM Glutamine and 100U/ml penicillin and streptomycin. 150,000 monocytes were seeded in each well of an adherent 24-well plate pre-coated with bovine collagen type I and cultured overnight (37C, 5% CO2). The following day 0, 25 or 50ng/ml of ac-LDL was added. Ac-LDL was kindly provided by Sophie Janssens, Ghent, Belgium who received the material from Wilfried Le Goff, Paris, France. 14 hours later cells were harvested and live F4/80+ cells were FACS-purified in RLT plus buffer containing 1% β-mercaptoethanol. RNA isolation, cDNA synthesis and qPCR were performed as described above.

85Bone marrow chimeras were generated as described previously. Briefly, 6-12 week old Clec4f-Dtr mice (CD45.1) were anaesthetized by intraperitoneal administration of Ketamine (150 mg/kg) and Xylazine (10 mg/kg). Mice were lethally irradiated with 8 Gy, with the livers being protected with a 3-cm-thick lead cover. Once recovered from the anesthesia, mice were reconstituted by intravenous administration of 5-10×106 BM cells from CD45.2 Acvrl1fl/fl or Fcgr1-CrexAcvrl1fl/fl mice. 4 weeks after reconstitution mice were administrated a single dose of 500ng DT via intraperitoneal injection to deplete KCs. Chimerism was assessed 7 or 13 days later by flow cytometry and compared with chimerism levels in blood Ly6Chi monocytes.

86Clec4f-Dtr mice were administered 10mg/kg ALK1Fc, TGFβRIIFc or appropriate isotype controls (hIgG1 and mIgG2a; Acceleron Pharma) by intraperitoneal injection on days -1, 2, 3 and 5. On Day 0 mice were also administered a single dose of 500ng DT i.p. to deplete KCs. Livers were harvested at day 7 to assess KC development.

87In all experiments, data are presented as mean ±SEM and/or individual data points are presented unless stated otherwise. Statistical tests were selected based on appropriate assumptions with respect to data distribution and variance characteristics. Details of the precise test used for each analysis can be found in the figure legends. Statistical significance was defined as p<0.05. Sample sizes were chosen according to standard guidelines. Number of animals/patients is indicated as ‘‘n’’. The investigators were not blinded to the group allocation, unless otherwise stated.

88The sc/snRNA-sequencing, CITE-seq FCS files and spatial transcriptomics datasets will be made available for visualization, analysis and download at www.livercellatlas.org.