Abstract 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% O[2]) 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. Keywords: Bronchopulmonary dysplasia, Hyperoxic lung injury, Progenitor cells, Alveolar type I and Alveolar type II cells, Lipid and matrix homeostasis Graphical abstract Image 1 [41]Open in a new tab Highlights * • Using scRNA-seq, a total of 45 clusters within 32 cell states are identified in mouse lung. * • Neonatal hyperoxia results in persistent alterations in the murine lung cells into adulthood. * • Neonatal hyperoxia causes persistent changes in each of the five Type II cell states with unique signatures. * • These alterations are also observed in premature infants requiring short- or long-term mechanical ventilation. * • Neonatal hyperoxia causes alveolar simplification in adult mice. 1. 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) [[42]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 [[43][2], [44][3], [45][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 [[46][5], [47][6], [48][7], [49][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 [[50]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 [[51]10,[52]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 [[53]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) [[54]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 [[55][14], [56][15], [57][16], [58][17], [59][18], [60][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 [[61][20], [62][21], [63][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 [[64]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 [[65]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 [[66]25]. In mice, lung alveolarization starts at pnd4. Hyperoxic exposure in newborn mice causes persistent lung injury into adulthood [[67]26]. This model is frequently utilized to investigate pathogenesis and to identify potential therapeutic targets for hyperoxic lung injury in BPD [[68]27,[69]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. 2. Methods 2.1. 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% O[2]) for 72 h in an A-chamber (BioSpherix, Parish, NY) [[70]29,[71]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). 2.2. 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 [[72]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. 2.3. scRNA-seq Single cell encapsulation was performed using the Chromium Single Cell Chip G kit on the 10 × Genomics Chromium Controller [[73]32,[74]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 [[75]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 [[76]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. 2.4. 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 [[77]34,[78]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 O[2]) (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. 2.5. 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 [[79]4]. Utilization of human lung samples was done in compliance with the Institutional Review Board guidelines of Women and Infants Hospital, Providence RI. 2.6. 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. 2.7. 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. 2.8. 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 [[80]29]. In brief, we inflated non-lavaged mouse lungs with 1% low melt agarose at a pressure of 25 cm H[2]O, 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. 2.9. 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. 3. Results 3.1. 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 ([81]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 [[82]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, O[2]/pnd7, Air/pnd60, and O[2]/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, O[2]/pnd7, Air/pnd60, and O[2]/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 ([83]Fig. 2A, [84]Supplementary Table 1). The top marker genes in each cluster was shown in [85]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 ([86]Fig. 2B–C, [87]Table 1, [88]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 ([89]Fig. 2B, [90]Supplementary Table 2). These results confirm known lung cell heterogeneity but from an unbiased framework of gene expression. Fig. 1. [91]Fig. 1 [92]Open in a new tab 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. [93]Fig. 2 [94]Open in a new tab 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, O[2]/pnd7, Air/pnd60 and O[2]/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, O[2]/pnd7, Air/pnd60 and O[2]/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