A single-cell atlas of mouse lung development

Abstract

Lung organogenesis requires precise timing and coordination to effect spatial organization and function of the parenchymal cells. To provide a systematic broad-based view of the mechanisms governing the dynamic alterations in parenchymal cells over crucial periods of development, we performed a single-cell RNA-sequencing time-series yielding 102,571 epithelial, endothelial and mesenchymal cells across nine time points from embryonic day 12 to postnatal day 14 in mice. Combining computational fate-likelihood prediction with RNA in situ hybridization and immunofluorescence, we explore lineage relationships during the saccular to alveolar stage transition. The utility of this publicly searchable atlas resource (www.sucrelab.org/lungcells) is exemplified by discoveries of the complexity of type 1 pneumocyte function and characterization of mesenchymal Wnt expression patterns during the saccular and alveolar stages - wherein major expansion of the gas-exchange surface occurs. We provide an integrated view of cellular dynamics in epithelial, endothelial and mesenchymal cell populations during lung organogenesis.

Keywords: Lung development; Mouse; Progenitor cells; RNA velocity; Single-cell transcriptomics; Type 1 pneumocyte.

Fig. 1.

Capturing the process of alveolarization at a single-cell resolution. (A) Schematic of time points sampled, overlaid onto histologically identified stages of murine lung development. (B) Experimental workflow: at least four mice were pooled per time point; lungs were harvested and single-cell suspensions were generated by enzymatic digest; viable, CD45− Ter119− cells were sorted by FACS for scRNA-seq library preparation. (C) UMAP embedding of all 102,571 cells that clustered, annotated by cell type. (D) Hallmark identifier genes for each cluster are plotted where the size of the dot indicates the proportion of cells within a cluster expressing that gene, and level of blue saturation indicates relative expression level.

Fig. 2.

Cell populations in the lung change markedly over time. (A-C) Relative proportion of cell types as a function of time, as detected by scRNA-seq, in epithelium (A), endothelium (B) and mesenchyme (C). Error bars indicate standard error of the mean from at least four mice from each time point after E12. (D-F) Genes in epithelial cells (D) with statistically significant expression changes as a function of time (q<0.05) by Monocle3 analysis are plotted in a heatmap. Gene expression was log transformed and normalized from 0 to 1 and arranged using hierarchical clustering before plotting. Saturation shading of yellow to blue: high to low expression-level spectrum. Expression patterns were separated into four groups (k=4). When hallmark genes (Fig. 1D) were included in the heatmap, their location is indicated, with matching color to stacked area plots and UMAP above. Genes without statistically significant expression changes over time are not indicated. Similar analysis was performed for the endothelial cells (E) and the mesenchymal cells (F).

Fig. 3.

Distinct epithelial cell identities are established by E18. (A) UMAP embedding of lung epithelial cells (n=10,918) colored by cell type. Although epithelial cells from early lungs (E12 and E15) do occupy distinct UMAP space (‘early epithelium’), their broad spread is consistent with substantial heterogeneity and priming towards several nascent lineages. (B) RNA-velocity vectors were calculated using CellRank and overlaid on the UMAP embedding, with line thickness indicating velocity magnitude. (C) UMAP embedding colored by time point. (D) Epithelial cells from prenatal time points (E12, E15, E16 and E18) are projected on a UMAP and colored by cell type. The transitional epithelial cell cluster is circled in a dotted line. Marker gene expression by cluster in the prenatal epithelium is displayed in a dot plot (below), in which higher expression is represented as a darker color. The size of the dot indicates the proportion of cells expressing that marker. (E) Epithelial cells from postnatal time points (P0, P3, P5, P7 and P14) are projected on a UMAP and colored by cell type. The transitional epithelial cell cluster is encircled in a dotted line. Marker gene expression by cluster in the postnatal epithelium is displayed as a dot plot (below). (F) A row-normalized Jaccard index was calculated between clusters identified in the current analysis and transitional epithelial clusters identified in other studies. The intensity of color indicates relative similarity, with darker colors indicating higher similarity. (G) Cell-trajectory inference was calculated with CellRank, and probabilities of becoming a ‘transitional’ cell plotted. (H) Terminal state probabilities of the transitional cells. The bars are colored based on fraction of cells in each cluster present prenatally (gray) or postnatally (white). Central bars represent the median and range bars represent the interquartile range. (I) RNA ISH of Mdk (red) and Epcam (white, epithelial marker) at selected time points. For additional images see Fig. S2H. (J) Quantification of ISH using HALO for percent of Mdk+ epithelial cells (Epcam+Mdk+/Epcam+) cells over total epithelial cells. A non-parametric Kruskal–Wallis test was performed (***P<0.001). (K) RNA ISH of Cdkn1a (red, transitional cell marker), Sftpc (green, type II cell marker) and Hopx (white, type I cell marker) at selected time points. Cells that express Cdkn1a, Sftpc and Hopx are outlined in a dotted line. For channel-separated images see Fig. S3. (L) Quantification of ISH for percentage of transitional cells (Sftpc+Hopx+Cdkn1a+) cells over total epithelial cells. (M) Confocal micrographs of immunofluorescence and RNA ISH of pan-cytokeratin (PanCK, red), Sftpc (green) and Cdkn1a (white) at selected time points. Cells that express Cdkn1a and PanCK are outlined in a dotted line. For additional time points see Fig. S2I. Insets show magnification of boxed areas. (N) Distance from each transitional cell (PanCK+Cdkn1a+) to the nearest AT2 cell (Sftpc+) was calculated and plotted for each cell. Box and whisker plots show median values (middle bars) and first to third quartiles (boxes); whiskers indicate 1.5× the interquartile range, and dots are points outside the interquartile range. Scale bars: 25 µm.

Fig. 4.

Endothelial cells undergo specification by E15, with subsequent emergence of Car4+ aCap cells by E18. (A) UMAP embedding of lung endothelial cells (n=22,237) colored by cell type. (B) RNA-velocity vectors of endothelial cells were calculated with CellRank and overlaid on the UMAP embedding. (C) UMAP embedding colored by time point. (D) UMAP embedding of endothelial cells colored by Kdr expression: darker colors indicate greater expression. (E) UMAP of Car4 expression. (F) Latent time analysis: latent time increases along the x-axis, shows expression of Kdr in aCap cells. (G) RNA ISH of Pecam1 (red; endothelial cell marker) and Car4 (white; aCap cell marker). A complete series is presented in Fig. S4F. (H) Quantification of ISH was performed using HALO; fraction of Pecam+ cells that are Car4+ was plotted. Boxplots represent pooled data, with each dot representing the percentage of Pecam+Car4+/Pecam+ cells per image (n=20 images per time point, from n=3 individual mice). Increases were observed from E15 to P0 and from P0 to P14. (I) RNA ISH of Kdr (red), and Car4 (white). Additional images are presented in Fig. S5. Insets show magnification of boxed areas. (J) Quantification of Kdr expression in Car4− (black) and Car4+ (red) cells. *P<0.05, ***P<0.001, ****P<0.0001: Mann–Whitney U-test. Box and whisker plots show median values (middle bars) and first to third quartiles (boxes); whiskers indicate 1.5× the interquartile range, and dots are points outside the interquartile range. Scale bars: 25 µm.

Fig. 5.

The mesenchymal cell population undergoes dynamic population shifts between P0 and P3. (A) UMAP embedding of lung mesenchymal cells (n=67,035) colored by cell type. (B) RNA-velocity vectors overlaid on the UMAP depict inferred differentiation/maturation trajectories. (C) UMAP embedding colored by time point. (D) RNA ISH of Tgfbi (red; myofibroblast marker) and Wnt2 (white; Wnt2+ fibroblast marker) at E15, E18, P0 and P3. Images are stitched scans. A complete time series of representative images is presented in Fig. S7C. (E) HALO analysis of the Tgfbi and Wnt2 ISH and assignment of cell identities based on marker expression. A Clark-Evans clustering index was calculated from individual cells and aggregated by image. Values >1.0 indicate that cells are more regularly spaced apart (ordered), whereas values <1.0 indicate progressively increased clustering. Box and whisker plots show median values (middle bars) and first to third quartiles (boxes); whiskers indicate 1.5× the interquartile range, and dots are points outside the interquartile range. (F) A row-normalized Jaccard index was calculated between clusters identified in this study and those from other studies. Increased saturation of blue box shading indicates higher similarity. (G) RNA ISH of Tgfbi (red; myofibroblast marker), Sftpc (green; type II cell marker) and Wnt5a (white) at E15, E18, P0 and P3. Images are stitched scans. A complete time series of representative images is presented in Fig. S8A. Insets show magnification of boxed areas. Scale bars: 100 µm.

Fig. 6.

Epithelial alveolar type I (AT1) cells express components of elastin and the basement membrane. (A) Basement membrane gene expression in AT1 cells over time was plotted, with dot size indicating proportion of cells within a cluster expressing that gene and color saturation indicating relative expression level. (B) Latent time analysis shows expression of Lama3 in AT1 cells. (C) Expression of elastin components in AT1 cells over time is plotted. (D) Latent time analysis shows expression of Fbln5 in AT1 cells. (E) Basement membrane gene expression was plotted in alveolar epithelial cells and other cell types grouped into ‘Other Epithelium’, ‘Endothelium’, and ‘Mesenchyme’. (F) Elastin component gene expression in indicated cell types was plotted. (G) RNA ISH of Lama3 (red) and Hopx (white, AT1 marker). A complete series is presented in Fig. S9A. (H) ISH quantification using HALO. Box plots represent pooled data, each dot representing the percentage of Lama3+Hopx+/Hopx+ cells per image. (I) RNA ISH of Col4a3 (red) and Hopx (white). For additional images see Fig. S10B. (J) ISH quantification representing percentage of Col4a3+Hopx+/Hopx+ cells per image. (K) RNA ISH of Fbln5 (red) and Hopx (white). For additional images see Fig. S10A. (L) ISH quantification representing percentage of Fbln5+Hopx+/Hopx+ cells per image (***P<0.001 by Mann–Whitney U-test). Box and whisker plots show median values (middle bars) and first to third quartiles (boxes); whiskers indicate 1.5× the interquartile range, and dots are points outside the interquartile range. Scale bars: 25 µm.

Fig. 7.

Functional development of the lung is a result of spatial and temporal organization of multiple cell types. An abbreviated summary model of the developing lung compared between two representative time points: mid-embryogenesis (E15) and postnatal (P3). At E15, the lung epithelium is composed of indistinct, uncommitted cell types (as defined by marker-gene analysis) that are very rare or absent in adulthood. The prenatal endothelium in the lung parenchyma is composed of general capillary (gCap) endothelial cells, whereas the majority of cells in the remaining tissue are Wnt2+ fibroblasts, with interspersed clusters of myofibroblasts. By P3, many cell types observed in mature lungs are already present, including AT1, AT2 and transitional epithelial cells. The gCap endothelium has differentiated into Car4+ alveolar capillary (aCap) cells, alongside the gCap endothelium, and mesenchymal cells have undergone extensive spatial reorganization. Representative cell clusters are highlighted in the inset UMAPs derived from Fig. 1C.