Temporal analyses of postnatal liver development and maturation by single-cell transcriptomics

Abstract

The postnatal development and maturation of the liver, the major metabolic organ, are inadequately understood. We have analyzed 52,834 single-cell transcriptomes and identified 31 cell types or states in mouse livers at postnatal days 1, 3, 7, 21, and 56. We observe unexpectedly high levels of hepatocyte heterogeneity in the developing liver and the progressive construction of the zonated metabolic functions from pericentral to periportal hepatocytes, which is orchestrated with the development of sinusoid endothelial, stellate, and Kupffer cells. Trajectory and gene regulatory analyses capture 36 transcription factors, including a circadian regulator, Bhlhe40, in programming liver development. Remarkably, we identified a special group of macrophages enriched at day 7 with a hybrid phenotype of macrophages and endothelial cells, which may regulate sinusoidal construction and Treg-cell function. This study provides a comprehensive atlas that covers all hepatic cell types and is instrumental for further dissection of liver development, metabolism, and disease.

Keywords: construction of metabolic zones; endothelial cells; functional maturation of liver; hepatocyte heterogeneity; interaction of macrophages; postnatal liver development; single cell transcriptomics.

Figure 1.. scRNA-seq identifies hepatic cell types in developing and adult liver

(A). UMAP visualization of liver cell types at D1, D3, D7, D21 and D56. Colors indicate cell types, including hepatocyte (Hep-neonatal from D1, D3 and D7; Hep-D21; Hep-D56), endothelial cell (EC), hepatic stellate cell (HSC), cholangiocyte, fibroblast, mesothelial cell (meso), megakaryocyte, erythroid cell (pro-erythroblast, erythroblast, erythrocyte), T cell, natural killer (NK) cell, B cell (Pro-B, large pre-B, small pre-B, B), dendritic cell (classical dendritic cell 1 – cDC1, classical dendritic cell 2 - cDC2, plasmacytoid dendritic cell - pDC, activating dendritic cell - aDC), monocyte, Dcn⁺ macrophage (Dcn⁺ Mac), Kupffer cell, neutrophils (immature neutrophil – iNP, intermediate mature neutrophil – imNP, mature neutrophil - mNP), basophil, granulocyte-monocyte progenitor (GMP) and hematopoietic progenitor cell (HPC). (B). Expression of selected markers for cell types. The dot size corresponds to the ratio of cells expressing the gene in the cell type. The color scales correspond to the averaged expression levels. See also Table S2. (C). Temporal UMAP visualization of hepatic cells from Figure 1A. The UMAP was separated and colored by time points (n = 2~4 for each time point). (D). tSNE map of hepatocytes from 5 time points (indicated by 5 colors). Cells in Hep-neonatal, Hep-D21 and Hep-D56 from Fig.1A were segregated and re-analyzed. (E). Ridge plots displaying the expression distributions of indicated markers in hepatocytes from Fig.1D at five time points. X-axis indicates log-normalized expression levels. (F). tSNE map displaying Mki67 expression in hepatocytes from Fig.1D. (G). Violin plots displaying expression levels of indicated markers in hepatocytes from Fig. 1D. Y-axis indicates log-normalized expression levels.

Figure 2.. scRNA-seq identifies distinct transcriptome profiles in hepatocytes at each time point

(A-E). tSNE map of hepatocytes from D1 (A), D3 (B), D7 (C), D21 (D) and D56 (E). Colors indicate subpopulations identified. (F). tSNE map displaying Scd2 expression in D1 hepatocytes from Fig.2A. (G-H). Violin plots displaying Scd1 (G) and Scd2 (H) expression in hepatocytes. Y-axis indicates log-normalized expression. (I-J). qRT-PCR analysis of Scd1 (I) and Scd2 (J) expression in isolated hepatocytes. Statistical analysis used student’s T-test. Values are presented as means ± SD. (*** p < 0.001, ** p < 0.01, * p < 0.05). See also Figure S1L for cell type marker expression. (K). Scatter plot displaying correlation between Malat1 expression (x-axis) and mtRNA percentage (y-axis). Pearson correlation coefficient (0.24) is shown above the plot. Malat1 expression is log normalized value. mtRNA percentage is the proportion of transcripts mapped to mitochondrial genes. Colors indicate D21 subpopulations in Fig.2D. (L). (Left) smFISH of Malat1 and Neat1 at D7 and D21. (Right) The quantitative results showing fluorescence signal intensity (MFI) and correlation coefficient between Malat1 and Neat1. Scale bar, 100 μm. (M). Density plot displaying distribution of predicted layers of hepatocytes. The predicted layer for each hepatocyte was normalized from layer 1 (pericentral) to layer 9 (periportal). See also Table S3 for gene signatures used for zonation prediction. (N). (Left) Immunostaining of E-cadherin and CYP2E1 in D3, D7 and D21 livers. Scale bar, 100 μm. The insets show loss of E-Cad around CV over time, establishing a clear zonation. (Right) Quantitative results showing the changing expression levels from PV to CV as the arrow indicates. (O). tSNE map displaying predicted layers for D56 hepatocytes from 2E. Color bar indicates predicted layers from Layer 1 (CV area) to Layer 9 (PV area).

Figure 3.. Dynamic changes of transcription factor activities and metabolic functions in hepatocytes

(A). Pseudotime analysis of hepatocyte development from D1 to D56 with Monocle 2. (B). Pseudotime analysis of data in 3A. Color indicates inferred pseudotime used for 3C-3K. (C). Density plot displaying distribution of inferred pseudotime (x-axis). (D). Heatmap representing trends of differentially expressed genes as a function of inferred pseudotime. Genes in row are grouped into 6 clusters based on expression patterns. See also Figure S10 and Table S4 for the gene list and enriched pathways. (E). Scatter and fitted plots of Bhlhe40 scaled expression and activity values along pseudotime. The inferred pseudotime was stretched from 0 to 100. (F-I). Scatter plot and fitted plots of enrichment scores for indicated pathways along inferred pseudotime. Only fitted plots are shown for 3H and 3K. (J). Heatmap displaying genes from the mevalonate pathway with differential expression along inferred pseudotime. (K). Fitted plots of enrichment scores for indicated pathway along inferred pseudotime.

Figure 4.. Development of liver endothelial and mesenchymal cells

(A). tSNE map of endothelial cells from D1 to D56. Cells in EC from figure 1A were segregated and re-analyzed. (B). Time point compositions of each EC subpopulation labelled in 4A. (C-D). UMAP visualization of cells from D1, D3, D7 and embryonic cells (including hemangioblasts, hematopoietic cells, HSEC and endothelium) from published data (Lotto et al., 2020). Colors indicate endothelial cell subpopulation labelled in figure 4A (C), or embryonic cell types (D). (E). Pseudotime analysis of cells from 4C with Monocle 3. Root cells are labelled in white circle 1; branch points are labelled in black circles with numbers. (F). Pseudotime analysis of ECs from D1 to D56 with Monocle 2. Color indicates inferred pseudotime. (G). Fitted plot of scaled expression and activity values along pseudotime of indicated genes. The inferred pseudotime from 4F was stretched from 0 to 100. (H). tSNE map of mesenchymal cells from D1 to D56. Cells in meso, HSC and fibroblast from figure 1A were segregated and re-analyzed. Colors indicate cell types assigned. (I). tSNE map of mesenchymal cells from 4H. (J). Time point compositions of each cell type labelled in 4H.

Figure 5.. Dynamic changes of hematopoietic and immune cell populations

(A). Force-directed graph (FDG) of HPCs, GMPs, erythroid cells, neutrophils, B cells, DCs, basophils, monocytes and Kupffer cells. (B). Cell type compositions of immune cells at each time point. (C). tSNE map of T and NK cells from D1 to D56. Cells in T and NK clusters from figure 1A were segregated and re-analyzed. (D). tSNE map displaying Foxp3 expression in cells from 5C. (E). Percentages of Treg cells out of all T cells at indicated time point. (F). Violin plots displaying differentially expressed genes in Treg cells at D7, D21 and D56.

Figure 6.. A subtype of macrophages emerges transiently around postnatal day 7

(A-B). UMAP visualization of Adgre1 (A) and Csf1r (B) expression levels in all cells included (Figure 1A). (C). tSNE map of Kupffer cells from D1 to D56. Cells in Kupffer clusters from figure 1A were segregated and re-analyzed. (D). Percentages of indicated cell types out of total immune cell population. (E). Dot plots displaying expression levels of selected markers for indicated cell types. Dot size corresponds to the ratio of cells expressing the gene in the cell type. The color corresponds to the scaled average expression level. (F). UMAP visualization of Lyve1 expression levels in all cells included in this study (Figure 1A). (G). Immunostaining of F4/80, ERG and LYVE-1. Representative image taken under confocal microscope. Dcn⁺ Mac cells were indicated by white arrowhead. Scale bar, 10 μm. (H). Quantified FACS data of CD146⁺ F4/80⁺ cells in isolated NPCs, showing an enrichment at D7. Statistical analysis was done with Kruskal-Wallis test, followed by Dunn’s test. (** p < 0.01). (I). RNA velocities of Kupffer cell, EC and Dcn ⁺ Mac from D1 to D56, visualized on UMAP. (J). Clustering analysis of cells from 6I. The Kupffer cell cluster showing differentiation potential to Dcn ⁺ Mac was labelled by asterisk (transitioning KC). (K). Time point composition of asterisk labelled cluster (transitioning KC) from 6J.

Figure 7.. Predicted hepatic cell-cell interactions

(A). Bar plot comparing numbers of significant ligand-receptor interactions between indicated cell types and Dcn⁺ Mac or Kupffer cells. (B). Dot plot of selected ligand-receptor interactions between Dcn⁺ Mac and indicated cell types (Dcn⁺ Mac|LSEC, Dcn⁺ Mac|Treg, Dcn⁺ Mac|HSC) or between Kupffer and LSEC (Kupffer|LSEC, Kupffer|Treg, Kupffer|HSC). For example, ‘moleculeA_moleculeB in cellC|cellD’ indicates the putative interaction between molecule A expressed by cell type C and molecule B expressed by cell type D. P values are indicated by circle sizes, and the means of the average expression levels of interacting ligand and receptor are indicated by color. (C-D). Dot plot of selected ligand-receptor interactions between hepatocytes and LSECs (C) or Kupffer cells (D). Interacting pairs in red indicate ligands expressed by hepatocytes; interacting pairs in blue indicate ligands expressed by LSECs or Kupffer cells.