Early human fetal lung atlas reveals the temporal dynamics of epithelial cell plasticity – Nature Communications

Cell composition of the human fetal lungs

A total of 19 freshly isolated fetal lung tissues from elective pregnancy terminations ranging from gestational weeks (GW) 10–19 were collected and processed for 3’ single-cell RNA sequencing (scRNA-seq, 10X Genomics). One tissue was sampled twice from different regions of the lung to ensure even capture and representation of the lungs (denoted with an asterisk in Supplementary Table 1). A total of 170,256 single cells were sequenced at ~ 60,000 read depth per cell. FASTQ files were generated, and the reads were aligned with the human reference genome hg38 (GRCh38) (Supplementary Fig. 1A). After the removal of low-quality cells, dead cells, and doublets, 156,698 cells were analyzed. After normalization and dimensionality reduction, reciprocal principal component analysis (rPCA) was used to identify anchors for integration to generate the fetal lung dataset. Twelve male and seven female lungs were collected based on the expression of the male sex-determining genes, SRY and DDX3Y, as well as the XIST (Fig. 1A). Leiden clustering informed by Clustree²³ was used to unbiasedly define the number of cell clusters (Supplementary Fig. 1B). Differentially expressed genes (DEG) based on non-parametric Wilcoxon rank sum test were used to annotate the cell populations which were broadly classified into 5 uniquely different cell populations (Fig. 1B and Supplementary Data 1) designated as stromal (N = 98,166; COL1A1, COL1A2), epithelial (N = 16,068; EPCAM), endothelial (N = 13,376; PECAM1), immune (N = 16,258; CD74, CD3, NKG7), and schwann cells (N = 512; NRXN1, S100B, PLP1) (Fig. 1C). Further sub-clustering of each main cell population was informed by Clustree and DEG which resulted in a total of 58 distinct cell types/states in the developing lungs (Fig. 1D). Assessing the cell proportions over gestational time showed stromal cells constituted the largest proportion of cells throughout the 10 weeks of development (Fig. 1E).

A Biological sex of the lungs was determined by SRY, XIST, and DDX3Y expression. The graphical image was created with Biorender.com. B UMAP visualization of the dataset highlighting the main cell types: stromal, epithelial, endothelial, immune, and Schwann cells. C Gene expression heatmap of the top 5 differentially expressed genes in each cell type. D UMAP projection of all 58 cell types/states identified within the integrated dataset. E Proportion of the major cell types across all gestational weeks. F Integrated UMAP projection of three published fetal lung datasets of similar developmental time points (He et al., Sountoulidis et al., and Cao et al., in gray) overlayed with our dataset (red). Source data are provided as a Source Data file.

Full size image

We performed an integrated analysis to identify similarities and differences between our dataset and three recently published datasets encompassing similar developmental lung tissue^10,11,12. An integrated UMAP overlay of our dataset (in red, Fig. 1F) with the published datasets by Cao et al., He et al., and Sountoulidis et al., (combined in gray) identified 54 of the 58 cell types (~ 93%) were common amongst all the datasets (Supplementary Table 2 and Supplementary Fig. 1C). However, 4 unique cell types were exclusively captured in our dataset and include basal cells, mature ciliated cells, and submucosal gland (SMG) basal cells in the epithelial cluster and monocyte 2 from the immune cluster. This may reflect the differences in cluster annotations, cluster resolution, and differences in DEG associated with each cell type in the datasets and/or sampling error. Additional comparisons between our dataset with each published fetal atlases are detailed in Supplementary Fig. 2 and Supplementary Note. 1. Future efforts to integrate all the datasets and harmonize cluster annotations will strengthen the collectively rich resource of fetal lung atlases available.

Herein, we will describe the epithelial and stromal cell compartments and additional analyses of the fetal lung endothelium in Supplementary Fig. 3 and Supplementary Note. 2, and immune cells in Supplementary Fig. 4 and Supplementary Note. 3.

We next compared our fetal lung dataset to a published human adult lung single-cell atlas using the MapQuery function in Seurat to identify cells of similar biological states of the adult lung⁷ onto our reference fetal lung dataset. As the percentage of anchors to query cells was less than 15% (~ 4%), only high confidence prediction scores (greater than 0.9) were used to predict cell identity. This resulted in 83% of fetal lung cell types captured in the adult lung (Supplementary Fig. 5A). Notable cell types that were not captured included tip cells and NKX2-1 + SOX9 + CFTR +  cells. We subsequently reversed the reference and query datasets and found only 56% of adult lung cell types captured in our fetal lung dataset (Supplementary Fig. 5B). Cell types that were not captured included alveolar-associated cell types (i.e., alveolar fibroblasts and alveolar macrophages), nasal cells (goblet (nasal) and multiciliated (nasal), and ionocytes.

We also compared our dataset to a previous mouse embryonic single lung cell atlas by Negretti et al.⁹ and found 62% common cell types (Supplementary Fig. 5C). It should be noted that the mouse lung dataset by Negretti et al captures tissues from embryonic day 12 and 15 (similar to pseudoglandular lung tissues in human). Therefore, the large number of cells not shared between our fetal lung dataset with the mouse embryonic lung datasets may be explained by the differences in developmental stage of the tissue as well as species differences.

We also used the 10X Xenium transcriptomic platform to resolve the spatial expression of a curated set of genes on archived GW15 and GW18 human fetal lung tissues. We used a custom panel consisting of 289 genes from the 10X genomics human lung gene panel and 50 custom genes informed by our single-cell RNA-sequencing dataset (Supplementary Data 1). Clustree was then used to inform Leiden clustering and cell populations (epithelial 1-2, stromal 1–5, endothelial 1, and immune 1) were determined by DEGs (Supplementary Fig. 1D–I). We identified 9 major clusters in GW15 tissue and 7 major clusters in GW18 tissue (Supplementary Fig. 1F, H). Annotated clusters were then spatially mapped to show the distinct regional distribution of each cluster (Supplementary Fig. 1J, K).

The developing lung stroma

Stromal cells formed the largest proportion of cells in all the lung tissues analyzed (Fig. 1B). We identified 10 stromal cell subtypes (Fig. 2A) based on DEGs (Supplementary Data 1). These included lipofibroblasts, lipofibroblast precursors, cycling fibroblasts, early fibroblasts, airway fibroblast progenitors, uncommitted cells, airway smooth muscle cells (SMC), vascular SMC, chondrocytes, and an unknown stromal cell population. Lipofibroblasts and lipofibroblast precursors were the most abundant stromal cell types across all gestational weeks (Fig. 2B). Both cell types express TCF21, a marker previously shown in both mouse fetal and adult lung lipofibroblasts²⁴. Lipofibroblast precursors differentially expressed higher levels of several canonical mesenchymal associated genes, including PLEKHH2 and MACF1 (Fig. 2C), similar to mouse lipofibroblasts cell lineages²⁵. Regulon activity of the top differential transcription factors identified JUN, FOS, and STAT3 in lipofibroblasts (Fig. 2D), which are known to regulate the synthesis of structural extracellular matrices such as collagen and elastin, and are important in alveolar development²⁶. While alveolar development occurs later in gestation, it is conceivable that lipofibroblasts in early fetal lung development play a role in laying down the ECM proteins required to form airway structures. Differential regulon activity of the transcription factors driving lipofibroblast precursors include NR2F1 and ETV1, which regulates cell growth, proliferation and differentiation^27,28 (Fig. 2D). Airway fibroblasts expressed high levels of TWIST1, NFIX, and ZEB1, a core epithelial-to-mesenchymal (EMT) transcription factor¹² (Fig. 2D).

A UMAP visualization of the fetal stromal cell subtypes. B Proportion of stromal cell subtypes across gestational week. C Gene expression heatmap of DEGs representing each cell subtype. D Heatmap measuring regulon activity of the top differentially expressed transcription factors (TF) genes based on regulon specificity score (RSS) via SCENIC. E Xenium spatial plots of EGR1, ATF3, IRF1 in lipofibroblasts, TCF21 and MACF1 in lipofibroblast precursors, and SERPINF1, FBN1, IGFBP5 in airway fibroblast progenitors. Colors in (B) indicate cell type as in (A). Numbers in the top bar of (C) and (D) indicate cell type as in (A). Source data are provided as a Source Data file.

Full size image

We next used Xenium, high resolution multiplexed in situ spatial profiling, to determine the spatial localization of the expressed DEG in each stromal cell subtype. We identified genes associated with airway fibroblast progenitors (SERPINF1, FBN1, IGFBP5), lipofibroblast precursors (MACF1 and TCF21), and lipofibroblasts (EGR1, ATF3, IRF1) in the stromal regions of the fetal lung tissues (Fig. 2E). Airway fibroblast progenitors were primarily concentrated around the large airways (epithelial 1), while lipofibroblast precursors were seen throughout the stroma. In addition, lipofibroblasts were observed scattered among the lipofibroblast precursors surrounding the epithelium.

Cycling fibroblasts and early fibroblasts were also present at a relatively high proportion in the early lung tissues (GW 10–14; ~ 6–10% and ~ 9–12% respectively) but gradually reduced in later gestational lung tissues (GW15 onwards; ~2–5% and ~3–7% respectively) (Fig. 2B and Supplementary Data 2). Other than the common collagen genes (COL1A1, COL1A2), cycling fibroblasts expressed high levels of genes associated with cellular proliferation TOP2A, CENPF, MKI67 (Supplementary Data 1 and Fig. 2C). Early fibroblasts expressed high levels of genes associated with chromatin remodeling such as HELLS and TYMS²⁹ (Fig. 2D), previously shown to be highly expressed in mouse embryonic fibroblasts. The transcription factors with differentially high regulon activity in early fibroblasts include E2F7, E2F8, FOXM1 (known regulators of the cell cycle in fibroblasts³⁰), TCF7, and E2F1.

There was also a significant proportion of uncommitted cells that were closely associated with cycling fibroblasts and lipofibroblast precursors based on the UMAP projections (Fig. 2A). These uncommitted cells shared some transcriptional similarities to cycling fibroblasts and lipofibroblast precursors but also expressed genes associated with other cell types in general (Fig. 2C). Airway SMC differentially expressed higher levels of TAGLN, DES, and MYH11, while vascular SMC expressed IGFBP7, MEF2C, and EGFL6, both similarly reported in other studies¹¹. The transcription factors and associated regulons of TBX5, MYEF2, and KLF9 were differentially active in airway SMCs, while ZFN41, PITX3, and HMGA1 were found in vascular SMCs. Chondrocytes differentially expressed many of the collagen genes, including COL2A1, COL11A1, COL9A1, and RUNX1, as well as FOXC1 associated regulons known to regulate chondrocyte development^31,32. An unknown small cluster of stromal cells (cluster 10) with distinct DEG expressed high levels of AHSP, ALAS, and BPGN, genes correlated with erythroid cells. These unknown stromal cells also expressed abundant levels of TMEM158, a gene enriched in cells undergoing EMT transition through STAT3 signaling³³. Their role in the fetal lung remains undetermined.

Spatial localization of the top DEGs associated with cycling and early fibroblasts showed these cells were generally found in the stromal regions of lung tissue (Supplementary Fig. 6A). Airway SMCs were found surrounding the developing epithelium while vascular SMC genes concentrated in a distinct region presumed to be the site of a developing vessel (more evident in endothelial cluster of GW18 lung tissue). Colocalization of THBS1 and SOX9 was found in a cartilaginous region of the lung tissue. Immunofluorescence staining of MEF2C, DES, and ACTA2 was used to mark vascular SMC, airway SMC, and pan-SMC, respectively, and was localized to emerging blood vessels and surrounding the epithelium respectively (Supplementary Fig. 6B).

To infer developmental lineages of the fibroblast clusters, we used LatentVelo to estimate RNA velocities on a subset of the stromal cells without airway SMC, vascular SMC, and chondrocytes. A new UMAP was generated for this subset, and we visualized these velocities on the UMAP embedding and with PAGA (Fig. 3A, B). Inferred trajectories suggest that cycling fibroblasts form early fibroblasts that then contribute to airway fibroblast progenitors and lipofibroblasts. We further validated these trajectories with Slingshot, using cycling fibroblasts as the root cells or cell origins, and found a gradual change in gene expression as the cells moved along a differentiation trajectory to lipofibroblasts or airway fibroblast progenitors (Fig. 3C, D).

A LatentVelo analysis revealing inferred cell state trajectories (arrows) projected onto the UMAP. B Partition-based graph abstraction (PAGA) plot identified lineage trajectories informed by LatentVelo velocities. Line thickness indicates inferred transition strength. C Slingshot trajectory analysis with a root at cycling fibroblasts identifies a trajectory to lipofibroblasts (UMAP colored by pseudotime, dark blue to yellow). The trajectory heatmap plot shows the progression of cell types along the trajectory and the change in significantly varying genes. D Slingshot trajectory analysis with a root at cycling fibroblasts identifies a trajectory to airway fibroblast progenitors (UMAP colored by pseudotime, dark blue to yellow). Numbers and colors in (B) indicate cell type as in (A). Source data are provided as a Source Data file.

Full size image

A dynamically evolving fetal lung epithelium

In the fetal lung epithelium, analysis of the top DEG coupled with high expression of canonical epithelial genes associated with specific cells revealed 19 subtypes (Fig. 4A, B and Supplementary Data 1). Many of the epithelial cell types were found in all lung tissues except for ciliated cell precursors, which emerged around GW13 onwards, mature ciliated cells, and SCGB3A2 + FOXJ1 + cells which emerged around GW14. Based on our sampling of tissues for sequencing (see “Methods”), basal cells and SMG secretory cells were not detected until around GW15, and only represented a very small fraction (~0.2–1% and ~1–5%) of the total epithelial cell population (Supplementary Data 2 and Supplementary Fig. 7A). PNEC and SOX2^highCFTR + cells were found in higher abundance in early fetal lungs (up to GW13) which then decreased in proportion in later gestational tissues.

A UMAP visualization of the fetal epithelial subtypes. B Gene expression heatmap of DEGs representing each cell subtype. C Dotplot of average scaled gene expression for CFTR, SOX9, and SOX2 in each epithelial cell type. D Xenium spatial plots of CFTR, NKX2-1, SOX9, SOX2, SCGB3A2, SFTPB in the developing airways of GW15 (top row) and GW18 (bottom row) fetal lung tissues. E Immunofluorescence staining shows the localization of NKX2-1, SOX2, SOX9, CFTR SCGB3A2, and SFTPB positive cells in GW15 (top row) and GW18 (bottom row) fetal lung tissues. At least 4 representative images were captured and analyzed for each gestational time point. The Right column panels are negative controls (no primary, secondary antibodies only). DAPI marks all nuclei. White arrowheads and orange arrowheads demarcate areas of SOX2low and SOX2high cells, respectively. Scale bar = 10 microns. Numbers in the top bar of (B) indicate cell type as in (A). Source data are provided as a Source Data file.

Full size image

Leveraging the expression of SOX2 and SOX9 to determine proximal and distal airway cells respectively³⁴, we identified subclusters of epithelial cells that expressed high levels of SOX2, which included mature ciliated cells, a CFTR-expressing cell type that co-expressed SCGB3A2 and SFTPB (SCGB3A2 + SFTPB + CFTR + ³²), SCGB3A2 + FOXJ1 + and SOX2^highCFTR + progenitor cell populations (Fig. 4C). On the contrary, SOX9 expressing cells included budtip progenitors, tip cells, NKX2-1 + SOX9 + CFTR + progenitor cells, and stromal-like 1 cell populations. Cells that expressed relatively lower levels of both SOX2 and SOX9 included SMG secretory cells, PTEN + STAT3 +, and NRGN + cells. Expression of SOX2 in the budtip progenitors in our dataset was relatively low compared to SOX9, which is consistent with previous findings^35,34,36.

Proliferating progenitors were identified by the elevated expression of genes associated with cellular proliferation MKI67 and cell cycle genes encoding the microtubule-binding protein CENPF, NUSAP1, and CDK1 (Fig. 4B). These cells also expressed abundant levels of BRCA1, MYBL2 and FOXM1 (Supplementary Fig. 7B). Assessment of the top DEG (Supplementary Data 1) associated with these cells suggested these proliferating progenitors were relatively uncommitted.

An interesting finding in our analysis was the identification of four epithelial cell types that expressed relatively higher levels of CFTR (Fig. 4C, average scaled expression greater than 1). These included NKX2-1 + SOX9 + CFTR +, SOX2^highCFTR +, SOX2^lowCFTR +, and SCGB3A2 + SFTPB + CFTR + cell populations. Previous studies in porcine fetal lungs suggest the role of CFTR in branching morphogenesis³⁷. However, the pleiotropic roles of CFTR in the developing lung remain to be determined. The SOX2^lowCFTR+ shared many common DEG with SOX2^highCFTR + cells, but SOX2^lowCFTR + cells also expressed GPRC5A, a gene involved in retinoic acid signaling instrumental to lung development³⁸ (Fig. 4B and Supplementary Data 1), and exhibited high regulon activity in transcription factors NFKB1 and IRF1 (Supplementary Fig. 7B). On the contrary, SOX2^highCFTR + cells contained elevated regulon activity in transcription factors TFDP2 and FOXA2, which has been shown to mediate branching morphogenesis and differentiation³⁹. Interestingly, these cells also expressed high levels of CD47 (Supplementary Data 1), a cell surface molecule used to enrich multipotent human induced pluripotent stem cells (hiPSC)-derived lung progenitors⁴⁰.

NKX2-1 + SOX9 + CFTR + cells expressed high levels of CPM, CLDN18, FGFR2, and MUC1. Interestingly, CPM is also a cell surface molecule used to enrich fetal lung epithelial progenitors generated from hiPSC that can generate alveolospheres⁴¹. The association between these cells would be an interesting future study. NKX2-1 + SOX9 + CFTR + cells expressed high levels of the transcription factor CEBPA, required for lung development and activation of the alveolar development program⁴², THRA, which activates alveolar type II cell proliferation during regeneration⁴³, and SOX5, required for branching morphogenesis in mouse lungs⁴⁴.

Spatial Xenium profiling showed NKX2-1 + SOX9 + CFTR + were restricted to the epithelium (Epithelial 2, blue shade) with enriched expression of all three genes in regions of bifurcations or developing bud areas of GW 15 and 18 lung tissues (Fig. 4D). In contrast, SOX2 + CFTR + cells were found in the stalk regions of the airways. SOX2 expression is known to regulate the cellular proliferation of epithelial cells⁴⁵, and therefore, SOX2^highCFTR + cells may be involved in the expansion of cells to form the developing airways.

Using immunofluorescence staining, we confirmed spatial expression of NKX2-1, SOX9, and CFTR triple-positive cells in the developing lung bud region (Fig. 4E). On the contrary, NKX2-1 + SOX2 + CFTR + cells were found in the fetal airways with distinct clusters of SOX2^high (orange arrowheads) and SOX2^low (white arrowheads) cells. We found no statistical significance when comparing the proportions of SOX2^high and SOX2^low cells at different time points (Supplementary Fig. 7C).

The SCGB3A2 + SFTPB + CFTR + cells referred to as triple positive (TP) herein, differentially expressed higher levels of CXCL17, CP, and CYB5A (Fig. 4B) and differentially high regulon activity of transcription factor HES1, a major NOTCH target gene, and ASCL2, previously shown to be a WNT/CTNNB1 transcriptional target in the gut epithelium⁴⁶ (Supplementary Fig. 7B). These TP cells were found abundantly scattered throughout the airways (Fig. 4D, E).

Unlike in adult lungs^3,47, fetal secretory, ciliated, and SMG cells expressed relatively low levels of CFTR compared to the progenitor cells. The rare ionocytes that express the most abundant CFTR transcripts in postnatal lung tissues⁴⁷ were not detected in any of the fetal lung tissues examined here. However, this does not preclude the development of these cells later in fetal lung development. Fetal club cells expressed high levels of the canonical genes SCGB1A1, SCGB3A2, and CYB5A (Supplementary Data 1) typically found in mature club cells. On the contrary, SCGB3A2 + FOXJ1 + cells expressed abundant genes associated with ciliogenesis (RSPH1, DNAH5, Fig. 4B) suggesting these may be transitory cells originating from a secretory precursor.

Ciliated cell precursors and mature ciliated cells shared the expression of several key genes including DNAH5, a protein-coding gene for microtubule assembly, FOXJ1, a key regulator of cilia gene expression, and RSPH1, a gene encoding a protein that localizes cilia (Clusters 10 and 11, Fig. 4B). Both cell types also expressed high levels of regulon activity from TP73 (Supplementary Fig. 7B), a transcriptional activator of the FOXJ1⁴⁸. Ciliated cell precursors expressed high levels of FOXN4 (Supplementary Data 1), a transcription factor required for the expression of motile cilia genes, and observed in cells undergoing ciliated cell differentiation⁴⁹. Mature ciliated cells, on the other hand, expressed prominent levels of the TUBA1A required for ciliogenesis⁵⁰. Interestingly, we found a unique cell cluster that co-expressed SCGB3A2 and FOXJ1 (SCGB3A2 + FOXJ1 +) which were found dispersed in the airways and found after GW14 (Supplementary Fig. 7D and 7E). The relationship between these cells with ciliated cells remains to be experimentally determined. However, assessment of the proximity scores of the cells from spatial transcriptomics suggests SCGB3A2 + FOXJ1 + cells are closely associated with other proximal cells including SCGB3A2 + SFTPB + CFTR +, SOX2 + cells, and club cells (low-resolution Visium Fig. 5A, high-resolution Xenium Fig. 5B), which further suggests a lineage relationship between these cells.

A Correlation between RCTD weights for Visium spots. Showing the colocalization of cell types within the same Visium spots. Hierarchical clustering shows groups of frequently colocalized cell types. B Xenium-informed RCTD neighborhood enrichment, showing cell types that frequently neighbor. Hierarchical clustering shows groups of frequently enriched neighbors. Source data are provided as a Source Data file.

Full size image

Tip cells and budtip progenitors both expressed genes related to surfactant protein production, such as SFTPC, SFTPA1, ABCA3, ETV5, and TESC (Fig. 4B). Budtip progenitors differentially expressed higher levels of CCN2, ERRFL1, MIR222HG, and genes associated with cellular proliferation (e.g., FOSB). Budtip progenitors were found to contain differentially high regulon activity from transcription factors NFKB1 and SOX9 (Supplementary Fig. 7B), both involved in distal epithelial cell differentiation⁵¹. Transcription factors including ELF5, CEBPA, and USF2 (important in regulating SFTPA expression⁵²), were differentially active in tip cells. While ETV5 was previously shown to regulate human distal lung identity¹² and mouse adult lung AT2 identity⁵³, we identified that its homolog, ETV1⁵⁴, was also abundantly expressed in both human fetal tip and budtip progenitor cells. Notably, we found that our tip cells and late tip cells identified by Lim et al.⁵⁵ shared similar DEGs, suggesting they may be the same or related cells. However, these cells also showed similarity to early and mid-tip cells, and therefore we chose to label them as tip cells and not define their stage explicitly. Co-localization analysis with spatial transcriptomics showed budtip progenitors and tip cells shared common genes (Fig. 4B) and were closely spatially correlated (Fig. 5A), suggesting these cells may also share regionally the same space and may represent the same cell type at different phenotypic states. Indeed, expression of CA2 and ETV5 was highly enriched in budtip and tip cells, while SFTPC was found to be enriched in budtip cells compared to neighboring tip cells (Supplementary Fig. 6C).

In our dataset, basal cells were defined by KRT5, KRT17, TP63, and CXCL14 expression (Fig. 4B and Supplementary Data 1) and were rare in our dataset, especially in younger gestational age lung tissues (often found <10 cells). This may represent insufficient sampling of the larger airways, such as tracheal tissues (which was not collected for sequencing). Immunofluorescence staining for deltaNP63 (basal cell), FOXJ1 (ciliated cell), and SCGB3A2 (secretory cell) markers showed the presence of basal cells in the developing airways to be sparse in early gestational tissue but appeared more prominent in later tissues. Furthermore, SCGB3A2 + FOXJ1 + cells were in GW15 but not in GW12 tissues, supporting the observation in our scRNA-seq analysis that these cells emerged later. The presence of FOXJ1 continues to increase over time (Supplementary Fig. 7D). Transcription factors abundantly active in basal cells included PITX1 and LEF1, both previously shown to regulate basal cell proliferation and differentiation (Supplementary Fig. 7B)^56,57.

Airway SMG basal cells expressed high levels of IGFBP, TIMP3, LGALS1, LGALS3, S100A4, and S100A6, all previously associated with SMG basal cells¹⁰. Similarly, SMG secretory cells expressed high levels of UPK3B, RARRES2, C3, and ALDH1A2 (Supplementary Data 1). Both SMG basal and SMG secretory cells shared elevated activity of the HOX family transcription factors HOXB2, HOXB4, and HOXA5 (Supplementary Fig. 7B), all previously shown to regulate proximal-distal patterning and epithelial differentiation in the developing mouse airways⁵⁸.

Club cells expressed high levels of SGCB3A1, SCGB3A2, SCGB1A1, and CYB5A (Fig. 4B). The TEAD4 transcription factor was differentially active in fetal club cells (Supplementary Fig. 7B), which has previously been shown to regulate club cell development/homeostasis⁵⁹. Club cells spatially correlated with differentiated/differentiating cells including ciliated cells, PNEC, and TP cells (Fig. 5A, B).

Similar to previous studies^10,11, PNEC emerged early in the developing lungs, but the proportion of PNEC appeared to decrease by the end of the pseudoglandular stage ~ GW16 (Supplementary Fig. 7A). PNEC expressed high levels of the canonical genes including CHGB, CALCA, and ASCL1 (Fig. 4B) and the transcription factors NKX2-2, NEUROD1, and ASCL1 were most abundantly expressed in these cells (Supplementary Fig. 7B), confirming their cell identity as previously described^10,11,60. Spatial localization for PNECs markers, ASCL1 and CALCA, showed clusters of these genes in the early GW15 airways but absent in the GW18 airways (Supplementary Fig. 8A). The development of clusters of PNECs, or neuroepithelial bodies, was previously described to play a role in bombesin-mediated lung growth as these PNEC cells were found adjacent to epithelial branching points⁶¹. Similar to previous finding⁹, two PNEC subpopulations separated by differential expression of ASCL1 and CALCA, respectively, were observed in our dataset (Supplementary Fig. 8C).

An unknown epithelial cell cluster was found that expressed high levels of mTOR-related STAT3 and PTEN (PTEN + STAT3 +). No other cell-type-specific gene was differentially expressed in these cells. In the developing mouse lungs, mTOR signaling is required for Sox2 acquisition and conversion of Sox9 + distal progenitors to Sox2 + proximal cells⁶². It is unclear if these PTEN + STAT3 + cells are equivalent to transitional cells of distal to proximal progenitors. However, these cells expressed DEG associated with cell cycle progression such as NFKB1, an NFKB subunit (REL), and IRF1 (Supplementary Fig. 7B). Spatial transcriptomic analysis of the expression of the PTEN, STAT3, and NFKB1 was used for identify these cells using Xenium (Supplementary Fig. 8B). Many stromal cells expressed all three genes but, in the epithelium, all three genes were found colocalized in the same cell and were more prominent in the stalk regions of the developing epithelium.

Other epithelial cell populations that were less defined and made up a smaller proportion of the total epithelium included the NRGN + cells that expressed abundant levels of PPBP (or CXCL17), PF4 (or CXCL4), and ELN, which would suggest these cells may play a role in inflammation and laying down extracellular matrices. Spatial expression of NRGN in GW15 fetal lungs showed a high amount of NRGN molecules in pockets of stromal cells surrounding the developing epithelium with few cells in the epithelium expressing NRGN (Supplementary Fig. 8B). In GW18 lung tissue, NRGN expression was most abundant in the developing endothelium with sporadic NGRN in some epithelial cells. It is unclear if NRGN + cells represent a specific epithelial cell state. Two stromal-like epithelial cell populations were also found in all fetal lung tissues examined. These cells were annotated based on the DEG of stromal associated genes, including COL3A1, which is expressed in squamous cell types⁶³. Stromal-like cells 2 also expressed MEOX2⁶⁴ (Supplementary Data 1) which regulates TGFβ signaling pathway and epithelial-mesenchymal cell transition. These cells also expressed abundant WNT2, previously shown to orchestrate early lung morphogenesis⁶⁵. Future studies will be needed to define the functional role of the stromal-like cells in the developing lung.

Inferred trajectories of epithelial subtypes

Using LatentVelo, we estimated RNA velocity to predict the lineage relationships of the epithelial cell types (Fig. 6A). We also used a separate pseudotime method to arrive at these same results for validation (Supplementary Fig. 9A). We performed slingshot analysis on the late-stage tissues to determine the origins and trajectories of the differentiated cells and found trajectories from budtip progenitors through the TP cells to mature ciliated cells, PNEC, and basal/club cells (Supplementary Fig. 9B). However, subsequent validation with LatentVelo and Palantir pseudotime was unable to resolve the origin of basal cells, likely due to the low numbers of these cells in our dataset.

A LatentVelo analysis revealing inferred cell state trajectories (arrows) projected onto the UMAP. B Terminal states are identified by using velocities with CellRank, shown on the UMAP. C Partition-based graph abstraction (PAGA) analysis reveals multiple connections with the fetal epithelia. Line thickness indicates inferred transition strength. D Analysis of cell type proportion over time (gestational weeks). Each point indicates a tissue sample. The best fit out of a constant, linear, or quadratic model is chosen for each cell type, shown with 95% credible intervals. E LatentVelo PAGA subset into Early (GW10-13), Mid (GW14–16), and Late (GW17-19). The weakening of PNEC transitions were observed, and the emergence of transitions to ciliated cells with increasing GW. F Based on the PAGA, the early trajectories to PNEC from TP cells are analyzed with Slingshot (UMAP colored by pseudotime, dark blue to yellow). Trajectory heatmaps show the progression of cell types (top colored bar) and significantly varying genes along the trajectory. List of GO, KEGG, and Reactome terms and pathways enriched in each gene cluster (row colors). Pathways are shown in text boxes for the cell type that most express the genes. G Mid-time point Slingshot trajectories highlighting the emergence of mature ciliated cells from TP cells are shown. H Based on PAGA, basal cells contribute to late club and PNEC cells. Slingshot trajectories show the development of PNEC cells from basal cells and the development of Club cells from basal and TP cells. Source data are provided as a Source Data file.

Full size image

We then analyzed LatentVelo velocities with CellRank⁶⁶ and identified three terminal states: mature ciliated cells, PNEC, and budtip progenitor cells (Fig. 6B and Supplementary Fig. 9C, D). Using slingshot, we identified a bifurcation in the early budtip progenitor cells towards either the late budtip progenitor cells or tip cells and NKX2-1 + SOX9 + CFTR + cells (Supplementary Fig. 9E).

Based on PAGA using LatentVelo velocities, we observed trajectories that suggest a high degree of cellular plasticity and changes in cell lineage relationships (Fig. 6C). We hypothesized that epithelial lineage development may be temporally regulated. Therefore, we dissected the analysis based on developmental vignettes to unmask these dynamic events. Based on the emergence of several canonical cell types, such as mature ciliated cells and the significant decline of PNEC cells at ~ GW14 (Fig. 6D and Supplementary Fig. 4A), we grouped the tissues into 3 vignettes: Early (GW10-13), Mid (GW14–16), and Late (GW17-19). The late vignette (GW17-19) marks the early canalicular stage of lung development, a period in which differentiation and emergence of many epithelial cell subtypes emerge⁶⁷. A few noteworthy trajectories were observed. In early gestational tissues, there were multiple sources of PNEC contributions including a trajectory from budtip progenitors through tip cells to SOX2^highCFTR to TP and then PNEC. Changes in the DEG reflected the gradual cell fate change to PNEC differentiation (expression of ASCL1 and GHRL) (Fig. 6F). Enriched GO/KEGG/Reactome terms or pathways for the cells along the trajectory reflected the change in phenotype. During the transition from TP cells to PNEC, genes involved in FGF, NOTCH, and TGFβ signaling pathways were upregulated. All the PNEC trajectories weakened in the late stages, and the proportion of PNEC cells at these later stages decreased (Fig. 6D).

Several trajectories converged onto TP in the mid-gestational tissues (Fig. 6E). These included cell sources from SOX2^lowCFTR + cells, SOX2^highCFTR+ cells, stromal-like cell 2 cells, proliferating progenitors, and basal cells. A strong connection (thick arrow line) between NKX2-1 + SOX9 + CFTR + and SOX2^lowCFTR + or SOX2^highCFTR + cells, and the latter two with TP were observed in all time points, suggesting these cells are developmentally related. Interestingly, the connection between budtip progenitors and tip cells appeared to reflect only a subpopulation of budtip progenitor cells that emerged later in development (mid-stage). Cells that expressed abundantly higher levels of SOX9 included budtip and tip cells (Fig. 4D). Our inferred trajectory predictions (Fig. 6F, G) and Xenium spatial localization suggest a relationship between SOX9 + cells. The spatial expression suggests these cells all reside close to one another in the bud or bud-adjacent region tip and next to stalk regions. We, therefore, predict that NKX2-1 + SOX9 + CFTR +, while they do express low levels of SOX9, may represent a transitional state from tip cells as they differentiate or acquire stalk identity. In assessing their DEGs, these cells are distinct from tip and SOX2 + cells but had similarities with both (Supplementary Fig. 9F), further supporting their lineage relationship.

During the mid (GW14–16) time points, the emergence and trajectories for mature ciliated cells were observed, stemming from budtip progenitors through tip cells, NKX2-1 + SOX9 + CFTR + to ciliated precursor cells or through SOX2^highCFTR + cells, TP cells, and SCGB3A2 + FOXJ1 + cells to mature ciliated cells (Fig. 6G). During the transition to SCGB3A2 + FOXJ1 + cells from TP cells, genes involved in interleukin and TGFβ signaling were upregulated. Fetal club cells appeared to contribute to the ciliated cell precursor (Fig. 6E), supporting the notion that fetal club cells may have greater developmental potential as has been observed in mouse lung injury models⁶⁸. However, this developmental trajectory was lost at later time points. Interestingly, there were no connections between fetal basal cells and club cells before GW17, which may be explained by the low numbers of basal cells captured for sequencing. At later GW, trajectories originating from basal cells contributed to PNEC and Club cell lineages (Fig. 6H), suggesting the basal cells at this stage have acquired developmental plasticity similar to fetal TP cells and adult basal cells⁶⁹.

The plasticity of TP cells in PNEC and ciliated cell lineages

Conchola et al.¹³ previously showed that TP cells can functionally contribute to PNEC and ciliated cells in the lower airways through ex vivo lineage tracing studies. In our dataset, we found a temporal contribution of these cells with a decline in PNEC in later tissues and a concomitant increase in ciliated cells around GW14 (Fig. 7A). Our trajectory analysis also suggests TP contributions to basal/club cells in later gestational tissues (>GW16). Future studies will be needed to determine the developmental potential of TP by using similar lineage tracing tools to validate these inferred trajectories¹³.

A Heatmap showing the gene-normalized expression of the top 20 differentially expressed genes for the PNEC and mature ciliated cells within the TP cells versus gestational week. B CellRank matures ciliated and PNEC fate probabilities for TP cells versus gestational week (distributions shown as violin plot). C Immunofluorescence staining for ASCL1 (PNEC marker, yellow) and FOXJ1 (ciliated cell marker, magenta) in GW12, 15, and 18 fetal lung tissues. At least 4 representative images were captured and analyzed for each gestational time point. The Right panel shows the negative controls (no primary or secondary antibodies only). Scale = 25 microns. D FGF signaling interaction probabilities towards TP cells from spatially neighboring cell types with significant FGF interactions. E Co-occurrence probability of deconvolved cell types within spatial neighborhoods of TP cells of a given radius with Xenium. F Significant ligand-receptor (L-R) interactions towards TP cells for the FGF pathway, subset to signaling cell types in proximity. Both color and point size indicate interaction probability. Red values show that L-R interactions significantly increased early. Blue values show interactions significantly increased mid/late. Significance is determined using a permutation test randomly permuting cell-type labels, significant interactions are chosen with p < 0.01. Source data are provided as a Source Data file.

Full size image

CellRank was used to estimate the PNEC and mature ciliated cell fate probability for the TP cells, showing a decline in PNEC probability and increase in mature ciliated cell probability with increasing GW (Fig. 7B). Immunofluorescence staining for PNEC marker ASCL1 (yellow) and ciliated cell marker FOXJ1 (magenta) showed an abundance of PNEC cells in early GW tissue which was lost in GW18 fetal lung (Fig. 7C). Co-occurrence analysis with spatial transcriptomics showed both cell types in proximity to TP cells and SOX2 + cells, suggesting these cells reside in the developing stalk regions and contribute to the lineage composition of the airways (Fig. 5A, B). Together, these results demonstrate a temporal regulation of PNEC and ciliated cell differentiation through TP cells, supporting the dynamic plasticity of these progenitor cells in the developing lung.

Temporal changes in FGF signaling to TP

We revealed several trajectories from mid and late fetal lung tissues which suggests that cell fate acquisition is more dynamic than previously understood⁷⁰. One explanation for these dynamic cell fate changes is that they reflect differences in the spatial microenvironment of the cell and the various signaling mechanisms involved (paracrine, autocrine, distant signaling)⁷¹. To determine the cellular interactions in regulating the change in gene expression and developmental potential of TP cells, we performed CellChat⁷² analysis on our scRNA-seq dataset. Informed by the spatial relationships with spatial transcriptomics, we identified unique signaling pathways targeting the TP cells. Importantly, we observed changes in signaling when assessed through each vignette: early (GW10-13), mid (GW14–16), and late (GW18-19) (Supplementary Fig. 10A, B). Signaling pathways that specifically targeted TP cells relative to other epithelial cells, included fibroblast growth factor (FGF, Fig. 7D), WNT (Supplementary Fig. 11and Supplementary Note. 4), and NOTCH (Supplementary Fig. 12and Supplementary Note. 4) signals. FGF and NOTCH signaling pathway genes were also found to be upregulated during the transition from TP cells to PNEC with trajectory analysis (Fig. 6F).

We then focused our investigation on the FGF signaling pathway (red arrow, Supplementary Fig. 10A) as it is a known pathway involved in lung branching morphogenesis⁷³. Moreover, several FGF recombinant proteins are used in directed differentiation protocols to generate fetal lung epithelia from human iPSC^16,18. Therefore, we reasoned this pathway might play a prominent role in regulating TP cell fate. We specifically focused on cell-cell interactions through FGF signaling with TP cells (Fig. 7D). Informed by Xenium, we identified cells that were spatially in “close proximity” to TP cells which included airway fibroblast progenitors, SOX2^highCFTR + cells, PNEC, club cells, and basal cells in GW15 lungs, as well as lipofibroblast precursors, SCGB3A2 + FOXJ1 +, airway SMC, and mature ciliated cells in GW18 lungs (Fig. 7E). We did not include PNEC, SCGB3A2 + FOXJ1 +, or mature ciliated cells (which are also spatially close to TP cells) as potential source cell types to TP, as we aimed to find the signals that may contribute to the dynamic contribution of TP cells to PNEC/ciliated cells.

We then estimated the interaction probability between specific ligand-receptor (L-R) interactions involved in FGF signaling across GWs between TP cells and FGF signal senders. We specifically looked at airway fibroblast progenitors and airway SMC as they were close in proximity with TP cells based on co-occurring cell types (Fig. 7E). We found significantly increased signaling between L-R interactions involving the ligands FGF2 and FGF7 expressed in airway fibroblast progenitors during in early (GW10-13) gestation, and we found significantly increased signaling for L-R interactions involving the ligand FGF18 expressed in airway SMC during mid (GW14–16) or late gestation (GW17-19) (Fig. 7F). Gene expression of the specific FGF ligands and receptors showed elevated expression of FGFR2 and FGFR3 in TP cells (Supplementary Fig. 10C). We also confirmed the high expression of FGFR2 and FGFR3 using Visium in regions with high proportions of TP cells at GW15 and 18 (Supplementary Fig. 10D). These results suggest greater contribution of FGF2 and FGF7 signaling by airway fibroblast progenitors during early development, and FGF18 signaling by airway SMCs later in development, which may explain the temporal switch from PNEC to ciliated cells contributions.

Similar analyses of the temporal changes in WNT and NOTCH signaling can be found in Supplementary Fig. 11 and 12, respectively, with WNT signaling from airway SMCs, lipofibroblast precursors, airway fibroblast progenitors, and basal cells, and NOTCH signaling from airway SMCs, SOX2^highCFTR + cells, and basal cells.

Overall, future studies will need to confirm the explicit role of these signaling pathways in specific cell lineage development through TP or other progenitor cell types/states and how the cell interprets the wide array of signals it receives to alter cell fate changes.

hPSC-derived fetal lung differentiations capture cellular heterogeneity and trajectories found in the native tissue

Directed differentiation protocols of hPSC aim to capture developmental milestones ensuring robust development of bona fide cell types in cell cultures. To do this, differentiation protocols must reflect developmental processes that are observed in the primary tissues. Many differentiation protocols, including for the lung, capture “fetal-like” states, but few have benchmarked the developmental stage of these cells to fully understand the phenotype of the cells and the potential limitations of the protocols. The human “fetal” stages of tissue development span nearly 9 months of gestation, during which a small batch of early embryonic cells must proliferate, differentiate, and migrate to form a complex tissue and eventual functional organ. As the differentiations become more advanced, with the capability to generate hPSC-derived tissues-specific organoids, these models represent an exciting research tool to study fundamental mechanisms of development. This is especially important when access to primary fetal tissues for research is limited. Here, we sought to determine the developmental stage of the cells we considered to be “fetal” lung cells and organoids generated from hPSC using our most recent lung differentiation protocol¹⁶ (Fig. 8A). We performed scRNA-seq on the hPSC-derived fetal lung cells, considered to contain mostly undifferentiated lung epithelial cells and fetal organoids, which were subjected to the next step along the differentiation pathway.

A Schematics of stage-specific differentiation of PSC towards mature proximal airway epithelia. A graphical image was created with Biorender.com. B UMAP projection of hPSC-derived fetal lung cells. C DEG overlaps between primary fetal cell types and hPSC fetal lung clusters. D UMAP projection of hPSC-derived fetal lung organoids. E DEG overlaps between primary fetal cell types and hPSC organoid clusters. F Spearman’s rank correlation coefficient (red to blue) analysis between hPSC-derived fetal lung models and primary-derived fetal lung epithelia and Travaglini et al. adult lung dataset using the top 50 DEG for each gestational week and adult lung. Source data are provided as a Source Data file.

Full size image

After clustering and performing DEG analysis, we compared DEGs for these hPSC cell clusters to the DEGs found in the primary fetal lung samples. In the hPSC fetal lung cells, we found significant overlap with proliferating progenitors and budtip progenitors in Clusters 0 and 1 (Fig. 8B, C). This was not surprising as the culture conditions were intended for fetal lung cell growth using retinoic acid, CHIR99021, and FGF7³⁶. Other cell type DEGs that emerged in these fetal lung cell cultures included PNEC, SOX2 +, SOX9 + progenitors, and TP cells in Cluster 2. On the contrary, hPSC-derived fetal lung organoids, which represent the next stage of differentiation, contained cells overlapping significantly with basal cells, club cells, PNEC, SMG secretory cells, ciliated cells, and TP cells (Fig. 8D, E). The high basal cell overlap in Clusters 0 and 2 was not surprising as our organoid expansion media contained dual TGFβ/SMAD signaling inhibitors: DMH1 and A83-01, which are intended to promote basal cell expansion⁷⁴.

To benchmark these hPSC differentiated cells to the gestational time point of fetal lung development, we used a Spearman correlation coefficient and found a relatively strong correlation to these early fetal lung tissues. We found that hPSC-derived fetal lung cells correlated the strongest to <GW12 lung epithelia tissue, and hPSC-derived fetal lung organoids correlated the strongest to >GW16 lung epithelia tissue (Fig. 8F). Importantly, none of the hPSC fetal lung models share a strong correlation with adult lung.

Using batch-balanced nearest neighbors, we integrated the hPSC-derived cells and the primary fetal lung epithelia (Fig. 9A). The integrated UMAP showed significant overlap of hPSC-derived fetal lung cells (purple) and organoids (green) with epithelial cells in the primary fetal lung tissue (Fig. 9B). To understand the extent to which these hPSC differentiation represent fetal lung epithelial development, we positioned these cells along our previously inferred trajectories from the native tissue. We considered the previously inferred trajectory for the primary fetal epithelium from budtip progenitors to TP cells, through tip cells, NKX2-1 + SOX9 + CFTR + cells, and SOX2^high/lowCFTR + cells (Fig. 9C). Using the integrated nearest neighbors graph we clustered the integrated cells and only assessed the clusters with at least 100 cells along the trajectory. Therefore, clusters with <100 cells overlapping cells between fetal and hPSC-derived cultures are not shown in the analysis (Fig. 9D).

A Integrated UMAP projection of hPSC-derived fetal lung cells and organoids integrated with primary fetal lung epithelia. B Integrated UMAP projection of hPSC-derived fetal lung cells and organoids cells. C Inferred primary fetal trajectory (black line with arrow) from budtip progenitors to TP cells overlaid on the integrated UMAP embedding (UMAP colored by pseudotime, dark blue to yellow). D hPSC-derived clusters overlapped on integrated UMAP projection. E Cell type scores based on the top 100 DEGs from the primary fetal epithelia reveal high proliferating progenitor scores in the hPSC fetal lung cells and high basal scores in the hPSC organoids. Scores are normalized so that the median score for a cell of its own type is 1. hPSC clusters from left to right increase in average pseudotime. F Trajectory heat map showing the change in gene expression of hPSC fetal lung cells and organoids within these clusters, ordered by increasing average pseudotime of the primary fetal cells within each cluster. G LatentVelo analysis revealing inferred cell state trajectories (arrows) projected onto the UMAP for hPSC fetal lung cells. H LatentVelo analysis revealing inferred cell state trajectories (arrows) projected onto the UMAP for hPSC organoid cells. I Schematics of the inferred trajectory of hPSC-derived cell types found in the fetal lung cells and organoids along the differentiation path (arrow pointing upwards). Source data are provided as a Source Data file.

Full size image

To unbiasedly determine the epithelial cell subtype in the hPSC-derived cultures along the trajectory, we created cell type scores from the top 100 DEG for each cluster identified in the primary epithelial tissue and using the average pseudotime of each of these clusters containing hPSC fetal lung cells and organoid cells, we positioned the clusters along the trajectory (Fig. 9E, arrow above points to the direction of trajectory). These scores indicated that the fetal lung organoids became increasingly differentiated along this trajectory with the acquisition of specific genes associated with various differentiated epithelial cell types/states such as basal and TP cells. Proliferating progenitor cell type scores seen with an early pseudotime in the hPSC fetal lung cells clusters, and early hPSC organoids clusters were lost with increasing pseudotime, while SOX2^high/lowCFTR⁺, TP, and basal cell type scores increased in the organoids. A clear pattern of the change in gene expression emerged along the trajectory sorted by average pseudotime, with the loss of proliferating progenitor marker genes seen in the hPSC fetal lung cells (e.g., TOP2A and CENPF) and acquisition of genes identifying the more differentiated cell types such as GPC5A for SOX2^high/lowCFTR + cells, CP for TP cells, and SPRR3 for basal cells (Fig. 9F).

To further confirm this analysis of the integrated cells, we performed an independent RNA velocity analysis on the hPSC fetal lung cells and hPSC organoids using LatentVelo. Similar to the integrated trajectory analysis, we found trajectories from hPSC fetal lung cluster 0 towards clusters 1 and 2 (Fig. 9G), and trajectories from hPSC organoid clusters 1 and 4 towards clusters 0 and 2 (Fig. 9H).

Overall, our hPSC differentiated fetal cell models captured some of the epithelial cell types/states and trajectories, specifically the budtip to TP cells, along several intermediary states as observed in the primary fetal pseudoglandular/canalicular lung tissue (Fig. 9I).

In summary, we provide a high-resolution spatiotemporal atlas of the human fetal lung with a focus on the lineage relationships of the epithelial cells and interacting stromal cells in forming the developing airways (Fig. 10).

A schematic of the spatial localization of the epithelial and stromal cell subtypes in the developing airways illustrating the spatiotemporal changes in the developing airways. The developing distal bud region is populated by undifferentiated progenitor cell types. As the proximal stalk region develops, proximal progenitors (located in regions between the bud and stalk) generate differentiated cells that populate the developing proximal epithelium with specialized epithelial cell types. Differentiation of the developing epithelia occurs along the distal (right) to the proximal (left) axis. Cell types were identified by scRNA-seq, and their spatial localization was defined by spatial transcriptomics and immunofluorescence staining. The graphical image was created with Biorender.com.

Full size image

• SEO Powered Content & PR Distribution. Get Amplified Today.
• PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
• PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
• PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
• PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
• Source: https://www.nature.com/articles/s41467-024-50281-5