A rat liver cell atlas reveals intrahepatic myeloid heterogeneity

Summary The large size and vascular accessibility of the laboratory rat (Rattus norvegicus) make it an ideal hepatic animal model for diseases that require surgical manipulation. Often, the disease susceptibility and outcomes of inflammatory pathologies vary significantly between strains. This study uses single-cell transcriptomics to better understand the complex cellular network of the rat liver, as well as to unravel the cellular and molecular sources of inter-strain hepatic variation. We generated single-cell and single-nucleus transcriptomic maps of the livers of healthy Dark Agouti and Lewis rat strains and developed a factor analysis-based bioinformatics analysis pipeline to study data covariates, such as strain and batch. Using this approach, we discovered transcriptomic variation within the hepatocyte and myeloid populations that underlie distinct cell states between rat strains. This finding will help provide a reference for future investigations on strain-dependent outcomes of surgical experiment models.


INTRODUCTION
The liver is a multitasking organ that performs a remarkably diverse set of functions including nutrient metabolism, regulation of immune responses, and protein synthesis.Despite its highly regenerative and tolerogenic nature, 1,2 inflammatory end-stage liver diseases such as druginduced liver injury, hepatitis infection, hepatocellular carcinoma, and autoimmune hepatitis are common. 3Despite the recent advancements in medical strategies to treat acute liver disease, [4][5][6] the development of therapeutic options is limited by our incomplete understanding of the cellular landscape of the liver in non-mouse animal models.The liver is composed of multiple cell types with complementary functions, including hepatocytes, biliary epithelial cells (cholangiocytes), mesenchymal cells (stellate cells and vascular smooth muscle cells [VSMCs]), myeloid cells, liver sinusoidal endothelial cells (LSECs) and multiple other immune cell populations. 7Hepatocytes make up the majority of liver volume and are involved in metabolism and drug detoxification, among other functions that are often zonated along the hepatic lobule. 8,9Myeloid cells are distributed throughout the liver and can adopt pro-inflammatory or anti-inflammatory roles, with phenotypic characteristics of recently recruited monocytic myeloid cells and more tissue-resident Kupffer cell-like populations, respectively. 7Current animal models used to recapitulate and study liver pathology include the porcine, murine, and rat models.A key advantage of the rat model (Rattus norvegicus) is its large size, which allows for better vascular access for disease models that include surgical interventions such as hepatectomies, 10,11 hepatic ischemia reperfusion-induced injury models, transplant injury, 12 and fibrocirrhotic bile duct ligation models. 13,14o date, our understanding of the rat liver has been informed by technologies such as bulk RNA sequencing (RNA-seq), [15][16][17][18] transcriptome microarrays, [19][20][21][22] immunohistology, 23,24 targeted qPCR, 19,23,25 and tandem mass spectrometry. 17These approaches have uncovered the presence of major expected hepatic populations in the rat liver; 19 however, the relatively low resolution and targeted nature of these approaches do not allow us to have a holistic understanding of how the interaction between diverse hepatic cells shapes the liver environment.Single-cell RNA sequencing (scRNA-seq) technology is a powerful tool for the unbiased profiling of heterogeneous tissues.While both human 7,[26][27][28][29][30] and murine [31][32][33][34] livers have been well studied at the single-cell level, the rat liver has remained poorly annotated.Studies using the rat model 13,35,36 demonstrate strain-associated differences in the liver and inflammatory disease severity.For example, while both Dark Agouti (DA) and Lewis (LEW) strains are prone to Th1-skewed responses in the joints of treatment-induced autoimmunity models, 37 DA rats appear to paradoxically have similar innate responses as autoimmune-resistant Albino Oxford (AO) rat strains. 38,39In the context of orthotopic transplantation, LEW recipient liver macrophages are better able to stimulate T cell proliferation in comparison to DA. 35 This provides a rationale for single-cell examination of rat strain-specific differences at baseline.We developed a single-cell transcriptomic atlas of the healthy liver, based on DA and LEW rats, to provide a reference atlas of the healthy rat liver.Our combined use of single-cell (sc) and single-nucleus (sn) RNA sequencing (RNA-seq) as well as spatial transcriptomics enabled us to discover cellular and molecular sources that drive the inter-strain variation and will be helpful for understanding strain-dependent hepatic disease models and rat liver biology in future studies.

The cellular landscape of the healthy rat liver
We generated a multi-strain single-cell transcriptomic map of the healthy rat liver to help examine the cellular complexity in this model system.Single-cell transcriptomes were generated from total liver homogenates (TLHs) of four 16-18 week-old healthy male rats following 2-step collagenase digestion (Figure 1A).Two livers from each of the DA and LEW strains were sampled, and a standard scRNA-seq mapping pipeline was applied (Figure 1B).In total, 226,270 single cells were called by the 10x Genomics Cell Ranger software and 23,036 passed additional quality control filters and were included in the final map (see STAR Methods, Figures 1C, 1D, and S1; Tables S1 and S2).Significant batch effects were evident while integrating the four rat samples; therefore, the Harmony 40 integration method was used to reduce the inter-sample technical confounding effects.After applying this batch correction, all clusters were represented by all animals, demonstrating that integration worked well (Figures 1E and 1F).Liver tissue from an additional two pairs of LEW and DA rats were processed for 10x Genomics snRNAseq to further inform parenchymal cell identities (Figure 1A).These samples went through standard quality control steps (see STAR Methods; Figure S2; Table S1) and were batch-corrected using the Harmony integration tool.The resulting map contained 12,497 nuclei.Cell populations were annotated, based on known marker genes, using top differentially expressed (DE) genes 41 (see STAR Methods; Figures 1G and 1H; Tables S3 and S4).To resolve the spatial distribution of rat hepatic cell populations, we conducted spatial transcriptomics on two healthy Wistar rat liver samples using 10x Genomics Visium technology.These samples were then quality controlled in a similar manner to scRNAseq data (Figures S3; Table S1).

Hepatocytes
Hepatocytes, organized in functional units referred to as lobules, make up the majority of the liver volume (Figure 1I).Many of their critical biological functions are zonated based on their spatial organization from the center of the lobule near the central vein (CV) to the outer regions near the portal triad.Data from both sc-and snRNA-seq protocols identify hepatocyte-like clusters, based on their correlation with hepatocytes of the mouse liver atlas (Figure S4), and expression of hallmark hepatocyte markers without high expression of immune endothelial and mesenchymal genes (Figures 1G and 2A-2E).
Comparative gene expression analysis of our data to a bulk RNA-seq dataset of laser capture microdissected zonated regions of the healthy mouse liver lobule 42 revealed poor zonated marker distribution in the scRNA-seq dataset compared to snRNA-seq (Figure 2C), as has been observed before. 27,43To further resolve hepatocyte cluster identity, the datasets were compared with a spatial transcriptomics map of the rat liver from two Wistar rats.Principal-component analysis (PCA) 44 of these samples revealed the largest dimension of variation was related to lobule zonation (Table S5).PC1 and PC2 in both samples indicate clear histological periportal to central venous zonation patterns (Figures 2D and S5).Additionally, key periportal human (Tf, Hmgscs1 27,43 ) and mouse (Ass1, Arg1 42 ) genes as well as periportal markers from rat proteomic studies (Gls2, Srd5a1, Orm1 45 ) are positively enriched in PC1 and PC2.Known pericentral markers (Oat, Sult1e1, Cyp2e1, Glul 27,42,43 ) are negatively enriched, reinforcing that these principal components represent zonation patterns 45 (Figures S6 and S7; Data Portal).Pathway enrichment analysis of the PCs was performed to further validate that PC1 and PC2 represent zonation features.Periportalbiased processes such as immunity, angiogenesis, lipid beta-oxidation, fatty acid catabolism, and gluconeogenesis regulation 8,46 are found in the positive side of PC1 while the negative side of PC2 is enriched in pericentrally biased metabolic processes, such as lipogenesis and various steroid and xenobiotic metabolic processes 9 (Figure 2E).Examination of key markers and correlation analysis between PC1 and PC2 of the spatial transcriptomics data and snRNA-seq hepatocyte clusters shows a clear presence of periportal and pericentral hepatocyte populations (Figures 2B-2D).These findings suggest that pericentral and periportal programming is well preserved across species.

Mesenchymal cells
The hepatic mesenchymal fraction includes populations such as hepatic stellate cells (HSCs), VSMCs, and fibroblasts (FBs). 27Mesenchymal cells anatomically reside between sinusoidal endothelial cells and hepatocytes and are involved in vitamin A storage, extracellular matrices (ECMs) synthesis, maintenance of hepatocyte function, 47 and regulation of sinusoidal circulation. 48These populations also help regulate immune responses during inflammation, 1 but upon activation can also be a source of maladaptive extracellular matrix deposition, as in the case with liver fibrosis. 49e annotated two clusters in our scRNA-seq map (scClusters 7 and 14) and one cluster (snCluster 24) in our snRNA-seq map, as mesenchymal-like based on DE genes including extracellular matrix proteins (Ecm1) and type III collagen alpha 1 (Col3a1) (Figures 1G, 1H, 2A, and  3A) which are essential to the role of HSCs in extracellular matrix deposition and have previously been described as mesenchymal genes 27,43 (expanded markers shown in Figure S8).

Endothelial cells
The hepatic endothelium consists of LSECs and vascular endothelium (portal and central venous endothelium).LSECs are a specialized endothelial population that line the hepatic sinusoids and contribute to the regulation of hepatic blood pressure, nutrient uptake, and the maintenance of HSC quiescence. 51,52Immunohistochemical staining in mice has described general endothelial cells in the liver as expressing high levels of Cd31 (Pecam) and Cd103 (Eng), periportal LSECs as expressing high levels of Cd36, with low levels of Lyve1, and central venous LSECs as expressing high levels of Cd32b and Lyve1. 53,54However, in rats, endothelial zonation has yet to be reported.
We identified two populations of Ptprc À cells in the scRNA-seq map (scClusters 3 and 11) and two populations in the snRNA-seq map (snCluster 11 and 30).These populations were annotated as endothelial-enriched based on the expression of Calcrl and Ramp2, which is .SnRNA-seq and spatial transcriptomic profiling of the rat liver resolves hepatocyte zonation Four additional rat liver samples were added and sequenced using snRNA-seq to better characterize hepatocyte and cholangiocyte populations and verify the strain variations identified based on the scRNA-seq TLH map.(A) UMAP projection of four snRNA-seq samples where cells are colored based on cell-type annotation.(B) Heatmap representing the average gene expression of zonated genes based on spatial data within the snRNA-seq clusters.(C) Pearson correlation between the average gene expression of the genes across snRNA-seq hepatocyte clusters and the nine layers of mouse liver cells was calculated (see STAR Methods).Mouse liver layer-9 is more periportal and layer 1 is pericentral.Red represents a positive correlation, and blue represents a negative correlation.(*: p value <0.05, **: p value <0.01, ***: p value <0.001).(D) Projection of zonation signature scores, captured by PC1, across the spatial transcriptomics spots of two healthy Wistar rat liver cryosections.The top negatively loaded genes in PC2 (and PC1) of both samples are enriched in pericentral markers, and the top positively loaded genes in PC1 (and PC2) factors are enriched in periportal markers.Red and blue represent high periportal and pericentral zonation scores, respectively.The two heatmaps represent the Pearson correlation between the zonation factors PC1 and PC2 and the average expression of snRNA-seq hepatocyte clusters.Both PC1 and PC2 are positively correlated with PP-like clusters Hep0 and Hep1 and negatively correlated with CV-like clusters Hep2 and Hep3.The asterisk and triangle symbols indicate the factors used for pathway enrichment analysis.(E) Pathway enrichment analysis using GSEA (gene set enrichment analysis) to examine active cellular pathways in periportal and central venous regions of the healthy rat liver based on spatial PC1 and PC2 loadings visualized as an enrichment map.The pathways enriched in the pericentral and periportal areas are based on PC2 (asterisk) of liver cryosections-A (left) and PC1 (triangle) of liver cryosections-B (right) respectively.Each circle represents a gene ontology (GO) biological process term.The size of the circles represents the number of genes in that pathway and blue lines indicate significant gene overlap.

ll OPEN ACCESS
involved in adrenomedullin signaling pathways 55 (Figures 1H, 2A, 3A, and S10A).ScCluster 3, the most abundant endothelial cell population, was characterized by enriched expression of Lyve1, Fcgr2b, Sparc, and Stab2 (Figures 1H and S10A) with little expression of Vwf (Figure S10A) suggesting an LSEC identity.By correlation analysis, both scClusters 3 and 11 were similar to mouse sinusoidal, inflammatory, and cycling endothelial populations (Figure S4).These clusters did not show differential expression of known zonated endothelial genes such as Rspo3 56 and Clec4g, and both clusters expressed high levels of Fcgr2b (known to be enriched in CV LSECs 54 ) and Aqp1 (known to be enriched in periportal LSECs 57 ) (Data Portal).Endothelial genes such as Lyve1 and Vwf found in scCluster 11 were also found snClusters 11 and 30 (Figure 2A; Data Portal).SnRNA-seq endothelial populations were subclustered for increased resolution, and comparisons were made to our spatial transcriptomic data (Figure S10B).The resulting subcluster 3 had a stronger expression of PC1-enriched periportal markers (Vwf and Ltbp4) with little expression of Lyve1, 58 while subcluster 1 and 0 expressed higher pericentral-associated genes such as Lyve1, Fcgr2b, and Bmp2 (Figures S10C and S10D; Data Portal).Further examination of other known periportal markers in our spatial transcriptomics data did not reveal clear endothelial zonation patterns (Figure S10D; Data Portal), perhaps due to the low capture of endothelial genes by the Visium spatial transcriptomics platform. 43

Biliary epithelial cells
Cholangiocytes are liver-specific biliary epithelial cells whose primary function is the production and modification of bile as it flows along the biliary tract. 59In line with previous literature, cholangiocytes were poorly captured with scRNA-seq and were only detected by our snRNA-seq map. 43SnCluster 29 of the snRNA-seq map was identified as being enriched in the expression of Epcam, Sox9, Anpep, and Anxa4, 43 resulting in a total of 108 cholangiocyte-like cells with Anxa4 and Epcam showing a periportal distribution on the spatial transcriptomics map (Figures 3A and 3C).

Myeloid cells
The liver contains more resident myeloid cells than any other solid organ in the body. 60Tissue-resident myeloid cells exhibit immense phenotypic plasticity and can perform a diverse set of functions.Depending on the local immune microenvironment and external stimuli, bone marrow-derived monocytes can be recruited to the liver, where they participate in both liver injury and tissue repair.In comparison, the primary function of sessile resident myeloid cells is to clear debris, in addition to mediating the tolerogenic environment of the liver in the steady state. 61,62ur single-cell analysis revealed multiple clusters of Cd68 + myeloid-enriched cells.Cd68 + myeloid scClusters 5 and 10 were characterized by enriched expression of Marco, Vsig4, Cd5l, Cd163, and Hmox1 (Figure 1H, see extended gene expression in Figure S11).These clusters appear to be more Kupffer cell-like due to the expression of key genes (Marco, Cd5l, Clec4f) which have been previously described to annotate more tissue-resident myeloid populations. 63Specifically, Vsig4 is a co-inhibitory ligand that has a hepatoprotective role in maintaining the intrahepatic tolerance required to suppress triggered immune responses 64,65 and has been shown to be highly expressed in murine Kupffer cells (KCs), 64,65 as well as being a core KC gene in pig and macaque KCs. 34These findings may suggest a tolerogenic role of Marco + Cd5l + Cd68 + cells, which are represented by snCluster 19 of our snRNA-seq map (Figure 3A).Our analysis of scCluster 9 revealed a mixed cluster of Ptprc+ immune cells enriched for Cd68 + myeloid cells (enrichment and subclustering of immune cells are discussed in the following).ScCluster 9 is enriched in described macrophage and monocyte markers (Cd68, Cd74, Lyz2, and major histocompatibility complex (MHC) class 1-related genes) without the expression of Vsig4 and Marco, suggesting that it is enriched for recently recruited macrophage/monocyte populations.However, ScCluster 9 also contains additional immune populations such as T cells (Cd3e), conventional denditic cells (cDCs) cDC1s (Clec9a, Xcr1, Batf3, Irf8) (Figure S11), cDC2s (Clec10a, Irf4, Sirpa), and plasmacytoid dendritic cells [pDCs] (Siglech).This macrophage/monocyte cluster was represented by snCluster 33 of our snRNA-seq dataset, but due to lower capture of non-myeloid immune cells by snRNA-seq technologies, 43 it contains only a minor Ighm + B cell population (Figure 3A).To resolve zonation, an examination of key myeloid markers (Cd68, Cd163) and key genes of KC-like myeloid cluster genes (Cd5l, Marco, Aif1, Hmox1, Clec4f) in PC1 was performed.The positive enrichment in PC1 suggests the presence of myeloid cells is skewed toward the periportal areas (Figures 3D-3F; Data Portal).To validate this enrichment, quantification of immunohistochemistry stainings of Cd163, Cd68, and Hmox1 was performed using a publicly available (G) Representative spatial distributions of CD68 + cells in the rat liver lobule.Rectangular layers 350 um wide were drawn from the portal tract (layer 1) to the central vein (layer 10) region.Digital images were scanned at 203 magnification.The scale bar represents 100 mm in the full image and 20 mm in the enhanced area.Each rectangular layer is referred to as a region of interest (ROI).(H) Quantification of CD68 + cell densities (#CD68+ cells/layer mm2) in the liver lobule for DA and LEW rats.30 ROIs were assessed per strain across three animals.A higher number of CD68 + cells were detected near the periportal area.No significant strain-specific differences in the spatial distribution of CD68 + cells were noted.(I) Representative spatial distributions of CD163+ cells in the rat liver lobule.(J) Quantification of CD163+ cell densities (#CD163+ cells/layer mm2) in the liver lobule for DA and LEW rats.The statistical analysis reflects similar strain and zonation patterns as CD68.Statistical significance was determined using a two-way ANOVA followed by Sidak's multiple comparisons test.QuPath-based image analysis protocol. 66,67This analysis confirms the periportal-biased nature of non-inflammatory myeloid cells (Figures 3G-3J and S12).

Varimax PCA analysis uncovers biological sources of variation between rat strains
To better understand strain-specific differences in our map, we applied varimax PCA, 44,68,69 a matrix factorization method, to separate DA and LEW signals (principal components, or factors) in the data from other signals for further interpretation (Figure 1B, Figure 4, Table S5).To identify factors that can explain strain-specific differences, we used a random forest to predict strain labels from the factors identified per cell and discovered the factors most important for the strain label classification (Figure 4A).We also identified principal components that explain celltype signals using correlation analysis (Figures 4B and S13).The resulting factors were interpreted using pathway and gene set enrichment analysis (see STAR Methods).Using this approach, two main strain-specific factors (varimax PC5 and 15) were identified (Figures 4A and S14) within the scRNA-seq TLH map.The strongest strain-specific signal is observed with varimax PC5, which affects all cells in the data (Figures 4C  and 4D).Genes with the strongest association with this factor are hepatocyte markers (Apoc1, Fabp1, and Cytochrome p450 genes), suggesting that this factor mainly represents strain variations within the hepatocyte populations (Figure 4E).The global association of this factor with all cells in the scRNA-seq dataset is likely a cell-dissociation procedure artifact caused by fragile hepatocytes leaking RNA into the cell homogenate before sequencing (Figure S15). 43DA strain-associated genes in this factor are enriched in nuclear receptors, such as Hnf4a, Pparg, and Esr1 (Table S6) (Figure S16).Pparg promotes de novo lipogenesis and fat accumulation in hepatocytes. 70,71This hepatocyte-specific strain signal was confirmed in the snRNA-seq dataset (Figure S17; Table S5).The second-strongest strain-specific signal is varimax PC15, which is mainly associated with myeloid populations of both rat strains (Figures 4F and 4G), as confirmed by the genes with the strongest association with this factor (Figure 4H), the expression pattern of Marco, Visg4, Cd68, and Lyz2 marker genes (Figure 4I), and the correlation with myeloid cells in our map (Figure 4J).
Comparing the expression level of the top varimax PC15 genes in the myeloid cells of the two strains confirms the strain-specificity of this factor (Figures 5A-5C).Pathway analysis identified higher activation of lymphocyte-mediated immune responses, lymphocyte migration and chemotaxis, response to interferon, and allograft rejection pathways in LEW compared to DA Marco-enriched myeloid cells (Cluster 5) (Figure 5D; Table S7).This factor is enriched in myeloid and T cell differentiation transcription factors (TFs) (Figures 5E and 5F).LEW-enriched TFs include Irf8, Irf1, Spi1, Pou5f1, Stat4, and Stat5a, which are mostly inflammatory process-associated genes present in chronic diseases like rheumatoid arthritis [72][73][74] (Figure 5F).3][74] The DA-specific TFs include PPAR-g, Nucks1, Runx1, Mitf, and Gata1, which have been described more broadly in the literature 73,[75][76][77] (Figures 5E and S18; Table S6).][80][81] No strong myeloid-specific strain-related varimax factors were discovered using the snRNA-seq map, which can be explained by the lower representation of non-inflammatory myeloid cells within the snRNA-seq map (276 cells) compared to the scRNA-seq TLH map (1,668 cells). .Varimax PCs capture rat hepatic cell identity signatures and strain-specific differences (A) Bar plot representing the feature importance scores (mean decrease Gini impurity) of the top 20 features (varimax factors) of the random forest model trained to predict the strain attributes of the rat hepatic cells.Varimax PC5 and 15 are the most informative features to differentiate cells of each strain from another, which indicates the two factors have captured strain-related variations within the map.(B) A correlation heatmap between the average gene expression of each cluster and the loading scores of varimax factors (capturing the contribution of all genes to a factor).Columns are varimax factors and rows are cell populations.Each cell-type cluster is defined by key marker genes, and dark red or blue indicates that the expression of a marker gene set is positively or negatively correlated, respectively, with a particular varimax factor.A high absolute correlation value indicates a match between a varimax factor and a cell-type cluster.(C) The projection of cells over varimax-1 and 5 indicates that the cells from each strain form distinct clusters over varimax-5.(D) Boxplot indicating the distribution of varimax-5 score over each strain.Cells from DA and LEW strains represent significantly different varimax-5 scores (Wilcoxon-test p value <2.2e-16), indicating that varimax-5 has captured strain differences.(E) The top 10 genes on the top (left table) and bottom (right table) of the varimax-5 loading list mainly contain known hepatocyte markers, indicating that varimax-5 has captured hepatocyte-specific strain differences.Genes with high positive scores (left table) are associated with the DA strain and genes indicating negative loading scores (right table) are LEW-related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(F) Projection of cells over varimax-1 and 15 indicates that a population of cells from each strain (dotted lines) forms distinct clusters over varimax-15.Annotation of the selected cells indicates that they are mainly from the Marco+ myeloid cluster 5. (G) Boxplot indicating the distribution of hepatic cells based on strain over varimax-15  2J).However, we were able to validate the myeloid-specific strain factor identified in the scRNA-seq TLH map.We selected the top 10 PC15-associated genes and calculated their enrichment within the myeloid cells of the snRNA-seq map.In line with our scRNA-seq results, varimax PC15 signatures show strain differences within the myeloid cells of the snRNA-seq samples (Figures S19A-S19F).We also evaluated the DE genes between DA and LEW within the myeloid population of the snRNA-seq samples using a generalized linear mixed model considering covariates, like sample and strain. 82The most DE genes include the top strain-related genes identified by the varimax analysis approach (Itgal, RT1-A1, Timp2, Lilrb3l, RT1-T24-3) along with some expected ambient-RNA transcripts (Figure S19G).These results suggest that the baseline hepatic microenvironment in the LEW rat is more pro-inflammatory compared to the DA strain and highlight myeloid cells as potential drivers of the enriched inflammatory pathway activation in LEW rats.We then considered whether myeloid cell frequency in the DA and LEW livers may be influencing the inflammatory status of LEW rats.Alterations in cell-type frequencies in scRNA-seq data are confounded by sample-specific dissociation efficiency.Therefore, we employed immunohistochemistry to compare the frequency of CD68 + cells between the two strains.However, quantification of CD68 staining showed no significant difference in the frequency of CD68 + cells in LEW vs. DA (Figures 3G and  3H).These results suggest that the variations in inflammatory potential are not likely caused by differences in the frequency of intrahepatic CD68 + cells.
Immune enrichment maps rat lymphocyte and myeloid populations at higher resolution Our scRNA-seq TLH map and snRNA-seq map contained hepatocyte-derived ambient RNA, as expected 43 (Figure S20), which interfered with immune cell marker identification and resulting immune cell annotation.To provide a more detailed resource of rat hepatic immune cells, two additional immune-enriched samples were mapped (Figure 6A).These samples underwent additional washing steps and red blood cell depletion to reduce the hepatocyte-released ambient RNA (Figure S20).The percentage of cells annotated as hepatocytes decreased from 71.14% in the scRNA-seq TLH map to 49.11% in the immune-enriched map.The general immune cell marker, Ptprc, was expressed in 24% of the total cells in the immune-enriched map compared to 4% within the initial map (Figures 6B and 6C).Unfortunately, the scRNA-seq TLH and immune-enriched maps could not be integrated computationally, presumably due to the technical differences in their generation (Figures S21A-S21C).The varimax-based pipeline was also ineffective to deconvolute the sources of variation in the merged dataset of both sets of samples (Figure S21D).Consequently, the immune-enriched samples were analyzed separately.In total, 3,830 (1,161 + 2,669) single cells from the DA and LEW samples were integrated into the immune-enriched map after quality control (see STAR Methods; Figure S22).Similar to the TLH map, the immune-enriched samples were batch-corrected and the final clusters represent cells from both DA and LEW rats (Figures 6D and 6E).The clusters were annotated based on the same approaches used for the initial samples (Extended results; Table S8).The immune-enriched map has captured a more diverse set of liver-resident immune cells (Figures 6A and 6F), enabling a more detailed description of these cell populations (Figure 6G) compared to the scRNA-seq TLH and snRNA-seq maps.A comparison of the scRNA-seq TLH and immune-enriched maps using correlation analysis confirmed that the immune-enriched map provides a higher resolution of lymphocytes and myeloid cells (Figure 6H).As a refinement to the immune annotations in the TLH map, individual populations of Cd3 + T cells (clusters 10), natural killer (NK)-like cells (cluster 7), B cells (cluster 12), and pDCs (cluster 17) were identified (described in the following) in the immune-enriched map (Figures 6G and 6H).Cluster 10 was characterized by enriched expression of Cd3 + T cell markers (Cd3g, Cd3e, Cd3d, Coro1a) (Figure S23).Cluster 12 identified a subset of cells enriched for B cell genes Cd19, Ms4a1 (Cd20), Ighm, Cd74, Cd79b, and Fcmr, with no expression of Ighd or Ighg, suggesting that this cluster might be Cd19 + Cd20 + IgM + IgD -immature B cells 83 (Figure S24).The correlation heatmap (Figure 6I) indicated high gene expression similarity with the mouse 84 B cell populations, supporting that this is a B cell population.Enriched gene expression in cluster 17 correlates with both monocyte-like macrophages (Cd74 and Tyrobp), and pDCs (Siglech, 34 Ptprcap, 85 and Ptcra 85 ) (Figure S25).When comparing the expression of this cluster to the mouse liver cell atlas, we see a high correlation with pDCs (Figure 6I), suggesting that the predominant cellular population of this cluster may be pDC enriched. 84Cluster 11 displays a correlation with monocytic macrophages and dendritic cells and similar DE genes as scCluster 9 suggesting it is a mixture of recently recruited immune  populations (Figure S26).DE genes in cluster 7 include Tbx21 [aka T-bet], Ncr1, Prf1, Nkg7, Ccl5, Cd8a, Gzmk, Klrd1, and Cd7, with low expression of Cd3d, suggesting it is an NK-like population 7,[26][27][28][29][30]86 (Figure S27). Theexpression of top genes in this cluster correlated with the NK cell population in the mouse dataset (Figure 6I) reinforcing that this cluster is an NK-enriched cluster.The Ptprc + clusters of the immune-enriched map were subclustered for further evaluation (Figures S28-S33; Table S9).Upon subclustering of the Ptprc + clusters, cDCs (cDC1: Clec9a, Xcr1, Batf3; cDC2: Clec10a, Tmem176b [87][88][89] ), which were mixed with other immune populations in the TLH, formed a separate subcluster indicating a higher resolution result (Figures S25, S26, and S33).Analysis of subcluster 5 (77 cells) (Figure S33) revealed enriched expression of recently recruited monocyte/macrophage markers Cst3 and Cd74, as well as cross-presenting DC markers Xcr1, Clec9a, and Tlr3.90 When looking at expression of these DC markers in our uniform manifold approximation and projections (UMAPs), we see that a subpopulation on the right side of this subcluster had enriched expression of cDC1 genes (Xcr1, Clec9a 34 ) and the subpopulation on the left is enriched in cDC2 markers (Clec10a, 34 Tmem176b 87,89 ), suggesting that this subcluster may contain a mixture of cDC1-and cDC2-like cells.
Comparison of previously published mouse liver data with the rat single-cell atlas indicates high consistency of the majority of the cell types between these two species (Figure 6I).We also attempted to determine if we could capture the strain-specific factors identified based on the TLH scRNA-seq map (Figure S34) in the scRNA-seq immune-enriched map.Similar to the snRNA-seq validation, we selected the top 10 genes which represented each factor and evaluated their enrichment pattern within the immune-enriched map.In line with our previous predictions, both varimax PC5 and 15 signatures show strain differences within the immune-enriched samples and are specific to hepatocyte and myeloid populations, respectively.Using immunohistochemistry, we then examined if the presence of infiltrating T cells (CD3, CD8) correlates with the differences in inflammatory potential.A periportal-biased presence of T cells was detected, but no significant frequency differences between strains were observed (Figure S12).In summary, the immune-enriched map represents a more detailed evaluation of the immune landscape of the healthy rat liver and provides additional information on B cells, DCs, Cd3 + T cells, and NK-like populations in comparison to the TLH scRNA-seq map.

Validation of computationally inferred strain-specific inflammatory differences
To functionally validate the computationally identified strain-specific differences in the inflammatory potential of hepatic myeloid cells, we performed ex vivo LPS stimulations followed by intracellular cytokine staining.In these assays, we LPS-stimulated fresh non-parenchymal cells separated by differential centrifugation from flushed, enzymatically dissociated LEW and DA rat livers.We examined cytokine secretion from tissue-resident myeloid cells via intracellular cytokine staining for tumor necrosis factor alpha (TNFa) (see STAR Methods; Figure S35).We found a higher frequency of LEW intrahepatic myeloid cells (CD45 + CD68 + CD11b + ) secreting TNFa in response to LPS stimulation compared to DA liver-resident myeloid cells (Figures 7A-7C), which suggests, in agreement with the computational findings (Figures 5C-5F), that the inflammatory potential of the hepatic myeloid cells in LEW rats (% TNFa positive = 35.25 G 3.18 (SEM)) is higher than that of DA rats (%TNFa positive = 22.25 G 1.45 (SEM)).However, despite the overall higher per-cell TNFa response in LEW myeloid cells, the overall difference in the TNFa + mean fluorescence intensity (MFI) did not reach significance (Figure 7D).In the computational analysis, the higher inflammatory potential of LEW liver myeloid cells was accompanied by the relative enriched expression of Itgal transcripts (Figure 4H), which corresponds to the protein Integrin Subunit Alpha L (ITGAL).ITGAL is a component of Lymphocyte function-associated antigen 1 (LFA-1), the expression of which is associated with inflammation and several autoimmune conditions. 91Further examination of the post-stimulation intracellular cytokine data revealed that the strain-specific pro-inflammatory differences rested primarily within Figure 6.An immune-enriched scRNA-seq rat liver map provides a higher resolution of lymphocytes and myeloid populations (A) UMAP projection of immune-enriched samples where cells that share similar transcriptome profiles are grouped by colors representing unsupervised cell clustering results.As opposed to the total liver homogenate map, B cells and plasmacytoid dendritic cells (pDCs) have been well-captured in the immuneenriched map, and Cd3+ and NK-like cells form distinct populations.The legend indicates the unique color representing the cell-type annotation of each cluster.The cluster number is shown within the curved brackets.(B) Expression distribution of Ptprc, a general immune cells marker, over UMAP projection of total liver homogenate cells.(C) Ptprc expression over UMAP projection of immune-enriched map's cells.Comparison with Figure 4B indicates that the immune-enriched map provides a better representation of the immune population compared to the total liver homogenate map.(D) Bar plot indicating the relative contribution of input samples to each cell population.Both samples have been represented in each of the clusters (cell types).(E) Labeling UMAP projection of cells based on the input sample indicates that cells from different samples have been well-integrated and clusters represent celltype differences rather than sample-specific variations.(F) The number of cells in each major population colored by the contribution of each input sample.(G) Dot plot indicating the relative expression of marker genes in each population.The size of the circle indicates the percentage of cells in each population which express the marker of interest.(H) Comparison of total liver homogenate and immune-enriched rat liver maps.Rows and columns of the correlation heatmap represent the clusters within total liver homogenate and immune-enriched maps, respectively.The color of the heatmap cells indicates Pearson correlation values between the cluster average gene expressions.The top 500 highly variable genes in each map were used for correlation calculation.The dotted box indicates that the total liver homogenate map's cluster 9 has a high correlation with B cells, pDCs, myeloid cells, and Cd3+ cell population of the immune-enriched map, confirming that immune-enriched map provides a higher resolution of lymphocytes and myeloid cells.Non-immune cell types of the two maps are consistent.(I) Comparison of rat healthy liver immune-enriched map and mouse healthy liver map [https://pubmed.ncbi.nlm.nih.gov/33106666/].The rows and columns of the correlation heatmap represent the rat and mouse clusters, respectively.The color of the heatmap cells indicates Pearson correlation values between the cluster average gene expressions.The one-to-one orthologs in the top 2,000 highly variable genes of the two maps were used for correlation calculation (see STAR Methods).The comparison indicates a high consistency between the gene expression pattern of hepatic cell types between rats and mice.ITGAL + myeloid cells, reflecting bioinformatic analysis that the LEW liver possesses a more inflammatory CD68 + CD11b + myeloid population (Figures 5, 7E, and 7F).We also observe a lack of strain-specific differences in the frequency of either CD68 + ITGAL + or CD68 + myeloid cells in the flow cytometry analysis (Figure S36).This finding is consistent with previous studies showing that DA liver myeloid cells exhibit less inflammatory characteristics, and a muted ability to stimulate T cell proliferation in comparison to LEW myeloid cells in mixed lymphocyte alloreaction assays. 92To expand on the characterization of myeloid function in these strains, CD68 + magnetic bead-based myeloid cell purification was performed on three pairs of LEW and DA rat liver TLH cell suspensions (Figures S37 and S38).The pro-inflammatory cytokine production of these cells in response to a series of LPS concentrations was then measured via a multiplexed rat cytometric bead array (CBA).Although hepatic myeloid cells from both strains displayed a dose-dependent cytokine response to LPS stimulation (Figure S39), LEW myeloid cells secreted significantly more interleukin-18 (IL-18), a LEW-enriched gene in varimax PC15 (Figure 5B), compared to DA myeloid cells (Figure 7G).Moreover, inflammatory cytokines (IL-6, IL-1a, GM-CSF, CXCL1) that are regulated by TFs positively enriched in varimax PC15, such as PU.1, 93,94 Irf8, Irf1, 95 C/EBP-b, 96,97 and Stat4 98,99 (Figures 5E and 5F), are also elevated in the stimulated LEW versus DA myeloid cells (Figures 7G and S39).Examination of these strain-specific inflammatory potential differences may serve as a point of focus for further investigation of the mechanisms behind immune-regulated hepatic disease susceptibility such as hepatic neoplasia, and liver transplant rejection.

DISCUSSION
In this study, we used a multi-platform approach to create a multi-strain atlas of the healthy rat liver.This resource helps identify rat hepatic cell types and serves as a useful baseline for hypothesis generation or to identify cellular alterations in liver disease models.We identified key immune and parenchymal populations in the healthy rat liver and their marker genes and examined their zonation tendencies within hepatic lobules.We also identified in silico strain-specific differences in hepatic myeloid populations isolated from DA and LEW rats, findings which we validate using ex vivo assays.This study illustrates cellular and molecular sources that may contribute to strain differences and highlights the potential role of myeloid cells in contributing to the baseline inflammatory state in the LEW model liver compared to the DA liver.
Tissue dissociation is a major challenge in single-cell studies of the liver, as different cell types respond differently to dissociation protocol conditions. 43We mitigated this challenge by using a combination of multiple dissociation conditions, multiple single-cell mapping technologies (scRNA-seq and snRNA-seq), and spatial transcriptomics to capture populations not well represented by either technology individually.This combined approach enabled us to better capture the diverse set of liver cell types and their zonation signatures.
Matrix factorization methods, such as varimax-rotated and standard PCA, enabled us to identify cellular identity and strain-related differences within our scRNA-seq dataset, in addition to identifying lobule zonation signatures within the spatial transcriptomics data.We found that myeloid cells from LEW livers have higher inflammatory potential than those from DA livers.We demonstrated this at the transcriptional level via scRNA-seq and confirmed this with snRNA-seq, an approach which is more resistant to dissociation-induced biases.These findings were functionally verified in vitro through intracellular cytokine staining and via the measurement of secreted cytokines following LPS stimulation of DA and LEW myeloid cells.We speculate that there is a baseline higher inflammatory milieu in the LEW rats that drives the strain-specific differences in these animals.
To examine the relevance of these identified strain differences, future rat liver atlasing efforts should include disease states such as fibrosis, ischemia reperfusion injury, or transplant rejection.Longitudinal atlasing of the liver microenvironment in these scenarios will provide valuable insights into disease-promoting populations, potentially leading to new targets to limit hepatic inflammation.Our data support the notion that reprogramming hepatic myeloid cells may be an attractive avenue to target and modulate inflammation in the rat liver. 100Taken together, our transcriptomic maps of the rat liver microenvironment contribute to our understanding of the cellular basis of the rat liver function in addition to uncovering hepatic differences between rat strains.They also provide a framework to investigate new therapeutic options in this model animal, which can be ultimately transferred to humans to cure and prevent hepatic inflammation.
Figure 7.The inflammatory potential of myeloid cells found in LEW rats is greater than that found in DA rats Myeloid cell inflammatory potential was evaluated after lipopolysaccharide (LPS) stimulation of freshly isolated liver-resident non-parenchymal cells.LPS-induced TNFa secretion was measured via intracellular cytokine staining (ICS).The non-parenchymal liver cell dissociate was obtained via a gentle enzymatic perfusion process and differential centrifugation.The resulting cells were plated in 12 well plates for 3.5 h before being stimulated for 6 h under a concentration of 1 ng/mL of LPS in the presence of 1:1000 concentration of Monensin and Brefeldin.

Limitations of the study
We recognize several limitations in our study.First, increasing the scRNA-seq datasets' sample size could provide higher statistical power.We opted to use an independent snRNA-seq dataset to increase the robustness of our findings.Compared to the former, snRNA-seq is less prone to dissociation-sensitive bias and can better capture sensitive parenchymal populations such as cholangiocytes.Rat studies are generally limited by a lack of immunological tools available, which limits the scope of in vitro validation strategies.This issue can be improved upon by testing and optimizing tools from other model systems for cross-reactivity and producing rat-specific antibodies.As well, ambient RNA is a major technical issue for studying liver tissues using sc/snRNA-seq technologies and can mask liver biological signals.This background noise was prominent in rats compared to humans, possibly due to the smaller vasculature, leading to more challenging tissue dissociation.As a result, we relied on factorization approaches, such as varimax PCA, to identify and separate biological and technical signals.Current computational methods for ambient RNA removal are limited 101 and were unable to remove the technical contamination while preserving the biological signal.Improvements in ambient RNA removal methods in the future will be beneficial to liver single-cell studies.We annotate the key DE genes expressed in each cluster and acknowledge that contamination by additional populations cannot always be excluded.Refinement of cell population and strain variation annotation, including rare cell populations, is of interest for future studies.Single-cell-resolution spatial transcriptomics methods will be useful for this.The integration of single-cell sequencing assay for transposase-accessible chromatin (scATACseq) and cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq), which can capture epigenomics, transcriptomics, and protein expression, will also lead to more refined annotations of rare cell populations.Our map only includes male samples.The inclusion of female samples will be important to understand sex-related differences in the liver.Doublet detection was performed using the scDblFinder (1.10.0)package. 104Predicted doublets had a uniform distribution within the maps and were not removed (Figure S1H).

Cell clustering, differential expression, cluster annotation
Seurat's shared nearest neighbor Louvain clustering algorithm was used to cluster the cells, based on the Harmony-corrected principal components.Differentially expressed (DE) genes associated with each cluster were identified using Seurat's FindMarkers (logfc.threshold= 0, min.pct= 0, min.pct= 0, min.cells.group= 1) implementation of the non-parametric Wilcoxon rank-sum test.scClustViz 106 was incorporated into the clustering pipeline to help find the optimal clustering resolution manually, based on known cell annotations. 41Resolution 2.5 was chosen for the single nuclei map, and resolution 0.6 was chosen for both the single cell total liver homogenate and immune-enriched maps.The Ptprc + clusters of the immune-enriched map were subclustered to examine cell subtypes, and in this case resolution, 1.0 was used.Resolutions 0.4 and 1.0 were used for subclustering of the mesenchymal population of the TLH scRNA-seq and snRNA-seq maps, respectively.Endothelial populations of the snRNAseq map were subclustered using resolution 0.8.Manual cell annotation involved evaluating the top DE genes based on known markers according to the literature.Differential expression between the DA and LEW strains within the myeloid population of the snRNA-seq map was performed using a generalized linear mixed model implemented in the NEBULA (version 1.3.0)package. 82Strain and sample were modeled as fixed and random effects, respectively and library size was used as an offset.Top strain-related genes were then identified by scoring the output based on the -log10(p values)*log (Fold change) score (Figure S19).

Matrix factorization using varimax PCA
We used matrix factorization to separate out and study the hidden patterns (factors) within our scRNA-seq data, 119 which may represent factors such as cell type gene expression program or a technical factor.Matrix factorization decomposes the gene expression matrix into the product of two lower-dimension matrices: 1) the loading matrix, which defines the relationship between the genes and the factors and can be used for pathway analysis and gene expression marker discovery; and 2) the score matrix, which represents the relationship between the factors and the cells and can be used for cluster analysis and dataset visualization.Here, we used a matrix factorization method called varimax PCA 68 to identify the hidden factors within our healthy single cell RNA-seq rat liver maps, as it worked better than standard PCA.Standard PCA identifies orthogonal dimensions that capture the maximum amounts of variation in the data.Varimax PCA applies an orthogonality-maintaining rotation to the PCA loading matrix with the goal of improving the interpretability of the PCs.This higher interpretability is mathematically achieved by maximizing the variance of the squared loadings in each factor. 68arimax PCA was applied to the normalized total liver homogenate, immune-enriched and snRNA-seq gene expression matrices separately.Interpretation of the varimax factors starts with matching factors with cell clusters and known covariates of interest (e.g., strain, sex).Varimax factors were serially plotted against PC1, to create a two-dimensional plot to help visually identify whether a separation on the basis of a specific cluster or strain was evident.For instance, different distributions of DA and LEW-derived cells over the factor of interest indicate that it has captured strain-specific variations (Figures 4C and 4F).Other factors visually correlated with known cell types (Figures S13 and 4B).
We used correlation analysis and random forest (RF) binary classifiers to automate the factor interpretation process.The correlation between the average gene expression of each cell cluster and the loading scores of each varimax factor was calculated.The top 10 differentially expressed genes of each cell population (cluster) were used to calculate the Pearson correlation scores.The results were plotted as a heatmap (Figure 4B), which was used to match each cluster with one or more varimax factors with a high absolute correlation value.The resulting matched factor and cluster pairs were robust to the number of selected top DE genes (10, 20, 30, 50).A Random Forest model was used to identify the varimax factors that capture strain-specific variations.This classifier was trained to predict the strain attributes of each cell by using varimax factors as input features.Evaluating the feature importance of the trained model uncovers the most informative varimax factor to predict the strain of interest.The model was implemented by the randomForest 107 (version 4.6.14)package and evaluated using the caTools (version 1.18.2) [https://CRAN.R-project.org/package=caTools]library (Accuracy: 0.9995, Sensitivity: 0.9994, Specificity: 0.9996).The feature matrix (varimax factors) and the corresponding labels were split in a 3:1 ratio using the 'sample.split'function of the ca-Tools package into the train and test sets, and the feature importance of the trained models was assessed.The factor with the highest feature importance score (as measured by the mean decrease in Gini score 120 ) was chosen as the best-matched factor for the predicted covariate.
To deconvolve strain-associated biological variations from sample-related confounding factors, at least two samples per strain are required.Consequently, the immune-enriched map's strain-specific varimax factors were disregarded.Two strain-specific factors were identified from the total liver homogenate map.To assess whether the varimax PCs represent technical or biological signals, the correlation between each factor and three major technical covariates, including library size, number of expressed genes, and percentage of mitochondrial gene expression was calculated.All the strain-specific components indicated a near-zero correlation with these technical covariates (Figure S14C).The hepatocyte-specific strain variation captured by varimax PC5 was separately discovered in the snRNA-seq atlas (varimax PC16) after running the same matrix factorization pipeline on this dataset (Figure S17).Due to lower representation of non-inflammatory myeloid cells within the snRNA-seq map (276 cells) compared to the scRNA-seq TLH map (1668 cells), we did not re-capture a strong myeloid-specific strain related varimax factor within the snRNA-seq liver atlas.Alternatively, to further confirm the biological relevance of strain-specific factors found in the total liver homogenate map, we created a gene signature for each strain-related factor (PC5 and PC15) by selecting the top 10 positively and negatively loaded genes for each and used these to score each cell within each strain sample of the snRNA-seq and immune-enriched maps, using the UCell 108 package (version 1.0.0)(Figures S19E, S19F, and S34).

Figure 1 .
Figure 1.scRNA-seq profiling of rat liver reveals 17 distinct cell populations (A) Overview of single-cell RNA-seq pipeline, including both the experimental and analysis workflows.(B) Major steps of the standard and matrix factorization-based single-cell RNA-seq data analysis pipeline.(C) Viable cell selection for a Lewis rat liver sample (LEW-1) based on library size and mitochondrial transcript proportion shown as an example.High-quality cells were identified from the single-cell libraries having a minimum library size of 1500 transcripts and a maximum of 40% mitochondrial transcript proportion.(D) UMAP (uniform manifold approximation and projection for dimension reduction) plot of four rat samples including 2 samples from each Dark agouti (DA) and Lewis (LEW) rat strains.Cells are colored by the number of expressed genes, with lighter colors indicating higher gene counts.(E) Bar plot indicating the relative contribution of input samples to each cluster.All samples are represented in each cluster.(F) UMAP projection of cells labeled based on the input sample indicates that cells from different samples have been well-integrated and clusters represent celltype differences rather than sample-specific variations.(G) UMAP projection of four total liver homogenate rat samples (each point represents a single cell) where cells that share similar transcriptome profiles are grouped by colors representing unsupervised clustering results.The legend indicates the unique color representing the cell-type annotation of each cluster.The cluster number is shown within the curved brackets.(H) Dot plot indicating the relative expression of marker genes in each population.The size of the circle indicates the percentage of cells in each population which express the marker of interest, and the color represents the average expression value of the marker.(I) The number of cells in each major cell type population colored by the contribution of each input sample.RBC: red blood cell, PCA: principal-component analysis, DE: differentially expressed, QC: quality control, Mac: macrophage, Mo: monocyte, Endo: endothelial, Mes: mesenchymal, Hep: hepatocyte.

Figure 2
Figure2.SnRNA-seq and spatial transcriptomic profiling of the rat liver resolves hepatocyte zonation Four additional rat liver samples were added and sequenced using snRNA-seq to better characterize hepatocyte and cholangiocyte populations and verify the strain variations identified based on the scRNA-seq TLH map.(A) UMAP projection of four snRNA-seq samples where cells are colored based on cell-type annotation.(B) Heatmap representing the average gene expression of zonated genes based on spatial data within the snRNA-seq clusters.(C) Pearson correlation between the average gene expression of the genes across snRNA-seq hepatocyte clusters and the nine layers of mouse liver cells was calculated (see STAR Methods).Mouse liver layer-9 is more periportal and layer 1 is pericentral.Red represents a positive correlation, and blue represents a negative correlation.(*: p value <0.05, **: p value <0.01, ***: p value <0.001).(D) Projection of zonation signature scores, captured by PC1, across the spatial transcriptomics spots of two healthy Wistar rat liver cryosections.The top negatively loaded genes in PC2 (and PC1) of both samples are enriched in pericentral markers, and the top positively loaded genes in PC1 (and PC2) factors are enriched in periportal markers.Red and blue represent high periportal and pericentral zonation scores, respectively.The two heatmaps represent the Pearson correlation between the zonation factors PC1 and PC2 and the average expression of snRNA-seq hepatocyte clusters.Both PC1 and PC2 are positively correlated with PP-like clusters Hep0 and Hep1 and negatively correlated with CV-like clusters Hep2 and Hep3.The asterisk and triangle symbols indicate the factors used for pathway enrichment analysis.(E) Pathway enrichment analysis using GSEA (gene set enrichment analysis) to examine active cellular pathways in periportal and central venous regions of the healthy rat liver based on spatial PC1 and PC2 loadings visualized as an enrichment map.The pathways enriched in the pericentral and periportal areas are based on PC2 (asterisk) of liver cryosections-A (left) and PC1 (triangle) of liver cryosections-B (right) respectively.Each circle represents a gene ontology (GO) biological process term.The size of the circles represents the number of genes in that pathway and blue lines indicate significant gene overlap.

Figure 3 .
Figure 3. Comparisons between transcriptomic platforms and immunohistochemistry suggest zonation patterns of selected hepatic populations To provide information on the zonation of hepatic populations within hepatic lobules, key cluster markers of each cell type were examined in PC1 and PC2 of each spatial sample.(A) Dot plot indicating the relative expression of marker genes in each population of the snRNA-seq map.The size of the circle indicates the percentage of cells expressing the marker of interest, and the color represents the average expression value of the marker.

Figure 3 .
Figure 3. Continued Expression values of (B) mesenchymal marker Ecm1 (C) cholangiocyte marker Anxa4 (D) myeloid marker Cd163 (E) non-inflammatory myeloid marker Marco (F) myeloid marker Cd68.Red and dark blue indicate higher and lower expression values in each spot, respectively.(G)Representative spatial distributions of CD68 + cells in the rat liver lobule.Rectangular layers 350 um wide were drawn from the portal tract (layer 1) to the central vein (layer 10) region.Digital images were scanned at 203 magnification.The scale bar represents 100 mm in the full image and 20 mm in the enhanced area.Each rectangular layer is referred to as a region of interest (ROI).(H) Quantification of CD68 + cell densities (#CD68+ cells/layer mm2) in the liver lobule for DA and LEW rats.30 ROIs were assessed per strain across three animals.A higher number of CD68 + cells were detected near the periportal area.No significant strain-specific differences in the spatial distribution of CD68 + cells were noted.(I) Representative spatial distributions of CD163+ cells in the rat liver lobule.(J) Quantification of CD163+ cell densities (#CD163+ cells/layer mm2) in the liver lobule for DA and LEW rats.The statistical analysis reflects similar strain and zonation patterns as CD68.Statistical significance was determined using a two-way ANOVA followed by Sidak's multiple comparisons test.(*: p value <0.05, **: p value <0.01, ***: p value <0.001,****: p value <0.0001).Data are represented as mean G SEM.Each dot in H and J represents an ROI region (n = 30).ROI: region of interest, BD: bile duct, CV: central vein, PV: portal vein.
Figure 3. Continued Expression values of (B) mesenchymal marker Ecm1 (C) cholangiocyte marker Anxa4 (D) myeloid marker Cd163 (E) non-inflammatory myeloid marker Marco (F) myeloid marker Cd68.Red and dark blue indicate higher and lower expression values in each spot, respectively.(G)Representative spatial distributions of CD68 + cells in the rat liver lobule.Rectangular layers 350 um wide were drawn from the portal tract (layer 1) to the central vein (layer 10) region.Digital images were scanned at 203 magnification.The scale bar represents 100 mm in the full image and 20 mm in the enhanced area.Each rectangular layer is referred to as a region of interest (ROI).(H) Quantification of CD68 + cell densities (#CD68+ cells/layer mm2) in the liver lobule for DA and LEW rats.30 ROIs were assessed per strain across three animals.A higher number of CD68 + cells were detected near the periportal area.No significant strain-specific differences in the spatial distribution of CD68 + cells were noted.(I) Representative spatial distributions of CD163+ cells in the rat liver lobule.(J) Quantification of CD163+ cell densities (#CD163+ cells/layer mm2) in the liver lobule for DA and LEW rats.The statistical analysis reflects similar strain and zonation patterns as CD68.Statistical significance was determined using a two-way ANOVA followed by Sidak's multiple comparisons test.(*: p value <0.05, **: p value <0.01, ***: p value <0.001,****: p value <0.0001).Data are represented as mean G SEM.Each dot in H and J represents an ROI region (n = 30).ROI: region of interest, BD: bile duct, CV: central vein, PV: portal vein.

Figure 4
Figure 4. Varimax PCs capture rat hepatic cell identity signatures and strain-specific differences (A) Bar plot representing the feature importance scores (mean decrease Gini impurity) of the top 20 features (varimax factors) of the random forest model trained to predict the strain attributes of the rat hepatic cells.Varimax PC5 and 15 are the most informative features to differentiate cells of each strain from another, which indicates the two factors have captured strain-related variations within the map.(B) A correlation heatmap between the average gene expression of each cluster and the loading scores of varimax factors (capturing the contribution of all genes to a factor).Columns are varimax factors and rows are cell populations.Each cell-type cluster is defined by key marker genes, and dark red or blue indicates that the expression of a marker gene set is positively or negatively correlated, respectively, with a particular varimax factor.A high absolute correlation value indicates a match between a varimax factor and a cell-type cluster.(C) The projection of cells over varimax-1 and 5 indicates that the cells from each strain form distinct clusters over varimax-5.(D) Boxplot indicating the distribution of varimax-5 score over each strain.Cells from DA and LEW strains represent significantly different varimax-5 scores (Wilcoxon-test p value <2.2e-16), indicating that varimax-5 has captured strain differences.(E) The top 10 genes on the top (left table) and bottom (right table) of the varimax-5 loading list mainly contain known hepatocyte markers, indicating that varimax-5 has captured hepatocyte-specific strain differences.Genes with high positive scores (left table) are associated with the DA strain and genes indicating negative loading scores (right table) are LEW-related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(F) Projection of cells over varimax-1 and 15 indicates that a population of cells from each strain (dotted lines) forms distinct clusters over varimax-15.Annotation of the selected cells indicates that they are mainly from the Marco+ myeloid cluster 5. (G) Boxplot indicating the distribution of hepatic cells based on strain over varimax-15.(Wilcoxon-test p value <2.2e-16).The outlier data points (dotted lines) are mainly myeloid cells.(H) The top 10 genes with positive (right table) and negative (left table) varimax-15 loading scores are immune-response related.Genes with positive scores (right table) are associated with the LEW strain, and genes indicating negative loading values (left table) are DA related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(I) Expression pattern of known myeloid marker genes Marco, Vsig4, Cd68, and Lyz2 over UMAP.Dark green represents high expression values.The distribution of general myeloid markers (Cd68, Vsig4) and non-inflammatory myeloid marker (Marco) is consistent with the varimax-15 distribution (Figure2J).(J) The UMAP projection of cells colored based on the varimax-15 score shows the enrichment of varimax-15 over Marco + myeloid population (cluster 5).Darker colors represent higher values of varimax-15 scores.Data are represented as mean G SEM with each dot representing a single cell.Corrcoef.: correlation coefficient, Var: varimax PC. varimax PCs are referred to as PCs within the main text.
Figure 4. Varimax PCs capture rat hepatic cell identity signatures and strain-specific differences (A) Bar plot representing the feature importance scores (mean decrease Gini impurity) of the top 20 features (varimax factors) of the random forest model trained to predict the strain attributes of the rat hepatic cells.Varimax PC5 and 15 are the most informative features to differentiate cells of each strain from another, which indicates the two factors have captured strain-related variations within the map.(B) A correlation heatmap between the average gene expression of each cluster and the loading scores of varimax factors (capturing the contribution of all genes to a factor).Columns are varimax factors and rows are cell populations.Each cell-type cluster is defined by key marker genes, and dark red or blue indicates that the expression of a marker gene set is positively or negatively correlated, respectively, with a particular varimax factor.A high absolute correlation value indicates a match between a varimax factor and a cell-type cluster.(C) The projection of cells over varimax-1 and 5 indicates that the cells from each strain form distinct clusters over varimax-5.(D) Boxplot indicating the distribution of varimax-5 score over each strain.Cells from DA and LEW strains represent significantly different varimax-5 scores (Wilcoxon-test p value <2.2e-16), indicating that varimax-5 has captured strain differences.(E) The top 10 genes on the top (left table) and bottom (right table) of the varimax-5 loading list mainly contain known hepatocyte markers, indicating that varimax-5 has captured hepatocyte-specific strain differences.Genes with high positive scores (left table) are associated with the DA strain and genes indicating negative loading scores (right table) are LEW-related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(F) Projection of cells over varimax-1 and 15 indicates that a population of cells from each strain (dotted lines) forms distinct clusters over varimax-15.Annotation of the selected cells indicates that they are mainly from the Marco+ myeloid cluster 5. (G) Boxplot indicating the distribution of hepatic cells based on strain over varimax-15.(Wilcoxon-test p value <2.2e-16).The outlier data points (dotted lines) are mainly myeloid cells.(H) The top 10 genes with positive (right table) and negative (left table) varimax-15 loading scores are immune-response related.Genes with positive scores (right table) are associated with the LEW strain, and genes indicating negative loading values (left table) are DA related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(I) Expression pattern of known myeloid marker genes Marco, Vsig4, Cd68, and Lyz2 over UMAP.Dark green represents high expression values.The distribution of general myeloid markers (Cd68, Vsig4) and non-inflammatory myeloid marker (Marco) is consistent with the varimax-15 distribution (Figure2J).(J) The UMAP projection of cells colored based on the varimax-15 score shows the enrichment of varimax-15 over Marco + myeloid population (cluster 5).Darker colors represent higher values of varimax-15 scores.Data are represented as mean G SEM with each dot representing a single cell.Corrcoef.: correlation coefficient, Var: varimax PC. varimax PCs are referred to as PCs within the main text.

Figure 5 .
Figure 5. Strain-specific differences are found in intrahepatic myeloid cells (A) Expression pattern of the top DA-enriched genes (Ly6al, Cd163, Hmox1, Siglec5) over PC15 and 1. LEW and DA myeloid cells have been marked with dotted circles.Dark green indicates higher expression values.Comparison with Figure2Fconfirms that the selected genes have higher expression in the DA strain compared to LEW. (B) Expression pattern of the top LEW-enriched genes (Itgal, Il18, Ccl3, Timp2) over varimax-15 and 1.Comparison with Figure 2F confirms that the selected genes have higher expression in the LEW strain compared to DA. (C) Dot plot indicating the relative expression of strain-related genes within the myeloid fraction (clusters 5, 10, 9) of each strain.The top 10 genes with positive (LEW-associated) and negative (DA-associated) varimax-15 loading scores have been selected.The size of the circle indicates the percentage of cells in each population expressing the marker, and its color shows the average expression value.(D) Pathway enrichment analysis using GSEA (gene set enrichment analysis) to examine active cellular pathways in LEW vs. DA myeloid cells based on varimax-15 loadings visualized as an enrichment map.Each circle represents a gene ontology (GO) biological process term.The size of the circles represents the number of genes in that pathway, and blue lines indicate significant gene overlap.Since PC15 is positively correlated with the LEW strain and negatively correlated with DA, red circles represent activated pathways in LEW and blue indicates upregulated pathways in DA.No pathway was significantly upregulated in DA.Transcriptional factor (TF) binding site-based gene set enrichment analysis using gProfiler on the ChEA ChIP-Seq database identifies TFs which may be activated in (E) DA and (F) LEW myeloid cells.TFs are sorted based on their enrichment significance calculated as -log10(adjusted p value).Dark purple indicates higher significance.Purple boxes highlight TFs which are uniquely enriched in that strain.
Figure7.The inflammatory potential of myeloid cells found in LEW rats is greater than that found in DA rats Myeloid cell inflammatory potential was evaluated after lipopolysaccharide (LPS) stimulation of freshly isolated liver-resident non-parenchymal cells.LPS-induced TNFa secretion was measured via intracellular cytokine staining (ICS).The non-parenchymal liver cell dissociate was obtained via a gentle enzymatic perfusion process and differential centrifugation.The resulting cells were plated in 12 well plates for 3.5 h before being stimulated for 6 h under a concentration of 1 ng/mL of LPS in the presence of 1:1000 concentration of Monensin and Brefeldin.(A) Flow cytometry plots showing the gating strategy for macrophages.(B) Percentage of TNFa + secreting CD68 + CD11b + myeloid cells in the unstimulated control and stimulated conditions of Dark Agouti and Lewis macrophages.(C) Summary graphs of Lewis versus Dark Agouti total TNFa as a percentage of CD68 + CD11b + myeloid, (D) and of the mean fluorescence intensity (MFI) of Lewis vs. Dark Agouti TNFa.(E) Representative flow cytometry plot of TNFa secretion patterns based on ITGAL subpopulations.(F) and summary graph ITGAL expressing CD68 + CD11b + myeloid subpopulations.Plotted are the values from all 4 experimental replicates.Statistical significance for ICS was determined using a non-parametric 2 tailed Mann-Whitney test.(n = 4) (G) Cytometric bead array (LEGENDplex) was performed to quantify the level of cytokines (TNFa, Il-18, CXCL1) on culture supernatants of enriched CD68 + myeloid cells after 24 h of stimulation in various LPS concentration conditions (0, 0.05, 0.1, 1, 10 ng/mL).Three technical replicates were used per animal.Statistical significance of the CBA was determined using a two-way ANOVA and Sidak's multiple comparisons test (n = 3) Data are represented as mean G SEM with each dot representing a single animal.(*: p value <0.05, **: p value <0.01, ***: p value <0.001,****: p value <0.0001) DA: dark agouti, LEW: lewis, SSC-A: side scatter area, FSC-A: forward scatter area.

(
Continued on next page) . (Wilcoxon-test p value <2.2e-16).The outlier data points (dotted lines) are mainly myeloid cells.(H) The top 10 genes with positive (right table) and negative (left table) varimax-15 loading scores are immune-response related.Genes with positive scores (right table) are associated with the LEW strain, and genes indicating negative loading values (left table) are DA related.The absolute loading scores indicate the contribution of each gene to the corresponding factor.(I) Expression pattern of known myeloid marker genes Marco, Vsig4, Cd68, and Lyz2 over UMAP.Dark green represents high expression values.The distribution of general myeloid markers (Cd68, Vsig4) and non-inflammatory myeloid marker (Marco) is consistent with the varimax-15 distribution (Figure