Single-cell transcriptomics reveals lasting changes in the lung cellular landscape into adulthood after neonatal hyperoxic exposure.

Ventilatory support, such as supplemental oxygen, used to save premature infants impairs the growth of the pulmonary microvasculature and distal alveoli, leading to bronchopulmonary dysplasia (BPD). Although lung cellular composition changes with exposure to hyperoxia in neonatal mice, most human BPD survivors are weaned off oxygen within the first weeks to months of life, yet they may have persistent lung injury and pulmonary dysfunction as adults. We hypothesized that early-life hyperoxia alters the cellular landscape in later life and predicts long-term lung injury. Using single-cell RNA sequencing, we mapped lung cell subpopulations at postnatal day (pnd)7 and pnd60 in mice exposed to hyperoxia (95% O2) for 3 days as neonates. We interrogated over 10,000 cells and identified a total of 45 clusters within 32 cell states. Neonatal hyperoxia caused persistent compositional changes in later life (pnd60) in all five type II cell states with unique signatures and function. Premature infants requiring mechanical ventilation with different durations also showed similar alterations in these unique signatures of type II cell states. Pathologically, neonatal hyperoxic exposure caused alveolar simplification in adult mice. We conclude that neonatal hyperoxia alters the lung cellular landscape in later life, uncovering neonatal programing of adult lung dysfunction.

Introduction

Due to advances in neonatal and perinatal care, premature infants can survive at extremes of gestation (>22 weeks). Unfortunately, mechanical ventilation and supplemental oxygen used to save their lives can also impair the growth of their pulmonary microvasculature and distal alveoli. This results in continued dependence on supplemental oxygen beyond the 36 weeks corrected gestational age, and is referred to as bronchopulmonary dysplasia (BPD) [1]. This condition affects 10,000 to 15,000 premature infants annually in the US. The pathology of BPD is characterized by alveolar simplification and abnormal pulmonary vascular development [[2], [3], [4]]. Although most BPD survivors are eventually weaned off oxygen, many show evidences of pulmonary dysfunction and cardiovascular sequelae (e.g., pulmonary hypertension) as adolescents and adults [[5], [6], [7], [8]]. The mechanisms underlying alveolar simplification and dysregulated vascular development as a residual manifestation of BPD are not fully understood and it is not clear whether there is neonatal programing of adult lung dysfunction.

The lung contains more than 40 different cell types [9]. Spatio-temporal interactions between airway epithelial cells and other resident cells (e.g., alveolar epithelial cells, mesenchymal cells, and immune cells) play an important role in alveolar development. In the late fetus, pulmonary vascular and endothelial cell signals are essential for alveolar formation [10,11]. The gas exchange niche in the lung contains two major epithelial cell types, alveolar type I (AT1) and type II (AT2) cells. AT1 cells are thin, squamous cells that are in close contact with pulmonary endothelial cells lining the capillaries, creating a gas exchange region of oxygen and carbon dioxide [12]. AT2 cells are progenitors to ATI cells and also produce surfactant proteins and lipids that reduce surface tension in the alveoli, preventing alveolar collapse (atelectasis) [13]. Although previous studies using genomic, proteomic, and metabolomic analyses from whole lung tissue have furthered our understanding of the mechanisms underlying lung development, these provide limited spatial insights into the complex multicellular environment in development and with injury and repair [[14], [15], [16], [17], [18], [19]]. Single-cell RNA sequencing (scRNA-seq) allows for quantitative and unbiased characterization of cellular heterogeneity by providing genome-wide molecular profiles at the individual cell level. This approach has been applied recently to map cell populations in postnatal lungs from a term one day old infant and from neonatal mice on postnatal day (pnd) 1 to pnd10 [[20], [21], [22]]. These studies identified distinct cell states representing multiple populations of epithelial, endothelial, fibroblast, pericyte, smooth muscle, and immune cells and signature gene markers for each subpopulation. Using this technology, a recent study identified specific AT1 and AT2 progenitors during lung sacculation rather than bipotent cells giving rise to both AT1 and AT2 cells, contradicting previous dogma [23]. Whether and how these progenitor cells may be modified by early life injuries such as hyperoxia is not understood. Most recently, an analysis of the impact of prolonged early life hyperoxia (14 days, equivalent to 1.5 human years) on lung cell populations was performed in terminal experiments [24]. The authors demonstrated altered composition of all cellular compartments including alveolar epithelium and macrophage populations with hyperoxia, and suggested cellular crosstalk and inflammatory signaling as drivers of hyperoxic lung injury. The study did not evaluate whether there was persistence of these alterations once the injury subsided nor whether a shorter exposure could result in similar alterations. Our study is first to evaluate whether neonatal hyperoxic injury results in permanent changes in cellular compartments that explain later lung anatomical and physiologic dysfunction.

Mouse lungs at birth are structurally similar to human neonates born at 30–34 weeks of gestation, when the lung is in the saccular phase of development [25]. In mice, lung alveolarization starts at pnd4. Hyperoxic exposure in newborn mice causes persistent lung injury into adulthood [26]. This model is frequently utilized to investigate pathogenesis and to identify potential therapeutic targets for hyperoxic lung injury in BPD [27,28]. In this study, we employed scRNA-seq to characterize changes in lung cell subpopulations at pnd7 and pnd60 in mice exposed to hyperoxia for 3 days as neonates (<12 h old). We focused on changes in subpopulations of AT1 and AT2 cells as well as identification of potential biomarkers for these subpopulations to determine whether neonatal hyperoxia results in a program leading to adult lung dysfunction.

Methods

Mice and hyperoxic exposure

Newborn C57BL/6J mice (<12 h old, male and female) along with their mothers were exposed to room air or hyperoxia (>95% O2) for 72 h in an A-chamber (BioSpherix, Parish, NY) [29,30]. The dams were switched every 24 h between room air and hyperoxia to avoid injury. Following the 72 h of exposure, these pups were allowed to recover in room air until pnd7 and pnd60. No deaths occurred in mice due to hyperoxic toxicity. All animal experiments were reviewed and approved by the Institutional Animal Care and Use Committee of Brown University (IACUC: 1507000150).

Lung dissociation into single cells

Mice were anesthetized with ketamine and xylazine, the heart/lungs were exposed and perfused with 10 ml of sterile PBS via the right ventricle. Lung cells were dissociated as previously described [31]. In brief, lungs were inflated with 50 U/ml dispase in PBS (400 μl per 10 g body weight). Then 400 μl 1% low-melt agarose in PBS was added and allowed to set on ice. Lungs were incubated in 2 ml dispase for 45 min on a shaking platform, then incubated in a 37 °C water bath for 10 min. Lungs were shredded in a 10 cm dish with 7 ml complete DMEM/F-12 plus 50 U of DNase, using forceps and then placed on an orbital shaker at room temperature for 10 min. The lung tissue was washed through a 100 μm filter and then through the 40 μm filter. The dissociated lung tissue was rinsed thoroughly with medium, thereby generating a single cell suspension. Single cell suspensions from three mouse lungs were pooled for each condition. All lung samples were dissociated swiftly and at the same time to avoid batch effects, and used fresh and immediately for scRNA-seq.

scRNA-seq

Single cell encapsulation was performed using the Chromium Single Cell Chip G kit on the 10 × Genomics Chromium Controller [32,33]. Single cell cDNA and libraries were prepared using the Chromium Single Cell 3′ Reagent Kit v3.1 Chemistry. Libraries were sequenced by Genewiz on the Illumina Hiseq (2 × 150 bp paired-end runs). Single cell unique molecular identifier (UMI) counting (counting of unique barcodes given to individual transcript molecules), was performed using Cell Ranger Single Cell Software Suite 3.0.2 from 10 × Genomics. The Cell Ranger pre-build mouse transcriptome reference (mm10) was used for the analysis. Cell Ranger gene expression matrices were further analyzed using the R package Seurat v 3.2.1 [34]. The data were filtered using a mitochondrial cut off at 5%. Cells with at least 700 and at most 8000 expressed genes (features) were included in the downstream analysis. Using SCTransform, the four datasets were normalized and integrated to identify conserved cell populations across the samples [35]. The integration step was performed using 3000 features, the shared nearest neighbor graph was constructed using 50 dimensions, and clusters were identified using a resolution of 2. t-SNE dimensionality reduction was computed using 50 dimensions.

Computational analysis of scRNA-seq data

Cell states were identified by comparing our datasets to mouse cell atlas lung data using Seurat functions FindTransferAnchors and TransferData [34,36]. The mouse cell atlas data was also normalized using SCTransform and transfer anchors were determined using the mouse cell atlas data as the ‘reference’ and the experimental data as the ‘query’.

To indicate which genes were distinguishing different cell states, a Wilcoxon Rank Sum test was implemented in the Seurat function FindAllMarkers. This function compares each individual ‘identity’ in the Seurat object to all other identities in the object. The test was run on SCTransform normalized counts using the Seurat clusters as the identities in the Seurat object. To further determine specific gene markers for AT2 subpopulations, cluster-specific markers for clusters 3, 11, 13, 22, and 27 were tabulated by comparing the non-overlapping marker genes found in each cluster (e.g., cluster 3 marker genes not found in clusters 11, 13, 22, or 27).

A two-way ANOVA (anova_test function in rstatix) of the SCTransformed counts indicated an interaction between time points (pnd7 and pnd60) and exposures (Air and O2) (P < 0.05), so the FindAllMarkers function in Seurat was ran to further explore condition-specific markers. FindAllMarkers was run on the SCTransformed counts to discern condition-specific markers by using the experimental conditions as the identities in the Seurat object.

Lung tissues from premature infants

Human lung samples were obtained from premature infants between 23- and 29-weeks postmenstrual age, who lived 5–15 days and required mechanical ventilation (short-term), and controls were premature infants who were not mechanically ventilated and survived less than 24 h. Other human samples (long-term) consisted of preterm infants born at 23–29 weeks who had lived for more than 6 weeks, and had been ventilated for most of their life. Controls for this group consisted of stillborn term infants or age matched infants who had not been ventilated. This was described in a previous report [4]. Utilization of human lung samples was done in compliance with the Institutional Review Board guidelines of Women and Infants Hospital, Providence RI.

Immunofluorescence

Lung sections were de-paraffinized, rehydrated and subjected to heat-mediated antigen retrieval in a citrate buffer solution (Vector Labs). Samples were then exposed to hydrogen peroxide (3%) for 30 min to quench endogenous peroxidase activity. Non-specific binding was blocked by incubating the sections with 5% normal goat serum in PBS for 30 min. Lung section were stained overnight at 4 °C with antibodies against an AT2 cell marker pro-surfactant protein C (pro-SPC, sc-518029, Santa Cruz, 1:50 dilution), cystic fibrosis transmembrane conductance regulator (Cftr, ab59394, Abcam, 1:100 dilution), transketolase (Tkt, NBP1-87441, Novus Biologicals, 1:100 dilution), SPOCK2 (BS-11966R, ThermoFisher, 1:100 dilution) and 5′-aminolevulinate synthase 2 (Alas2, BS-9516R, ThermoFisher, 1:100 dilution). After incubation with secondary antibodies for 2 h at room temperature, sections were mounted in hard-set mounting medium containing DAPI (Vector Labs) and allowed to incubate overnight. Images were taken using a Zeiss Axiovert 200 M Fluorescence Microscope. Co-localization of pro-SPC with an interest target were counted in three randomly selected high-power fields (HPF) each sample, which was normalized into numbers of DAPI+ nuclei. These experiments were carried out in a blinded manner.

RNA in situ hybridization

RNAScope technology (ACDBio, Newark, CA) was used to perform Malat1 RNA in situ hybridization according to manufacturer's instructions. In brief, lung sections were deparaffinized using xylene following by incubation with hydrogen peroxide to block endogenous peroxidase activity. The slides were incubated with heat-mediated co-detection target retrieval solution to unmask the target RNA and protein. Lung sections were stained overnight at 4 °C with a pro-SPC antibody (AB3786, Millipore, 1:200 dilution). After fixing with 10% neutral buffered formalin and incubation with RNAscope protease plus, a negative control probe (Dap8), a positive control probe Ppib, and a Malat1 probe were added to lung slides with incubation at 40 °C for 2 h using the HybEZ™ oven (ACDBio, Newark, CA). Signals were amplified by three amplifiers with sequential hybridization AMP1 to AMP3. The amplified signal was detected using Opal 650 (1:1500 dilution) followed by incubation with Alexa Fluor-conjugated anti-rabbit secondary antibody 488 (1:250 dilution) for 30 min at room temperature. After the staining was completed, slides were washed with PBS containing 0.1% Tween 20, and mounted with Prolong Gold antifade mounting medium (Thermo Fisher Scientific). Images were taken using a Zeiss Axiovert 200 M Fluorescence Microscope. Co-localization of pro-SPC with Malat1 were counted, which was normalized into numbers of DAPI+ nuclei. These experiments were carried out in a blinded manner.

Lung morphometry

Mean linear intercept (Lm) and radial alveolar count (RAC) were measured in mouse lungs stained with hematoxylin and eosin (H&E) as previously described [29]. In brief, we inflated non-lavaged mouse lungs with 1% low melt agarose at a pressure of 25 cm H2O, and fixed them with 4% neutral buffered paraformaldehyde. Lung midsagittal sections with H&E staining were utilized to determine Lm. A perpendicular line was drawn from the center of the respiratory bronchiole to the distal acinus (as defined by the pleura or the nearest connective tissue septum). The number of septae intersected by each line was counted as RAC, and a minimum of 5 counts were performed for each animal.

Statistical analysis

Statistical analyses of immunofluorescence and lung morphology were performed using GraphPad Prism 9. The results were expressed as mean ± SEM. The student t-test was used for detecting statistical significance of the differences between means of two groups if the data are normally distributed. Welch-corrected t-test was used if the data are not normally distributed. The statistical significance of the differences among multiple groups was evaluated by using one-way ANOVA for overall significance, followed by Tukey-Kramer test. Statistical significance was considered when the p value was less than 0.05.

Results

Neonatal hyperoxic exposure alters the lung cellular landscape in later life

To generate a cell-type resolved map of lungs from air- and hyperoxia-exposed mice, we performed highly parallel genome-wide expression profiling of individual cells using the Drop-seq workflow (Fig. 1). We employed scRNA-seq to detect lung cell subpopulations at both pnd7 and pnd60 in mice exposed to hyperoxia as neonates and compared them to air-exposed controls. Cell Ranger gene expression matrices were further analyzed using the R package Seurat v3. To exclude low quality events, droplets with ≤500 or ≥8,000 genes, or with high mitochondrial transcripts (>5%) were removed [37]. This filtering resulted in an analytical dataset of 6,705 and 6,001 cells in air and hyperoxia groups, respectively, with an average detection of 1,200 genes per cell. The total numbers of lung cells detected in Air/pnd7, O2/pnd7, Air/pnd60, and O2/pnd60 groups were 2308, 1246, 4397 and 4755, respectively. This suggests that neonatal hyperoxia causes the loss of lung cells at pnd7 not at pnd60. The datasets from these 4 conditions (Air/pnd7, O2/pnd7, Air/pnd60, and O2/pnd60) were integrated to identify cluster markers across these conditions. We identified a total of 45 clusters (cell states) of lung cells from their corresponding marker genes (Fig. 2A, Supplementary Table 1). The top marker genes in each cluster was shown in Supplementary Fig. 1. This dataset contains 32 cell types with the four major cell types grouped as epithelial cells, endothelial cells, mesenchymal cells, and immune cells (Fig. 2B–C, Table 1, Supplementary Table 2). Within the epithelial cell types, AT1, AT2, ciliated, and Clara (Club) cells were identified, while three types of stromal cells with high expression of Acta2, Dcn, and Inmt were identified. Alveolar macrophage Ear2_high cells, T cell_Cd8b1 (A, B), dividing T cells and eosinophil granulocytes were major immune cells identified (Fig. 2B, Supplementary Table 2). These results confirm known lung cell heterogeneity but from an unbiased framework of gene expression.

Fig. 1

Experimental design. Whole lung single-cell suspensions from neonatal mice exposed to air or hyperoxia for 3 days. The samples collected after air recovery at both pnd7 and pnd60 were analyzed using the Drop-seq workflow.

Fig. 2

Drop-seq analysis identifies a diversity of cell types in mouse lungs. (A) t-SNE plot representing the integration of the four conditions (Air/pnd7, O2/pnd7, Air/pnd60 and O2/pnd60). Single cells are colored by cluster identity. Forty-five clusters were detected across the four conditions. (B) The dot plot shows (1) the percentage of cells expressing respective selected marker gene using dot color and (2) the average expression level of each gene based on unique molecular identifier (UMI) counts. Rows represent clustered cell types, demonstrating similarities of transcriptional profiles. (C) t-SNE plot representing the integration of the four conditions (Air/pnd7, O2/pnd7, Air/pnd60 and O2/pnd60). Single cells are colored by cell type identity. Thirty-two cell types are detected across the four conditions. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Table 1

Cell subpopulations and categories.

Type	Cell status	Air/pnd7	O2/pnd7	Air/pnd60	O2/pnd60
Epithelial cells	AT1 Cell	63	46	40	55
	AT2 Cell	432	234	453	530
	Ciliated cell	67	7	42	38
	Clara Cell	60	12	6	7
Endothelial cells	Endothelial cell_Kdr high	26	31	117	134
	Endothelial cell_Tmem100 high	227	154	618	644
	Endothelial cells_Vwf high	20	24	77	64
Stromal cells	Stromal cell_Acta2 high	210	98	40	46
	Stromal cell_Dcn high	28	16	21	34
	Stromal cell_Inmt high	157	134	281	288
Immune cells	Alveolar macrophage_Ear2 high	156	74	170	187
	Alveolar macrophage_Pclaf high	8	9	15	9
	Interstitial macrophage	24	19	67	56
	B Cell	482	152	1265	1274
	Ig−producing B cell	0	0	3	5
	Conventional dendritic cell_Gngt2 high	60	46	151	194
	Conventional dendritic cell_H2-M2 high	0	4	7	8
	Conventional dendritic cell_Mgl2 high	2	3	39	40
	Conventional dendritic cell_Tubb5 high	20	10	6	11
	Dendritic cell_Naaa high	7	5	34	40
	Dividing dendritic cell	0	0	6	5
	Plasmacytoid dendritic cell	30	21	126	162
	Basophil	11	8	11	11
	T Cell_Cd8b1 high	125	59	464	526
	Dividing T cells	22	17	13	19
	Neutrophil granulocyte	4	3	4	3
	NK Cell	16	10	181	189
	Nuocyte	16	18	109	114
	Eosinophil granulocyte	11	11	27	56
Progenitor	Monocyte progenitor cell	0	0	0	1
	Alveolar bipotent progenitor	3	1	2	0
Other	Dividing cell	21	20	2	5

Identifying specific changes in lung cell subpopulations caused by hyperoxic exposure in neonatal life

Using these four datasets, we looked at the changes in cell state abundance among Air/pnd7, O2/pnd7, Air/pnd60 and O2/pnd60 groups (Fig. 3, Fig. 4, Supplementary Figs. 2–3). Overall, the most abundant detected cell states were B cells, AT2 cells, endothelial cells (Tmem100 high), T cells (Cd8b1 high), stromal cells (lnmt high and Acta2 high) and alveolar macrophages (Ear2 high).

Fig. 3

t-SNE plots of lung cell subpopulations. t-SNE plots represent 45 clusters in four integrated conditions: Air/pnd7, O2/pnd7, Air/pnd60 and O2/pnd60.

Fig. 4

Percentage of cell states detected in each condition. Percentage of each cell status in Air/pn7, O2/pnd7, Air/pnd60 and O2/pnd60 groups.

Importantly, cell states that were most down-regulated by 3 days of hyperoxia vs air, at pnd7, were B cells, Club cells and ciliated cells. However, these cell states returned to normal levels vs air, by pnd60 (Fig. 5A and B). In contrast to the Club cells and the ciliated cells that were each represented by a single cluster (39 and 28 respectively), several clusters were identified as B cells (1, 2, 4, 6, 9, 21) (Fig. 5C). Overall, despite significant decreases at pnd7, cumulative B cells abundance across all representative clusters approximated air exposed values by pnd60. When evaluating differences in individual B cell clusters, cluster 4 was the most decreased after hyperoxia on pnd7 but reached comparable low levels at pnd60 vs air exposed. At pnd7, there were no significant differences in clusters 1, 2, 6, 9, and 21 by hyperoxia. At pnd60, cluster 1 was markedly upregulated while clusters 2, 6, 9 and 21 were reduced by hyperoxia (Fig. 5C and D). Cluster 1 has unique markers of mature regulatory B cells (CD19 and IgM), suggesting that hyperoxia enhances B cell maturation. These results suggest that neonatal hyperoxia causes persistent and specific changes in lung B cell subpopulations.

Fig. 5

Effect of hyperoxia on B cells, Club cells, ciliated cells, AT1 and AT2 cell heterogeneity. (A) Fold change in percentage of B cells, ciliated cells and Club cells in O2/pnd7, Air/pnd60 and O2/pnd60 groups compared to Air/pnd7. (B) Numbers of cells detected in each cell population per condition. (C) Fold change in percentage of clusters 1, 2, 4, 6, 9, and 21 in B cells from O2/pnd7, Air/pnd60 and O2/pnd60 groups compared to Air/pnd7. (D) Numbers of cells detected in each cluster per condition in B cells. (E) Fold change in percentage of AT1 and AT2 clusters in O2/pnd7, Air/pnd60 and O2/pnd60 compared to Air/pnd7. (F) Numbers of cells detected in each cell population per condition. (G) Fold change in percentage of clusters 3, 11, 13, 22, 25 and 27 in O2/pnd7, Air/pnd60 and O2/pnd60 groups compared to Air/pnd7. (H) Numbers of cells detected in each of these clusters per condition.

Lung AT1 and AT2 heterogeneity in mice with postnatal development

The lung cellular architecture changes as the animal matures [22]. We therefore investigated the states of AT1 and AT2 cells and their modulation with postnatal lung development. Four clusters of AT1 cells and twelve clusters of AT2 cells were identified, including three mixed AT1/AT2 cell clusters (12, 22 and 27) and one mixed alveolar macrophage/AT2 cell cluster 4 (Fig. 5E and F, and Supplementary Table 2). Of these, clusters 12, 29, 30, 39, 40 and 44 had too few numbers and were not further considered in our analyses.

As shown in Fig. 5G and H, under normoxia, the proportion of cells in clusters 3 increased at pnd60 compared to pnd7 suggesting amplification of this AT2 cell state with maturation. In contrast, cell proportions in AT2 clusters 11, 13, 22 and 27 were reduced by postnatal development under normoxic condition, as was the pure AT1 cluster 25, suggesting that these cells either matured to other cell states or were lost altogether. It is interesting to note that cluster 22 had both markers for AT1 and AT2 cells. The proportion of cells detected in this cluster was also highly reduced with postnatal development, which may suggest terminal differentiation into mature AT1 or AT2 cells, inability of some AT2 cells to serve as progenitors to AT1 cells, or the death of this cell state with postnatal lung development. Since the pure AT1 cell cluster 25 did not significantly change between air and hyperoxia at pnd60, we focused the remaining of our analyses solely on AT2 cell states.

We then wanted to better understand whether hyperoxia would perturb specific AT2 lung cell states, thereby altering lung architecture. Most AT2 clusters were increased by hyperoxia at pnd7, with the exception of clusters 22 and 27, which were decreased. The increase observed at pnd7 was reversed at pnd60 in cluster 3, whereas cluster 11 and, to a lesser extent, cluster 13 remained elevated (Fig. 5G and H). Cluster 22 was reduced at pnd7 and nearly eliminated at pnd60 after hyperoxia. Cluster 27 was eliminated by hyperoxia at pnd7 and this loss persisted at pnd60 (Fig. 5G and H). These results suggest that there is heterogeneity of AT1 and AT2 cell states during lung development, and that hyperoxic exposure differentially influences these subpopulations by promoting persistent proliferation of certain states and persistent suppression of others.

AT2 clusters show distinct and dynamic gene expression

To further identify the role of hyperoxia on modifications in AT2 cell states, and the implication this may have in lung function and architecture, cluster-specific markers for clusters 3, 11, 13 and 27 were tabulated by comparing non-overlapping marker genes found in each cluster (Supplementary Table 3). Compared to other clusters of AT2 cells, cluster 3 cells selectively expressed biomarkers of epithelial progenitor cells (Erbb3, Fgf1, Foxp2, Klf9, Atf3, Bmp4, Cdh16, Gata6, and Pard3) [23], suggesting a progenitor cell state. Importantly, this cluster also selectively expressed genes encoding ion channels including Cftr, Kcne2, Kcnj2, and Trpm6 which are also characteristic of lung epithelial progenitor cells [38,39]. Because this cluster was persistently reduced with hyperoxia at pnd60, it suggests that there was permanent loss of important AT2 cell progenitors needed for epithelial regeneration. Loss of lung Cftr+/pro-SPC+ cells (AT2 cells) further confirmed this in hyperoxia-exposed mice and in premature infants requiring mechanical ventilation (Fig. 6, Fig. 7).

Fig. 6

tSNE plots showing AT2 cell subpopulations and their unique gene expression. (A) tSNE plots show AT2 cell subpopulations in clusters 3, 11, 13, 22 and 27. (B) tSNE plots overlaid with SCTransform-normalized expression values of Cftr, Malat1, Tkt, Spock2 and Alas2 in AT2 cells among Air/pnd7, O2/pnd7, Air/pnd60 and O2/pnd60 groups.

Fig. 7

Expression of Alas2, Cftr, Spock2 and Tkt in AT2 cells from mice exposed to hyperoxia. C57BL/6J mice (<12 h old) were exposed to hyperoxia for 3 days, and allowed to recover in room air until pnd7 and pnd60. (A) Double immunofluorescence was performed to detect co-localization of Alas2, Cftr, Spcok2 or Tkt with pro-Spc (Spc) in lungs, and representative images are shown for each group. Arrows denote co-localization of Alas2, Cftr, Spcok2 or Tkt with pro-SPC. (B) Numbers of co-localized cells were counted in three randomly selected high-power fields (HPF) in each mouse. This was normalized to Dapi+ cells. Bar size: 50 μm. N = 3 mice including 2 males and 1 female per group. (C) Immunofluorescence was performed to detect Cftr+/pro-SPC+ cells in lungs of premature infants requiring mechanical ventilation. Numbers of Cftr+/pro-SPC+ cells were counted in three randomly selected HPF per sample. This was normalized to Dapi+ cells. Bar size: 50 μm. Data are expressed as mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001 vs respective air groups or control group. †P < 0.05, ††P < 0.01, †††P < 0.001 vs respective pnd7 groups.

The long non-coding RNA Malat1 is indispensable for mouse development [40,41] and is overexpressed in various lung cancers [42]. Malat1 is able to activate Nrf2, therefore protecting against oxidative injury to cells [43]. A small nuclear RNA Rn7SK is a noncoding RNA that plays a major role in regulating eukaryotic transcription [44] and is important in macrophage polarization for wound healing in the lung and other tissues [45]. Malat1 and Rn7sk are unique markers of Cluster 11. This cluster remained markedly increased after hyperoxia, suggesting that this may drive persistent overabundance of macrophages and a regional proliferative state in the lung alveoli. As shown in Fig. 8A and B, Malat1 was present abundantly in the mouse lung and Malat1+/pro-SPC+ cells were increased at air pnd60 compared to air pnd7, suggesting a developmentally regulated process. A further persistent increase of this AT2 cell cluster after hyperoxia at pnd60 vs air suggests that the lung remains primed for enhanced proliferation and wound healing after early exposure to hyperoxia, which could lead to increased potential for tumorigenesis and/or scarring in later life due to the resistance to oxidative damage. This remains to be determined. Similarly, Malat1+/pro-SPC+ cells were significantly increased in the lungs of premature human infants requiring mechanical ventilation compared to non-ventilated controls (Fig. 8C and D). This may indicate that assisted ventilation may alter the proliferative potential of premature lungs.

Fig. 8

Increased Malat1 expression in AT2 cells of mice exposed to hyperoxia as neonates and in premature infants requiring mechanical ventilation. (A, B) C57BL/6J mice (<12 h old) were exposed to hyperoxia for 3 days, and allowed to recover in room air until pnd7 and pnd60. RNAscope and immunofluorescence were performed to detect co-localization of Malat1 with pro-Spc in lungs. Representative images are shown for each group. Arrows denote co-localization of Malat1 with pro-SPC. Numbers of cells with co-localized Malat1 and pro-Spc were counted in three randomly selected high-power fields (HPF) for each mouse. This was normalized into Dapi+/cells. Bar size: 50 μm. N = 3 mice including 2 males and 1 female per group. (C, D) Immunofluorescence was performed to detect Malat1+/pro-SPC+ cells in lungs of premature infants requiring mechanical ventilation. Numbers of Malat1+/pro-SPC+ cells were counted in three randomly selected HPF per sample. This was normalized to Dapi+ cells. Bar size: 50 μm. N = 3 per group. Data are expressed as mean ± SEM. *P < 0.05, ***P < 0.001 vs respective air groups (B) or control subjects (D). †P < 0.05 vs respective pnd7 groups (B) or short-term ventilation group (D).

Custer 13 cells uniquely expressed genes involved in mitochondrial homeostasis, lipid metabolism, and the pentose phosphate pathway. This included Slc16a7 (metabolite transport), Slc25a39 (mitochondrial carrier protein), Higd1a (mitochondrial inner protein), Sgpp2 (sphingosine-1-phosphate phosphatase), Pcyt2 (phospholipid), Sqle (sterol biosynthesis), Dhcr24 (cholesterol biosynthesis), Fdft1 (cholesterol biosynthesis), Cyp51(sterol biosynthesis), Hmgcr (cholesterol biosynthesis), Srebf2 (cholesterol biosynthesis); Echs1 (β-oxidation), Acaa1a (lipid metabolism), and Tkt (pentose phosphate pathway). The latter regulates NADPH availability and therefore is important in regulating cellular redox homeostasis. This cluster of cells also selectively expressed genes involved in homeostasis of peptides and extracellular matrix, such as Hdc (forming histamine), Itih4 (stabilizing extracellular matrix), Tmprss4 (serine protease family), DDP9 (serine protease), Usp33 (ubiquitin Specific Peptidase 33), Col6a1 (collagen), and Cndp2 (histidine metabolism). These results suggest that cells in cluster 13 are metabolically active and serve to maintain matrix homeostasis. The proportion of this cell state increased at pnd7 and remained increased at pnd60 after neonatal hyperoxia, suggesting persistent active repair of the lung matrix into adulthood. This was not however fully corroborated by visualization of TKT+/pro-SPC+ cells by immunofluorescence, as these did not significantly increase at pnd60 after neonatal hyperoxia compared to air exposed controls at pnd60, as shown in Fig. 6, Fig. 7.

Cells from cluster 22 selectively expressed Cavin1, Pmp22, Vegfa, Spock2, Clic5, and Rtkn2, biomarkers expressed in AT1 cells [46], suggesting a mixed AT1/AT2 phenotype. This cell state diminished with maturation and was not altered further by neonatal hyperoxia at pnd60, likely indicating either cell death or terminal differentiation of these cells to other cell states. The loss of Spock2+/pro-SPC+ cells at pnd60 was documented by immunofluorescence, as shown in Fig. 6, Fig. 7, corroborating this loss of the AT1/AT2 cell state.

Cluster 27 uniquely expressed Alas2 (heme synthesis), Igfbp5 (epithelial-mesenchymal transition), Bpgm (control glycolytic intermediates), Bnip3l (apoptotic protector), and Gpc3 (a co-receptor or a receiver for Wnt). This cluster was prematurely lost after hyperoxia at pnd7 but was also lost with maturation (pnd60). In non-small cell lung carcinoma cells, inhibition of heme synthesis and mitochondrial function mediated by Alas2, suppresses proliferation, colony formation and migration [47]. Perhaps early loss of Alas2+/pro-SPC+ cells in hyperoxia impacts lung architecture by reducing the proliferative potential of key AT2 cells, in contrast to the proliferative advantage of others (cluster 11). This was corroborated in mouse lung slices where Alas2+/pro-SPC+ cells were reduced in pnd60 after neonatal hyperoxia (Fig. 6, Fig. 7).

Neonatal hyperoxia causes persistent lung alveolar simplification

To determine whether the persistent alterations in lung cellular landscape is associated with lung injury in mice exposed to hyperoxia as neonates, we measured Lm and RAC, parameters for lung injury and growth respectively. At pnd7, lungs of hyperoxia exposed mice showed a trend but no significant increase in Lm nor decreased RAC (Fig. 9). Neonatal hyperoxia increased Lm and reduced RAC in mice at pnd60 (Fig. 9). These results suggest that neonatal hyperoxia dictates persistent alveolar simplification in later life.

Fig. 9

Neonatal hyperoxia causes persistent lung injury in mice. (A) C57BL/6J neonatal mice (<12 h old) were exposed to air or hyperoxia (95% O2) for 3 days, and then allowed for recover in room air until pnd7 and pnd60. H&E staining was performed to measure mean linear intercept (Lm) and radical alveolar count (RAC) in mouse lungs. Bar size: 200 μm. Data are expressed as mean ± SEM. N = 5 mice including 3 males and 2 females per group. **P < 0.01, ***P < 0.001 vs respective air group. (B) Schematic figure showing hyperoxia-induced long-term changes in AT2 cell subpopulations and functions, leading to alveolar simplification.

Discussion

Several studies have identified lung cell subpopulations from neonatal mice and newborn infants under normoxic conditions [[21], [22], [23]]. Dynamic populations of lung cells are also observed in neonatal mice exposed to hyperoxia (up to 14 days) without air recovery using scRNA-seq [24]. When correlating life spans, 14 mouse days is approximately equivalent to 1.5 human years [48]. In the clinical circumstance, most infant survivors of BPD are weaned from oxygen in the first few months of life, nevertheless, they may suffer long-term consequences into adulthood. We explored whether a shorter hyperoxic exposure as a neonate persistently alters the lung cellular landscape in later adult life. This has considerable relevance clinically in devising early protective strategies to prevent this long-term impact. Here, we exposed newborn mice to hyperoxia for 3 days and allowed the animals to recover in room air for up to 60 days to determine whether hyperoxia in neonatal life results in programing of adult lung dysfunction through persistent effects on lung cell subpopulations using scRNA-seq. A total of 45 clusters were identified representing 32 cell types including epithelial cells, endothelial cells, mesenchymal cells, and immune cells. In particular, four clusters of AT1 cells and twelve clusters of AT2 cells were identified, including three mixed AT1/AT2 cell clusters (12, 22 and 27) and one mixed alveolar macrophage/AT2 cell cluster (cluster 4). We observed that cell numbers and gene transcription were altered with postnatal lung development, and that neonatal hyperoxia had persistent cluster specific effects which modified normal lung maturation and lung cell type composition. Neonatal hyperoxia has been shown to cause apoptosis and reduce stem cell function in the lung [49,50], which may contribute to the loss of lung cells at pnd7. In contrast, repair mechanisms after hyperoxic lung injury may persist into adulthood. This could explain no significant changes in total number of lung cells at pnd60 between air and hyperoxia groups.

The top five most abundant cells identified in this scRNA-seq analysis were B cells, AT2 cells, endothelial cells (Tmem100 high), T cells (cd8b1 high), and stromal cells (inmt high). A previous study using a FACS isolation protocol in neonatal lung samples demonstrated a higher proportion of epithelial cells [51]. Similarly, we found approximately 19% AT2 cells detected at pnd7 under normoxic conditions. This is corroborated by a recent report showing 16.8% AT2 cells in newborn mouse lungs [22]. B cells represented 20.8% of the total cell population in the present study, whereas there were only 1% B cells identified in newborn mouse lungs at pnd1 in a previous study [22]. These discrepancies may be due to the difference between lungs at the saccular stage (pnd1) versus at the alveolar stage (pnd7) of development. Recently, B cells were seen to be significantly altered by a prolonged neonatal exposure to hyperoxia [24], however no data were provided as to whether these changes persisted in recovery.

Previous studies have demonstrated that sample preparation, storage, and processing may alter gene expression signatures in scRNA-seq datasets [52,53]. We did not store cells before processing for that reason. Further studies are warranted to optimize these factors so as to acquire more cells for profiling scRNA-seq expression. The objective is to obtain no more than 10,000 cells since this could yield undesirable multiplet rate with the Chromium single cell chip.

Integrated growth factor and BMP signaling pathways promote lung epithelial progenitor cell expansion and maintenance [54]. Genes involved in these pathways were enriched in cluster 3 of AT2 cells. Neonatal hyperoxia increased the relative number of this cluster at pnd7. This may indicate lung AT2 cell expansion in mice exposed to hyperoxia as neonates [55]. In contrast, neonatal hyperoxia reduced the relative number of this cluster at pnd60, suggesting a long-term effect of neonatal hyperoxia and permanent attrition of this progenitor cell line. It is yet unclear whether a second injury would result in worsened outcome if lung epithelial progenitor expansion was curtailed permanently. It is known that depletion of AT2 cells leads to poor repair [56]. Interestingly, this cluster of cells expressed genes encoding ion channels including chloride (Cftr), potassium (Kcne2, Kcnj2), and manganese (Trpm6) channels. Indeed, as a genome-wide association study gene, Kcne2 was shown to be important for maintaining normal lung function [57]. In addition, embryonic stem cells-derived AT2 cells also express Cftr [58]. Persistent loss of this cluster, and thus Cftr expression with hyperoxia would bode poorly for maintenance of normal lung fluid balance and bacterial clearance even in adulthood, which could contribute to sustained tissue damage [59,60].

Cluster 11 uniquely expressed Malat1 and Rn7sk, both of which are non-coding RNA genes. Malat1 is a highly conserved nuclear non-coding RNA. It modulates mRNA and miRNA transcription by regulating alternative splicing and histone translational modification, as well as facilitating transcription factors to bind the promoter region of genes. Malat1 is ubiquitously expressed in all cells during early embryonic stages, while in adult tissues Malat1 expression varies among tissue types with high expression in the lung, brain, and pancreas [61,62]. Previous studies have shown that Malat1 expression is increased in peripheral blood samples from BPD patients and in rodent lungs of BPD models [63,64]. Malat1 was also increased in AT2 cells by neonatal hyperoxia. In several studies, knockdown of Malat1 significantly impaired tumor migration and invasion in vitro and in vivo [42,65]. It is unclear whether permanent expansion of this Malat1 expressing AT2 cell cluster after hyperoxia would disrupt redox homeostasis and promote or predict later lung tumorigenesis. This remains to be explored.

Using bulk mRNA sequencing, dynamic gene expression was analyzed from E16.5 to pnd28 in mouse lung tissue [22]. This revealed that surfactant synthesis and establishment of alveolar fluid balance increased during this developmental period. This was driven by AT2 cells, which are responsible for surfactant synthesis and require lipid metabolism. In our study, cluster 13 cells selectively expressed genes involved in lipid metabolism and matrix homeostasis. Tkt connects the pentose phosphate pathway to glycolysis, and contributes to the production of NADPH as a reducing equivalent. Indeed, the pentose phosphate pathway is activated and a redox cofactor NADPH is increased in lungs of mice exposed to hyperoxia as neonates [29]. Increased Tkt may actively engage fatty acid synthesis in AT2 cells by producing NADPH. Absolute numbers of this cluster increased during lung development. This may contribute to phospholipid synthesis to prevent alveolar collapse during lung mesenchymal expansion. The persistent expansion of this cluster after neonatal hyperoxia also suggests a protective response to matrix injury and maintenance of surfactant production but this expansion could also lead to excessive scarring if not balanced by enhanced epithelial cell proliferation.

It has been shown that the rs1245560 variant of the Spock2 gene is associated with BPD susceptibility in populations of European and African descent, and its expression level is correlated with alveolarization of normal rat lung [66]. Overexpression of Spock2 in mice leads to altered lung alveolar development and worsens lesions induced by hyperoxia [67]. A previous study also identified Spock2 as a marker of early transdifferentiation from AT2 to AT1 cells [68]. Spock2 was uniquely identified in cluster 22, which had a mixed AT1/AT2 genetic phenotype, and hyperoxia reduced Spock2 abundance in AT2 cells. A mixed AT1/AT2 phenotype was described by another group, albeit with different markers (Muc1, Pdpn, Egfr and Shh) [22]. This cluster was significantly reduced in adult mice compared to neonatal mice, suggesting less AT2/AT1 transdifferentiation once the lung is mature and likely denotes that cells from this cluster have terminally differentiated to other cell states. Trajectory inference, AT2/AT1 co-staining and lineage tracing would be required to definitively understand the fate of different clusters of AT2 cells with lung maturation and hyperoxic lung injury.

Alas2 utilizes a co-factor pyridoxal phosphate to catalyze the first step of heme synthesis in mitochondria. Heme plays important roles in mitochondrial oxygen phosphorylation, as it serves a prosthetic group for several proteins that constitute the complexes of mitochondrial electron transport chain. A previous study has shown that lung Alas2 gene expression is reduced in a baboon model of BPD [69]. This is corroborated by the present study showing that neonatal hyperoxia decreased Alas2 protein abundance in some AT2 cells. Reduced Alas2 abundance may lead to hyperoxia-induced reduction in mitochondrial respiration in lung epithelial cells [70]. This may drive cell metabolism towards glycolysis and abnormal cell proliferation in this cellular compartment [70].

Future studies using KEGG pathway enrichment analysis of genes in these clusters could reveal mechanistic insights into AT2 cells’ metabolism and function, including transdifferentiation and proliferation, during lung maturation and hyperoxic lung injury. Our model mimicked the clinical scenario where infants are briefly exposed to hyperoxia and weaned rapidly to room air. In other instances, lower concentrations of oxygen may be used to obviate hyperoxic injury in neonates. Whether lower concentrations of hyperoxic exposure also cause persistent changes in the lung cellular landscape and whether this occurs in a concentration-dependent manner remain to be investigated. Male sex is associated with increased severity of BPD [71]. Further investigation is warranted to explore sexual dimorphism of the transcriptome of lung cells at a single cell level during hyperoxic lung injury.

In conclusion, using scRNA-seq, we report a dataset describing the transcriptome of lung cells from the saccular stage to the mature lung, after exposure of mice to neonatal hyperoxia. We identified markers for all major lung cell types including multiple populations of mesenchymal, endothelial, epithelial and immune cells. We identified five subpopulations of AT2 cells that were dynamically altered with postnatal development and/or by neonatal hyperoxia. These subpopulations have distinct characteristics and likely different functions in the lung milieu (Fig. 9B). Persistent alterations in these specific cell states were identified into adulthood predicting lasting abnormalities in lung architecture and function. This may predispose to worsened injury upon a second insult in later life, or to persistently aberrant lung cell proliferation. This study reveals neonatal programing of adult lung dysfunction.

Author contributions

Conception and Design: AS, HY, PAD; Data acquisition and analysis: AS, HY, NO, JW, ALP, AR, MED; Data Interpretation: HY, PAD, AS, JW, GW, NO; Drafting of the manuscript: HY; Revision of the manuscript: HY, PAD, AS, JW, NO. GW provided technical support for single cell encapsulation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.