In Vitro and In Vivo Development of the Human Airway at Single-Cell Resolution

Abstract

Bud tip progenitor cells give rise to all murine lung epithelial lineages and have been described in the developing human lung; however, the mechanisms controlling human bud tip differentiation into specific lineages are unclear. Here, we used homogeneous human bud tip organoid cultures and identified SMAD signaling as a key regulator of the bud tip-to-airway transition. SMAD induction led to the differentiation of airway-like organoids possessing functional basal cells capable of clonal expansion and multilineage differentiation. To benchmark in vitro-derived organoids, we developed a single-cell mRNA sequencing atlas of the human lung from 11.5 to 21 weeks of development, which revealed high degrees of similarity between the in vitro-derived and in vivo airway. Together, this work sheds light on human airway differentiation in vitro and provides a single-cell atlas of the developing human lung.

Keywords: basal cell; human development; lung; lung development; lung organoid; organoid; single-cell RNA-seq.

Figure 1.. SMAD activation induces TP63 expression in bud tip progenitor organoids

A) Schematic of lung epithelial development. As the airways extend, bud tip progenitors are maintained as progenitors in the tips of branching buds, and leave cells behind that give rise to the intrapulmonary airways. Late in development, remaining bud tip progenitors differentiate into alveolar cells. No bud tip progenitors are present in the adult lung. B) Schematic of creation of epithelium-only bud tip progenitor organoids from 12 week fetal lungs. C) Feature plots for bud tip progenitor marker genes SFTPC, ID2 and HMGA1 from scRNA-seq of day 0 bud tip progenitor organoids maintained for 3 weeks in culture. Additional feature plots from this same data set are shown in Figure 4A. D) mRNA expression by QRT-PCR of basal cell marker TP63 in bud tip progenitor organoids treated for 3 days with serum-free basal medium supplemented with DMSO (control) or with signaling factors known to be important for lung development and cellular differentiation. DMSO (1:1,000 dilution), FGF7 (10 ng/mL), ‘3F’ (FGF7 10ng/mL, CHIR99021 3μM, ATRA 50 nM), Dexamethasone (25 ng/mL), IL6 (10 ng/mL), IL2 (50 U/mL), Smoothened Agonist (SAG; 500 nM), EGF (100 ng/mL), FGF10 (500 ng/mL), CHIR99021 (2 μM), IFNγ (10 ng/mL), TGFβ1 (100 ng/mL), BMP4 (100 ng/mL), DAPT (10 μM), Hydrocortisone (100 ng/mL). Gene expression is reported as arbitrary units. Treatment with TGFβ1 led to a significant increase in the expression of TP63 (one-way Analysis of Variance (ANOVA) (alpha=0.05, p<0.0001, F=14.7. Dunnett’s test of multiple comparison’s compared the mean of each group to the mean of the DMSO control group.) Estimated p values are shown on the graph. Error bars are plotted to show mean +/− the standard error of the mean. N=3 independent biological specimens. Data is from a single experiment. E) Bud tip progenitor organoids were treated with FGF7 (10 ng/mL), a permissive environment for TP63 expression compared to maintenance in bud tip progenitor medium (‘3F’), or with FGF7 (10 ng/mL) plus factors to inhibit SMAD signaling (A8301 [1 μM] and NOGGIN [100 ng/mL]) and TP63 gene expression was evaluated by QRT-PCR after 10 days in culture. A one-way Analysis of Variance was used followed by Tukey’s multiple comparison test to compare the means of each group to the mean of every other group. Estimated p values are reported on the graph. Error bars are plotted to show mean +/− the standard error of the mean. N=3 independent biological specimens. Data is from a single experiment. F) Bud tip progenitor organoids were treated for 3 days with SMAD activation or inhibition conditions and expression of TP63 was evaluated by QRT-PCR for all treatment groups. All media contained 3F components (FGF7 10ng/mL, CHIR99021 3μM, ATRA 50 nM), with individual groups containing combinations of: DMSO (1:1000 dilution), TGFβ1 (100 ng/mL), BMP4 (100 ng/mL), SB431542 (10 μM), LDN212854 (200 nM). One-way ANOVA alpha=0.05, F=21.19, p<0.0001; Tukey’s multiple comparisons of the mean of each group versus the mean in all other groups, estimated p values are reported on the graph. 3 days TGFβ1 and BMP4 is referred to as ‘dual SMAD activation’, or ‘DSA’. Data is plotted as arbitrary units. Error bars are plotted to show mean +/− the standard error of the mean. N=3 independent biological specimens. Data is from a single experiment and is representative of n=3 experiments. G) QRT-PCR for markers of canonical differentiated lung epithelial cell types showing DMSO (gray bars) and DSA treated (blue bars) organoids after 3 total days of treatment. Data is plotted as fold change over DMSO controls. Two-sided Mann-Whitney Tests were performed to compare the mean of the DMSO group to the Dual Smad Activation group (alpha=0.05). Error bars represent the mean +/− the standard error of the mean. n=3 independent biological specimens, and data is from a single experiment and is representative of n=3 experiments. H) Protein staining of DMSO treated (control) fetal bud tip progenitor organoids (top row) and 3 days of DSA treatment (bottom row) for TP63+ protein (green), KRT5 (pink) and DAPI (blue). Scale bar represents 50 μm. I) Quantification of (g). Total number of TP63+ cells were counted for 3–9 individual organoids across 3 biological replicates. n=3 independent biological specimens. A two-sided Mann Whitney test (alpha=0.05) was used to compare the means of each sample. For all graphs, p values are reported as follows: * p<0.05; ** p<0.01, *** p<0.001, **** p<0.0001.

Figure 2.. SMAD signaling in the native fetal lung and expansion of TP63+ organoids by SMAD inhibition.

A) Protein staining for TP63 (green), phospho-SMAD2 (pSMAD2; pink), E-Cadherin (ECAD; white) and DAPI (blue) of 12 week fetal lungs along the proximal-distal axis. Representative images shown from n=3 independent specimens. Scale bars represent 50 μm. B) Protein staining of 12 week fetal lungs along the proximal-distal axis for TP63 (green) Keratin 5 (KRT5; pink), E-Cadherin (ECAD; white) and DAPI (blue). Representative images shown from n=3 independent specimens. Scale bars represent 50 μm. C) Quantification of TP63+ cells, defined as TP63+ cells within an airway region divided by 100 DAPI+ cells within that airway region, throughout the bronchial tree from 10–20 weeks gestation. 100 DAPI+ cells were counted for each replicate for each region. Biological replicates are plotted using cyan, red, and black dots. Trachea (T), Large Cartilaginous Airway (LC), Small Cartilaginous Airway (SC), Non-Cartilaginous Airway (NC), Bronchiole (BO), Bud Tip (BT). Data is plotted as mean +/− the SEM. D) Overview of experimental design. Bud tip progenitor organoids were treated with 3 days of DSA and then treated with Dual Smad Inhibition (DSI) to expand the organoids. E) Brightfield images of bud tip progenitor organoids at day 0, day 3 of DSA treatment, and after 3 days of DSA and 7 days of DSI (day 10). Scale bar represents 200 μm. F) QRT-PCR for TP63 in bud tip progenitor organoids cultured with 3F bud tip progenitor maintenance medium for 10 days compared to organoids either treated with 3 days of DSA followed by DSI or with 3F directly followed by DSI. Bud tip progenitor organoids were derived from 3 separate 12-week fetal specimens. Data is plotted as the mean +/− SEM. Statistically significant variation in means was calculated using a one-way ANOVA (alpha=0.05) followed by dunnett’s test of multiple comparisons of the mean of each group to the mean of every other group. Data is from a single experiment. Estimated p values are reported as follows: * p<0.05; ** p<0.01, *** p<0.001, **** p<0.0001. G) After 21 days in culture, organoids were interrogated for protein staining of the bud tip progenitor marker SFTPC (green), fetal secretory cells (SCGB3A2 [green], SFTPB [white]), club cells (SCGB1A1 [pink]), goblet cells (MUC5AC [white]), multiciliated cells (AcTUB [white]), basal-like cells (TP63 [green] and KRT5 [pink]), and neuroendocrine cells (CHGA [pink], Synaptophysin [green]). Scale bars represent 50 μm.

Figure 3.. Defining cell signatures of human fetal lung epithelial cells.

A) Schematic of the experimental setup and samples included in the analysis. B) A total of 8,443 EPCAM+ cells were computationally isolated and cellular transcriptome heterogeneity was visualized using a tSNE plot revealing 12 clusters of cells. C) Left panel: boxplots (interquartile range with minimum and maximum, outliers removed from plot) show expression level distributions of canonical cell type markers in each cluster. Right panel: heatmap shows z-transformed cluster average expression levels of top 50 cluster markers ranked by log-transformed fold change in expression levels in each cluster compared to other clusters in any of the clusters. Clusters were identified by expression of markers canonically associated with cell populations based on published literature, or given a new name (e.g. secretory progenitor, bud tip adjacent). A list of the upregulated genes in the bud tip progenitor cluster (BP) and the basal cell cluster (BC) with top fold change in expression levels relative to all other clusters is shown and provides a gene signature for these cell populations. Cell cluster names and colors correspond to those in Figure 1B. D) Feature plots of canonical bud tip progenitor markers SFTPC and AGER, a canonical marker of alveolar epithelial type 1 cells. E) Protein staining for SOX9 (green), SFTPC (pink), and DAPI (blue) in bud tips of a 16 week fetal lung specimen. Protein staining for SOX9 (green), PDPN (pink), RAGE (also known as AGER; white), and DAPI (blue) in cells adjacent to the bud tips in a 16 week fetal lung sample. Scale bars represent 50 μm. F) Feature plots showing expression of canonical basal cell markers TP63 and KRT5, as well as KRT15 and IL33, which are highly specific to cluster 7. G) Protein staining of tracheal TP63+ basal cells in 17 week fetal lungs for other basal cell markers identified in cluster 7 by scRNA-seq: KRT15 (pink), EGFR (pink), IL33 (pink), F3 (white) and PDPN (pink).

Figure 4.. In vitro-derived airway organoid cells recapitulate molecular features of native fetal airway cells.

A) Feature plots showing expression of bud tip progenitor markers SFTPC and SOX9 and basal cell marker TP63 from scRNA-seq of in vitro cells at day 0 (bud tip progenitor maintenance), day 3 (after 3 days of DSA) and day 21 (3 days DSA, 18 days DSI). Additional feature plots of the Day 0 scRNA-seq data are shown in Figure 1C. B) tSNE projections from day 0, day 3 and day 21 in vitro derived cells showingthe sum of z-transformed expression levels of 331 fetal basal cell cluster markers, as identified from human fetal basal cells in Fig. 3, in day 0, day 3 and day 21 in vitro samples. C-E) Heatmap showing Pearson’s correlation coefficients (PCCs) of transcriptome between each of 2,000 randomly selected cell from each time point (day 0 [C], day 3 [D] and day 21 [E]) and the 12 fetal epithelial sub-clusters. PCCs were calculated using log-normalized expression levels of genes that are top 50 markers in any of the fetal epithelial sub-clusters and detected in in vitro cells. Colors represent range of PCCs for all cells in each time point. F) 21 day organoids were subjected to scRNA-seq and visualized by tSNE. Colors represent the fetal cell type identity as described in Fig 3 with the highest PCC to each individual 21 day cell. Top right corner: feature plot showing expression pattern of EPCAM across day 21 cells. C-F) MCP=multiciliated precursor, MC=multiciliated, NE=neuroendocrine, SP=secretory progenitor, BP=bud tip progenitor, BA=bud tip adjacent, BC=basal cell, SMG=submucosal gland, SMB=submucosal basal, GS=goblet-like secretory, IC=intermediate ciliated, CS=club-like secretory.

Figure 5.. Isolated in vitro-derived basal-like cells give rise to clonal airway organoids.

A) Overview of experimental design. B) Protein stain for TP63 (green), EGFR (pink), F3 (white) and DAPI (blue) in organoids treated for 3 days of DSA followed by 7 days of DSI. Scale bar represents 50 μm. C) EGFR/F3 and control FACS plots for 1 representative biological replicate (n=3 biological replicates). D-E) Protein staining for TP63 (red, nuclear) and quantification of TP63+ cells/total cells on sorted cells spun onto a glass slide using a Cytospin. Error bars represent the mean +/− the standard error of the mean. N=3 biological replicates, shown as 3 separate colors on the plot. F) Quantification of the percentage of cells in each organoid containing positive protein staining for canonical cell type markers in unsorted organoids versus organoids that were grown from isolated EGFR+/F3+ cells. N=3 biological replicates per group and n=3 technical replicates (individual organoids) per biological replicate. A total of 198 cells were counted for the EGFR-/F3- group, 110 cells counted for EGFR+/F3- group, 38 cells counted for the EGFR-/F3+ group, and 88 cells counted for the EGFR+/F3+ group. G) Protein staining of organoids derived from EGFR/F3 sorted cells for secretory marker SCGB1A1 (pink), multiciliated markers AcTUB+ and FOXJ1+ (white), epithelial marker KRT8 (green), basal cell marker TP63 (pink), goblet cell marker MUC5AC (white) and DAPI (blue). Data shown from a single biological replicate and is representative of N=3 biological replicates. Scale bars represents 50 μm. H)GFP (green), RFP (red) and brightfield images of whole organoids 11 days after re-plating and mixing EGFR/F3 sorted cells from GFP and RFP expressing groups. N=1 biological replicate with n=6 technical replicates (wells of mixed organoids). Data shown from a single experiment and is representative of n=3 experiments. Scale bar represents 200 μm (top row) and 100 μm (bottom row).