Space-time logic of liver gene expression at sub-lobular scale

Abstract

The mammalian liver is a central hub for systemic metabolic homeostasis. Liver tissue is spatially structured, with hepatocytes operating in repeating lobules, and sub-lobule zones performing distinct functions. The liver is also subject to extensive temporal regulation, orchestrated by the interplay of the circadian clock, systemic signals and feeding rhythms. However, liver zonation has previously been analysed as a static phenomenon, and liver chronobiology has been analysed at tissue-level resolution. Here, we use single-cell RNA-seq to investigate the interplay between gene regulation in space and time. Using mixed-effect models of messenger RNA expression and smFISH validations, we find that many genes in the liver are both zonated and rhythmic, and most of them show multiplicative space-time effects. Such dually regulated genes cover not only key hepatic functions such as lipid, carbohydrate and amino acid metabolism, but also previously unassociated processes involving protein chaperones. Our data also suggest that rhythmic and localized expression of Wnt targets could be explained by rhythmically expressed Wnt ligands from non-parenchymal cells near the central vein. Core circadian clock genes are expressed in a non-zonated manner, indicating that the liver clock is robust to zonation. Together, our scRNA-seq analysis reveals how liver function is compartmentalized spatio-temporally at the sub-lobular scale.

Extended Data Fig. 1. scRNA-seq pre-processing.

(a) Histogram of number of UMIs per cell barcode for each mouse. Red patches mark the cells used for background estimation (100-300 UMI/cell barcode), gray patches mark the cells used for downstream analysis (1000-10000 UMI/cell barcode). (b) Histogram of fraction of all UMIs mapping to mitochondrial genes. Filter used for downstream analysis in grey (0.09-0.35). (c) smFISH staining of a liver lobule with probes against Cyp2e1 (red) and Cyp2f2 (green). CV = central vein, PN = portal node. Overall, the data combine 10 images from from two mice. (d) Expression of Cyp2e1 and Cyp2f2 in cells with different fraction of mitochondrial expression. Three different filters for the fraction of UMIs mapping to mitochondrial genes (0-0.09, 0.09-0.35, 0.35-1) were applied, the data of all mice merged and the resulting datasets visualized as t-SNE plots. (e) Violin plots for the correlations between Cyp2e1 and Cyp2f2 expression in single hepatocyte populations with different filters for fractions of mitochondrial expression. Each dot represents one mouse (n=10 mice for each distribution) and the shape of the violin represents the density of points. (f) Comparison of the zonation profiles of Z and Z+R genes obtained in our current study and the previous reconstruction from Halpern et al. 2017. Profiles were interpolated to fit 15 layers, where 1 is pericentral and 8 is periportal. Dots indicate the center of mass (expression-weighted lobule layer) of the Z and Z+R genes computed in both datasets, for gene having an average expression of at least 10^-5 in Halpern et al. r is the correlation coefficient, p the corresponding p-value from a standard linear regression.

Extended Data Fig. 2. Log-transformed reconstructed profiles, pre-filtering of the genes and comparison with external datasets.

(a-c) Expression levels of the reconstructed profiles for the genes from Figure 1F-H after log-transformation (Methods). Dots in represent data points from the individual mice. Lines represent mean expression per time point. Shaded areas represent one standard deviation (SD) across the mice (n=2 or n=3 depending on the time point, Methods). (d) Biological variability of gene profiles across independent replicate liver samples, quantified in terms of the average relative replicate variance. 0 shows perfectly reproducible profiles while 1 the most variable genes (Methods). Genes inside the bottom-right box (x-cutoff at 10^-5; y-cutoff at 0.5) are selected and contain all but one of the reference genes. Colored dots show reference zonated genes (blue) and reference rhythmic genes (orange). (e) Comparison of the peak times for rhythmic genes in R and Z+R, with the dataset from Atger et al, PNAS 2015. Circular correlation coefficient is 0.746 (Methods). (f) Boxplot of the mRNA half-lives (data from Wang, J. et al, 2017) shows that R genes as a group (median, orange line) are the shortest-lived. Box limits are lower and upper quartiles, whiskers extend up to the first datum greater/lower than the upper/lower quartile plus 1.5 times the interquartile range. Remaining points are marked. MW stands for the two-sided Mann-Whitney test, and KW stands for Kruskal-Wallis test.

Extended Data Fig. 3. Z+R and ZxR transcripts with corresponding rhythmic protein accumulation in bulk mass spectrometry data.

(a-b) Rhythmic proteins corresponding to Z+R (a) and ZxR (b) transcripts were selected from Robles et al., 2014 (from original Table S2), and fitted with harmonic regression (p-value of rhythmicity from F-tests are indicated above the plot). Only proteins having a p-value<0.01 are shown. (c) Scatter plot of the phase of the fits from the transcripts (x-axis) against the phase of the fits from the proteins (y-axis). The diagonal is indicated with a dashed grey line, the theoretical upper bound (6h) for the delay between mRNA and protein is indicated with a dashed red line. All rhythmic proteins (q<0.2 in the original analysis) are represented.}

Extended Data Fig. 4. Additional validations for the Z+R category.

(a-b) smFISH (Stellaris, Methods) for Elovl3 (Z+R) and Arg1 (Z+R). smFISH quantifications were made for ZT0 and ZT12 (Methods). Left: representative images at ZT0, ZT12 for Elovl3 (a) or Arg1 (b). Pericentral veins (CV) and a periportal node (PN) are marked. Scale bar - 20μm. Right: quantified profiles for each gene at the two time points from smFISH (top, line plot is the mean number of mRNAs, shaded area indicate SD across twelve images), and scRNA-seq data (bottom, line plot is mean expression, shaded area is SD across mice, n=2 or n=3 depending on the time point, Methods).

Extended Data Fig. 5. The core circadian-clock is not zonated.

(a) Spatial and temporal profiles and fits for circadian core-clock genes. Peak times are indicated on the temporal representation. For the genes Cry1 and Clock, additional dashed lines represent fits for the R model, as the Schwartz BIC weights from the R and Z+R models were close (Supplementary Table 2). (b) Amplitudes and peak times of the core-clock circadian genes in a polar coordinate representation (clock-wise ZT times are indicated, distance from the center corresponds to the amplitude) show the expected organization of core clock transcript in the liver.

Extended Data Fig. 6. Spatio-temporal mRNA expression profiles for heat-responsive genes. Represented genes correspond to the top 200 targets from the ChIP-Atlas list for mouse HSF1.

(a) Polar plot representation of the transcripts that are R, Z+R, or ZxR among the HSF1 targets from the ChIP-Atlas (http://chip-atlas.org/). Genes involved in chaperone functions (chaperones, co-chaperones, or chaperone facilitators) are named. Color indicates zonation, while grey dots show purely rhythmic genes. (b) Spatial representation of the transcripts corresponding to proteins involved in chaperone functions, separated in central (left, cytoplasmic function) and portal (right, endoplasmic reticulum function) zonation.

Extended Data Fig. 7. Rhythmicity of Wnt targets in bulk RNA-seq, and proteomics liver time series data.

(a) Enrichment of rhythmic genes (R, Z+R and ZxR) among the targets of the Wnt pathway, computed on the bulk dataset (Atger et al., 2015). Targets above a given percentile (x-axis) of Apc-KO fold change are considered. The percentage of rhythmic genes in the whole Atger et al. dataset is indicated by a dashed blue line. (b) Bulk mRNA (coming from Atger et al. dataset) rhythmicity profiles of Wnt targets among the top-50 targets with highest Apc-KO fold change. Gene profiles are centered around their mean. An enrichment of the phases around ZT8-14 is observed, in agreement with Figure 6a. (c) Polar plot representation of the individual gene phases and amplitudes represented in panel b (bulk data). (d) Temporal representation of selected genes profiles from the scRNA-seq (top, n=2 or n=3 animals depending on the time point, Methods) and bulk proteomics (bottom, data from Robles et al, 2014, n=2 replicates per time point sampled every 3h) data. Represented profiles are ones with (1) the highest Apc-KO fold change, (2) a significantly rhythmic protein (p < 0.05, standard harmonic regression, F-test), and (3) belonging to the Z+R or ZxR category. (e) Enrichment/depletion at different times (window size: 3h), of both positive and negative Ras (N=31 and N=33, respectively) and Hypoxia (N=73 and N=41, respectively) targets (background: all R and Z+R genes). Colormap shows p-values (two-tailed hypergeometric test): red (blue) indicates enrichment (depletion).

Figure 1. A scRNA-seq approach to space-time gene expression in mouse liver.

(a-e) Global gene expression varies in both space and time, as shown using t-SNE visualizations of the scRNA-seq (n=19663 hepatocytes examined over ten independent animals). Each dot represents one cell. Individual cells are colored by the (a posteriori assigned) lobule layer (a), Zeitgeber time (b), expression levels of the zonated genes Cyp2f2 and Cyp2e1 (c-d), or the temporally-regulated and centrally zonated gene Elovl3 (e). (f-h) Reconstructed spatial profiles (lobule layers 1-8) of selected zonated genes (f, top: Glul pericentrally (PC) expressed, bottom: Ass1 periportally (PP) expressed); rhythmic but non-zonated genes (g, top: core-clock gene Bmal1 peaking at ZT18-0, bottom: clock-controlled Dbp, peaking at ZT6-12); zonated and rhythmic genes (h, Top: Elovl3, bottom: Pck1). Expression levels correspond to fraction of total UMI per cell in linear scale. Log-transformed profiles are in Extended Data Figure 1. Dots in f-h represent data points from the individual mice. Lines represent the mean expression per time point. Shaded areas represent one standard deviation (SD) across the mice. For the scRNA-seq we used n=2 (ZT6, ZT18) or n=3 mice (ZT0, ZT12) (Methods).

Figure 2. Space-time mRNA expression profiles categorized with mixed-effect models.

(a) Spatial and temporal variations for mRNA transcript profiles, calculated as standard deviations (SD) of log2 expression along spatial or temporal dimensions. Colored dots correspond to reference zonated genes (orange) and reference rhythmic genes (blue) (Methods). (b) Extended harmonic regression model for spatio-temporal expression profiles describing a static but zonated layer-dependent mean μ(x), as well as layer-dependent harmonic coefficients (a(x) and b(x)). All layer-dependent coefficients are modeled as second order polynomials; i denotes the biological replicates. Temporal dependency is modelled with 24h-periodic harmonic functions. μ_i are random effects needed due to the data structure hierarchy (Methods). (c) Schema illustrating the different categories of profiles. Depending on which coefficients are non-zero (Methods), genes are assigned to: F/N (flat or noisy, not represented), Z (zonated), R (rhythmic), Z+R (additive zonation and rhythmicity), ZxR (interacting zonation and rhythmicity). Graphs emphasize either zonation (top), with the x-axis representing layers, or rhythmicity (bottom), with the x axis representing time (ZT). Right side of the panel: two examples of fits (Elovl3 and Cyp7a1, respectively Z+R and ZxR). (d) Number of transcripts in each category. (e) Boxplot of the mean expression per category shows that zonated genes (Z, Z+R and ZxR) are more expressed than rhythmic (R) or flat/noisy (F/N). ZxR genes are the most expressed according to median expression (orange line). Box limits are lower and upper quartile, whiskers extend up to the first datum greater/lower than the upper/lower quartile plus 1.5 times the interquartile range. The number of genes per category is indicated. Remaining points are outliers. KW stands for the Kruskal-Wallis and MW for the Mann-Whitney (two-sided) tests.

Figure 3. Properties of dually zonated and rhythmic mRNA profiles.

(a) Proportions of pericentral (green) and periportal (blue) transcripts are similar in Z and Z+R. Mid-lobular genes (orange) are rare (<2%). The number of genes are N = 484, 14, 628, and 157, 6, 184 for the portal, mid-lobular, and central genes in Z and Z+R, respectively. KS stands for the two-sided Kolmogorov-Smirnov test. (b) Peak time distributions of rhythmic transcripts are similar in R and Z+R categories (two-sample Kuiper test). (c-d) Effect sizes of zonation (slope) vs. rhythmicity (fold-change amplitude, log2 peak-to-trough) in Z+R genes (c). Magnitude of time shifts (delta phase, in hours) vs. fold-change amplitude gradient (delta amplitude, in log2) along the central-portal axis in ZxR genes (d). Genes for which the protein is also rhythmic (p<0.05, harmonic regression, F-test) in bulk data are indicated in dark red with the corresponding label (represented in Extended Data Figure 3).

Figure 4. smFISH analysis of rhythmic and zonated transcripts.

(a) smFISH (RNAscope, Methods) of the core clock genes Bmal1 and Per1 (both assigned to R) in liver slices sampled every 3 hours. Left: representative images at ZT0, ZT06, ZT12 and ZT18 for Bmal1. Central vein (CV) and a portal node (PN) are marked. Scale bar is 50μm. Endothelial cells lining the PC and cholangiocytes surrounding the PP were excluded from the quantification. mRNA transcripts and nuclei were detected in PN and PC zones (Methods). Right: temporal profiles of Bmal1 and Per1 from smFISH at 8 time points from ZT0 to ZT21, every 3 hours (the line shows the mean number of quantified mRNAs, shaded area indicate SD across 5-8 images from one animal per time point), and scRNA-seq (same representation as in Figure 1f-h, PC is layer 1, PN is layer 8), both in PN and PC regions. (b-c) smFISH (Stellaris, Methods) for Pck1 (Z+R) and Acly (ZxR) in liver slices sampled every 12 hours. Left: representative images at ZT0, ZT12 for Pck1 (b) or Acly (c). Central vein (CV) and a portal node (PN) are marked. Scale bar - 20μm. Right: quantified profiles for each gene in the two time points from smFISH (the line shows the mean number of quantified mRNAs, shaded areas indicate SD across at least 10 independent images taken from at least 2 mice per time point) and scRNA-seq (same representation as in Figure 1f-h). As above, the scRNA-seq used n=2 (ZT6, ZT18) or n=3 mice (ZT0, ZT12).

Figure 5. Space-time logic of compartmentalized hepatic functions for Z+R genes.

(a) KEGG analysis of the Z and Z+R genes. For clarity, only KEGG pathways (second column) with adjusted p-value < 0.01 (standard Enrichr output test) are presented (Supplementary Table 3 for all enriched functions). The percentage of central genes is represented by a blue-red gradient. *: appears central due to the KEGG annotation system; however, fatty acid elongation is biased portally (see text). (b-c-d) KEGG analysis of Z+R genes. Representations of genes in central (b-c) and portal (d) enriched Z+R categories (Supplementary Table 3). Polar representation: peak expression times are arranged clockwise (ZT0 on vertical position) and amplitudes (log2, values indicated on the radial axes) increase radially. The radius coordinate of genes with an amplitude >0.9 is halved (indicated with a black circle around the colored dot).

Figure 6. Rhythmic activity of Wnt signaling.

(a) Enrichment/depletion at different times (window size: 3h), of both positive (N=471) and negative (N=149) Wnt targets (background: all R and Z+R genes). Colormap shows p-values (two-tailed hypergeometric test): red (blue) indicates enrichment (depletion). (b) Heatmaps representing scRNA-seq profiles of the top 50 Wnt targets (according to the Apc-KO fold change, Lgr5 was also added) showing rhythmic mRNA in bulk (p<0.01, harmonic regression, F-test, data from ref.). The profiles are computed in three different zones of the central-portal axis: central (layers 1-2), mid-lobular (layers 3-4-5) and portal (layers 6-7-8). Gene profiles (log2) are mean-centered in each zone. An enrichment of the phases around ZT8-14 can be observed, in agreement with Figure 6a. (c) Nuclear protein abundance from ref. of the Wnt effector TCF4 (encoded by the Tcf7l2 gene) in mouse liver shows a rhythm (p=0.003, harmonic regression, n=2, sampled every 3h) peaking at ZT7.5, consistent with the accumulation of mRNA targets a few hours later (panel a). (d) mRNA profiles from bulk RNA-seq (top) and scRNA-seq (bottom). Top five targets with the highest Apc-KO fold change, and rhythmicity in the bulk data (p<0.01, harmonic regression, F-test, data from ref.). Rhythmicity (indicated above each panel) is also computed in three different zones for the scRNA-seq data. As above, the scRNA-seq used n=2 (ZT6, ZT18) or n=3 mice (ZT0, ZT12).

Figure 7. Wnt targets could be explained by rhythmically expressed Wnt ligands from non-parenchymal cells.

(a) Representative smFISH images of Wnt2, Dkk3, and Rspo3 expression at ZT0 (left) and ZT18 (right), shown in green. Markers of non-parenchymal cells (NPCs) are shown in red (Methods). Nuclei are stained in blue (DAPI). Scale bars, 2 μm. (b) Violin plots representing quantitative analysis of smFISH images (n=1420 cells from 189 central veins of at least two mice per time point). Wnt2, Dkk3, and Rspo3 transcripts are quantified in NPCs lining the central vein (Methods). mRNA expression is in smFISH dots per μm³. P-values are obtained from Kruskal-Wallis test.