Single-cell spatial reconstruction reveals global division of labour in the mammalian liver

Abstract

The mammalian liver consists of hexagon-shaped lobules that are radially polarized by blood flow and morphogens. Key liver genes have been shown to be differentially expressed along the lobule axis, a phenomenon termed zonation, but a detailed genome-wide reconstruction of this spatial division of labour has not been achieved. Here we measure the entire transcriptome of thousands of mouse liver cells and infer their lobule coordinates on the basis of a panel of zonated landmark genes, characterized with single-molecule fluorescence in situ hybridization. Using this approach, we obtain the zonation profiles of all liver genes with high spatial resolution. We find that around 50% of liver genes are significantly zonated and uncover abundant non-monotonic profiles that peak at the mid-lobule layers. These include a spatial order of bile acid biosynthesis enzymes that matches their position in the enzymatic cascade. Our approach can facilitate the reconstruction of similar spatial genomic blueprints for other mammalian organs.

Extended Data Figure 1. Low magnification images of the six landmark genes.

smFISH examples of our 6 landmark genes - the pericentral genes Glul and Cyp2e1 and the periportal genes Ass1, Asl, Alb and Cyp2f2. Bright cells have high mRNA content. Scale bar-100μm. CV-central vein, PN-portal node. Micrographs are representative of at least 10 lobules and at least two mice per gene.

Extended Data Figure 2. Single cell RNAseq reveals three distinct liver cell populations.

(a-c) t-SNE plots. Each dot is a cell colored according to the aggregated expression of (a) the hepatocyte marker genes Apoa1, Glul, Acly, Asl, Cyp2e1, Cyp2f2, Ass1, Alb, Mup3, Pck1, G6pc, (b) the Kupffer cell marker genes Irf7, Spic, Clec4f and (c) the endothelial cell markers Ushbp1, Myf6, Oit3, Il1a, F8, Bmp2, C1qtnf1, Mmrn2, Pcdh12, Dpp4. (d) smFISH of a liver lobule demonstrating the antagonistic expression of Glul (red dots), and Ass1 (green dots). Scale bar-20μm. Inset of highlighted area demonstrates the minimal co-expression of these two genes. Scale bar-5μm. Blue is dapi-stained nuclei. CV – central vein, PN – portal node. Micrographs are representative of at least 10 lobules and at least two mice per gene. (e) Ass1 and Glul are mutually exclusive in the single cell RNAseq data. Each dot is a cell, X-axis is the expression of Glul in fraction of total cellular UMI, Y-axis is the expression of Ass1.

Extended Data Figure 3. Prior and posterior probability computation of the scRNAseq data.

(a-b) Zone-dependent probabilities of observing different expression levels for each of the six landmark genes. Horizontal lines denote the UMI observed for each gene for cell 609 (a) and cell 629 (b). The posterior probabilities for each layer suggest that cell 609 most probably originated from the periportal layer 8, whereas cell 629 originated from the pericentral layer 1. Note that the posterior probabilities incorporate the lower sampling of cell 609, which had a total of 3,358 UMI, compared to cell 629, which had 9,913 UMI. (c-d) Prior probability of observing a hepatocyte from each of the 9 layers. (c) Hexagonal lobule geometry and 9 concentric circles at constant radii increments spanning the central vein and portal node. (d) Area of each layer is the intersection with the hexagon.

Extended Data Figure 4. Sensitivity of reconstruction to number and features of landmark genes.

(a) Zonation properties of our landmark genes. –log10 of the Wilcoxon rank-sum p-values for the comparison of smFISH expression measurements in cells in sequential lobule layers for each of the 6 landmark genes. Stars denote layer pairs with significant changes in expression (p-value<0.05). (b) Mean reconstruction accuracy for different subsets of landmark genes (defined as 1 minus the mean of the sum-squared differences between the profiles predicted using each subset and the profiles obtained when using the entire panel of landmark genes). Dot color represents the number of landmark genes in each subset. (c) Reconstruction accuracy approaches saturation when using 6 landmark genes. G – Glul, Cf- Cyp2f2, Ce – Cyp2e1, As – Ass1, Al – Asl, Ab - Alb. Dotted line connects the most accurate partial subset of landmark genes for each panel size. (d) The contribution of each landmark gene to zonation reconstruction is strongly dependent on its spatial non-uniformity among different layers but less on its intra-layer variability. X-axis shows the entropy of the summed-normalized average profile of each landmark gene. Y-axis shows the average among all layers of the coefficient of variation of the intra-layer expression levels. The score for each gene is the average of the ratio in reconstruction error with and without the gene among all combinations that include the gene. High scores indicate that the landmark gene strongly improves reconstruction when added to combinations of other landmark genes.

Extended Data Figure 5. Validation of reconstructed zonation profiles using smFISH.

(a) Reconstructed zonation profiles based on the scRNAseq data (blue) and smFISH measurements (red) for our landmark genes. (b) Validation of reconstructed profiles for 20 genes not used for the inference algorithm. Profiles are normalized by the mean, patches show s.e.m. Note that Pck1 has a broader profile compared to other studies, since our measurements were done on fasted mice, a metabolic state in which gluconeogenesis becomes substantial in pericentral layers.

Extended Data Figure 6. Reconstructed zonation profiles capture a wide dynamic range.

Reconstructed zonation profiles of the pericentral gene Oat (magenta) and the pericentral progenitor marker Axin2 (red), and the periportal urea cycle enzyme gene Arg1 (green) and the periportal progenitor marker Sox9 (blue). Dashed box highlights a blow-up of Axin2 and Sox9, genes with 100-fold lower expression than Oat and Arg1. Expression values are the estimated fraction of the total cellular mRNA molecules. Lobule layers are numbered from the central vein (CV, layer 1) to the portal node (PN, layer 9). The slightly broader zonation profile of Sox9 mRNA compared to the Sox9 signal observed in using Sox9-GFP or Sox9-CreERT;R26RtdTomato may stem from differences in the dynamic range of the detection methods or from differences in the mouse metabolic states (fasted vs. ad-libitum).

Extended Data Figure 7. Porto-central ratio of reconstructed zonation profiles correlates with previous studies.

(a) Correlation between the pericentral bias computed from our data and the one from Braeuning et al.. Y axis is our pericentral bias (PC/PP), computed as the ratio between the median zonation profiles in layers 1-4 and the median in layers 6-9. X-axis is the ratio between expression in pericentrally and periportally enriched genes from. 79 of the 88 genes considered pericentral according to Braeuning et al. (red circles) are also pericentrally biased in our dataset (hypergeometric p<4E-9). 70 of the 81 genes considered periportal according to Braeuning et al. (green circles) are also periportally biased in our dataset (hypergeomtric p<1E-16). Black squares mark genes that we found to be significantly zonated. (b) Pericentral bias computed by our method correlates with previous bulk RNAseq studies that compared pericentral and periportal populations isolated via laser capture microdissection. RSpearman=0.74, p<E-80. Black dots represent genes considered significantly zonated Data in. (a-b) include genes with average expression across cells which is higher than 5E-6 of the total UMI counts. (c) Reconstructed zonation profiles for the genes found in Wang et al. to have higher expression in Axin2+ pericentral hepatocytes (qval<0.2). (d) Reconstructed zonation profiles for the genes found in Wang et al. to have lower expression in pericentral Axin2+ hepatocytes (qval<0.2).

Extended Data Figure 8. Zonation profiles of genes that significantly increased or decreased when perturbing signaling pathways.

(a-b) Zonation profiles of liver genes previously shown to increase (a) or decrease (b) in liver of mice exposed to chronic hypoxia (PO2 =11.5kPa, compared to PO2 =18.0kPa.20). (c-d) Zonation profiles of genes shown to significantly increase (c), log-fold>2) or decrease (d), log fold<-2) in expression in liver tumors harboring an activating Ha-ras mutation (FDR-corrected p-value<0.05). (e-f) Zonation profiles of genes that have higher (e, fold>1) or lower (f, fold<1)expression in hypopituitary mice (FDR-corrected p-value<0.05).

Extended Data Figure 9. KEGG pathways enriched for zonated genes.

(a) Average max-normalized zonation profiles of KEGG pathways enriched for zonated genes. Panel displays all pathways with more than 10 genes and hypergeometric q-value<0.1 (Table S5 exhibits full results). Each profile was normalized by its maximum along all layers. (b) Periportal bias of liver secreted proteins. Genes are based on an annotated liver secretome. (c) Pericentral bias for liver detoxification genes. Shown are genes for cytochrome P450, Uridine 5'-diphospho-glucuronosyltransferase and glutathione S transferase. Images in (b-c) include significantly zonated genes from each pathway (Kruskal-Wallis qval<0.01 and more than 5*10^-5 of the total cellular UMIs on average in at least one of the 9 layers), each profile normalized to maximum of 1.

Extended Data Figure 10. Non-monotonic zonation profiles of liver genes.

(a) Max-normalized zonation profiles of the concise set of 46 non-monotonic genes (Methods). (b) Mup3 is highly variable but exhibits a peak at layer 7 and a decreased expression in both the pericentral layer 1 and periportal layer 9. (c) Igfbp2 exhibits a non-monotonic zonation profile, peaking at layers 3-6. Scale bar-20µm. CV – central vein. Enlarged rectangles in (b-c) mark the periportal layers (blue), mid-lobule layers (green) and pericentral layers (magenta). Scale bar-8µm. Micrographs are representative of at least 10 lobules and at least two mice per gene. (d) Igfbp2, Igfbp1 and Igfbp4 are downstream to Igf1, the secreted protein that they bind and stabilize. (e) Non-monotonic zonation profiles observed do not stem from ploidy-specific regulation. Each dot represents a pair of adjacent cells (within 30µm of each other), x-axis is the higher ploidy cell, y-axis is the lower ploidy cell. Pvalues are Wilcoxon sign-ranked tests. Gene expression for Igfbp2, Hamp, and Cyp8b1 was quantified as the number of smFISH dots divided by the segmented cell volumes; Mup3 expression was quantified as the average intensity, due to its higher abundance. Quantification for each gene based on at least 60 pairs. (f) Immunohistochemistry of Cyp8b1 demonstrates higher protein concentration in layers 2-3 (dashed curve) compared to layer 1 (dotted curve). Scale bar-50μm. Micrographs are representative of 7 lobules in 2 mice.

Fig. 1. Strategy for spatially resolved single-cell reconstruction of liver zonation.

(a) Spatial barcode of zonated landmark (Lm) genes measured with smFISH. (b) scRNAseq provides the transcriptome of thousands of mouse liver cells. (c) Spatial barcode is used to infer the porto-central coordinates of each cell. For example, cell B, containing 50 UMI of the blue Lm1, and 4 UMI of the red Lmn is inferred to have resided in the outer periportal layers. Cell A, with 2 UMI of Lm1 and 34 UMI of Lmn is inferred to have resided in the inner pericentral layers. (d) These coordinates are used to reconstruct the spatial zonation profiles of all liver genes.

Fig. 2. A spatial barcode of zonated landmark genes.

(a) smFISH micrographs of our landmark genes. Gray dots - single mRNA, Red - phalloidin-stained membranes. Insets show pericentral regions of the spatially-discordant genes Glul (top) and Ass1 (bottom), demonstrating the dynamic range of smFISH. Scale bar-20µm. CV-central vein. Micrographs are representative of at least 10 lobules and two mice per gene. (b) Zonation profiles of the landmark genes. X axis is the scaled distance from the CV, Y axis is the max-normalized expression level. All landmark genes are highly significantly zonated (Kruskal-Wallis p<10^-98). Errorbars are s.e.m. based on at least 800 cells from 10 lobules and two mice.

Fig. 3. Genome-wide zonation patterns revealed by spatially-resolved single-cell RNAseq.

(a-b) t-SNE visualization of the scRNAseq data. Colors denote relative cellular expression of landmark genes Cyp2f2 (a) and Cyp2e1 (b). (c) Zonation profiles of the 3496 significantly zonated genes. Genes sorted by the profiles’ centers of mass. Left-examples of pericentral, periportal and non-monotonic profiles. (d-g) Signaling pathways affecting liver zonation. Shown are average zonation profiles, patches are s.e.m. (d) Green-Wnt-activated genes (n=810), red-Wnt-inhibited genes (n=193). (e) Green-hypoxia activated genes (n=95), red-hypoxia inhibited genes (n=45). (f) Green-Ras-activated genes (n=43) red-Ras-inhibited genes (n=50). (g) Green-genes reduced (n=81), red-genes increased (n=51) in expression in hypopituitary mice.

Fig. 4. Non-monotonic zonation and sequential order.

(a) Hamp, encoding the hormone hepcidin peaks in mid-lobule layers. Scale bar–20µm. CV–central vein, PN–portal node. Insets mark periportal (blue), mid-lobule (green) and pericentral (magenta) layers, inset scale bars–8µm. Micrograph representative of 12 lobules from three mice. (b) Spatial order of enzymes matches their position in the classic bile-acid biosynthesis cascade. Cyp7a1 and Hsd3b7 peak at layer 1, Cyp8b1 peaks at layers 2-3. Patches are s.e.m., based on 1,000 bootstrap iterations (Supplementary information). (c) Cyp8b1 peaks in layer 2. Dashed line-pericentral layer. Scale bars–20µm, inset scale bar-5µm. Micrographs are representative of 12 lobules and three mice.