Introduction

During the alveolar phase of lung development, alveolar walls protrude from secondary crests into the lumen of existing airways and form new alveoli to increase the lung’s surface area for gas exchange ¹. Many premature infants are born prior to the onset of alveologenesis and develop bronchopulmonary dysplasia (BPD), a chronic lung disease of prematurity caused by developmental arrest of secondary alveolar septation ². Over 50 years ago, Northway and colleagues described BPD as a heterogeneous pattern of severe lung injury with pathologic findings of severe inflammation, airway dysplasia, and fibrosis ³. Despite advances in neonatal care, the incidence of BPD remains unchanged as 30% of infants born before 30 weeks of gestation will develop BPD ^2,4. More infants are now surviving at earlier gestational ages, and the lungs of these infants have a “new” BPD phenotype characterized by a homogenous process of alveolar simplification, dysmorphic pulmonary vasculature, and mild airway thickening and fibrosis ^1,2. BPD is defined clinically by the need for prolonged respiratory support and supplemental oxygen, and infants with BPD suffer from long-term morbidity including severe respiratory infections, reactive airway disease, pulmonary hypertension, and neurodevelopmental impairment ^1,2.

The inflammatory response plays a significant role in BPD pathogenesis and results from a combination of perinatal infection, oxygen toxicity, and barotrauma from mechanical ventilation ¹. Several recent studies have suggested a role for the pro-fibrotic cytokine transforming growth factor-beta (TGFβ) in BPD pathogenesis ^4–9. Infants with BPD have increased TGFβ in serum and bronchoalveolar lavage fluid ^10–13. Experimentally, murine models of BPD show increased TGFβ signaling and activation of downstream fibrotic pathways ^{7–9,14–17}. At the same time, TGFβ plays an important role in normal development, and complete loss of TGFβ signaling impairs either embryonic lung development or alveolar septation ^18–29. TGFβ is a known regulator of myofibroblasts in fibrotic disease ³⁰, but it also drives mesenchymal cells to commit towards the myofibroblast lineage during embryonic lung development ²¹. Despite its essential role in normal lung development, the functional significance of increased TGFβ signaling in BPD pathogenesis remains unclear.

Alveolar myofibroblasts are contractile mesenchymal cells that are critical for alveolar septation ³¹. They are characterized by the expression of platelet-derived growth factor receptor alpha (PDGFRα), alpha-smooth muscle actin (α-SMA), and the production of elastin ^21,31–34. Myofibroblasts are derived from PDGFRα+ mesenchymal cells at birth and peak during alveologenesis when they colocalize with secondary septae ^{21,31,33–35}. Ablation of PDGFRα+ cells in neonatal mice causes alveolar simplification, while neonatal hyperoxia treatment causes loss of PDGFRα+ fibroblasts, dysregulated alveolar elastin deposition, and impaired PDGFRα+ cell contractility ^31,33,34,36. Inhibition of myofibroblast contraction has been shown to lead to alveolar simplification ³⁷. These studies support a fundamental functional role for myofibroblasts during alveologenesis, but whether the number of myofibroblasts is reduced in BPD or experimental models of BPD remains controversial ^{31,34,38–40}.

To identify conserved mechanism of injury in BPD pathogenesis, we implemented a novel strategy to compare two murine models of alveolar simplification by single-cell RNA-sequencing (scRNA-seq): neonatal hyperoxia exposure and loss of epithelial TGFβ signaling. Using flow cytometry and sequencing, we observed a dramatic reduction of myofibroblasts in both models of disease. Additionally, we found that increased TGFβ signaling, decreased PDGFRα signaling, and impaired proliferation are hallmark features of injured myofibroblasts during alveologenesis. Blocking TGFβ through several different approaches consistently worsened hyperoxia-induced disease, suggesting TGFβ plays an essential homeostatic function in both normal alveolar development and hyperoxia-induced disease, but that increased TGFβ signaling does not seem to be an important driver of alveolar simplification. We also demonstrated that PDGFRα+ cells undergo robust proliferation during normal alveologenesis, but this proliferation is impaired in both hyperoxia-treated and TGFβ-manipulated mice. We showed that modulation of PDGFRα+ cell proliferation is sufficient to cause alveolar simplification, even in the absence of hyperoxic injury. These studies demonstrate that PDGFRα+ cells are profoundly sensitive to injury during alveologenesis, they require TGFβ signaling, and their proliferation is essential for normal alveolar development.

Results

Neonatal Hyperoxia Exposure and Loss of Epithelial TGFβ Signaling Both Cause Alveolar Simplification

To identify conserved features of lung injury in bronchopulmonary dysplasia (BPD), we compared two different mouse models of lung injury known to cause alveolar simplification in mice. We first implemented the hyperoxia-induced model of BPD in which newborn mice are exposed to 75% oxygen from postnatal day (P) 0-10 and then recovered in room air until analysis (Figure 1A). This model takes advantage of the fact that mice are born during the saccular phase of lung development and undergo alveolarization from postnatal day P4 – P39 ^41–43. We harvested lungs to analyze alveolar architecture after completion of alveolarization at P40. Hyperoxia-treated mice showed significant alveolar airspace enlargement compared to controls (Figure 1B). Consistent with findings from other groups using similar protocols ^31,34, we did not see evidence of fibrosis or scarring. In parallel, we used the flexiVent rodent ventilator to measure respiratory mechanics. Hyperoxia-treated mice had increased total lung capacity, increased compliance, and decreased elastance compared to control mice (Figure 1C). These findings are consistent with those of emphysematous lungs, and these physiologic changes align with the enlarged, simplified alveolar structures observed by histology ^44–46.

Figure 1.

For our second model of alveolar simplification, we generated mice which lack the critical TGFβ receptor TGFBR2 in lung epithelium by crossing Nkx2.1-cre to the Tgfbr2 conditional allele (Tgfbr2^F/F;Nkx2.1-cre). These mice are viable at birth, have no respiratory distress, gain normal weight, and generate expected ratios of offspring with respect to each genotype (data not shown). Lungs from conditional knockout (cKO) mice show enlarged airspaces consistent with alveolar simplification, confirming earlier reports that TGFβ signaling in lung epithelium is required for normal secondary septation (Figure 1D-E) ²⁶. We then used the flexiVent system to analyze respiratory mechanics in cKO mice and observed increased compliance and decreased elastance comparable to what we saw in hyperoxia-exposed mice (Figure 1F). Previous work demonstrated that Tgfbr2^F/F;Nkx2.1-cre mice are protected from lung disease caused by transgenic overexpression of TGFβ1 or neonatal hyperoxia treatment ²⁶. To our surprise, Tgfbr2^F/F;Nkx2.1-cre mice were not protected from disease in our model of neonatal hyperoxia, but rather developed worse disease. These data suggest TGFβ signaling to the lung epithelium is required for both normal alveologenesis and that it plays a protective role in the response to hyperoxia-induced injury.

Profound Loss of PDGFRα+ Fibroblasts During Hyperoxia Treatment

To begin exploring the mechanisms underlying alveolar simplification in these models, we treated wildtype mice with hyperoxia from P0-P10 and analyzed dissociated lung cells by flow cytometry at P10. While we observed changes in many cell types of the lung with hyperoxia treatment, the most striking result was the profound loss of PDGFRα+ fibroblasts (Figure 2A-B). Using flow cytometry, we found that both the percentage and absolute number of PDGFRα+ fibroblasts were substantially reduced with hyperoxia treatment (Figure 2B). Next, we analyzed Tgfbr2^F/F;Nkx2.1-cre mice at P10 in normoxia and observed a similar phenotype of reduced PDGFRα+ fibroblasts (Figure 2C-D). While the magnitude of reduction was less impressive than we observed in the hyperoxia model, the most significant change in cell numbers in Tgfbr2^F/F;Nkx2.1-cre mice was the reduction in PDGFRα+ fibroblasts. Together, these results are consistent with earlier studies which showed a loss of PDGFRα+ cells with neonatal hyperoxia treatment and reinforce the critical role of PDGFRα+ fibroblasts in normal alveolar development ^31,36.

Figure 2.

scRNA-seq Reveals Loss of Myofibroblasts in Both Models of Alveolar Simplification

To identify conserved molecular and cellular mechanisms underlying pathologic alveolar simplification, we performed single-cell RNA sequencing on lungs from Tgfbr2^F/F (CTRL) and Tgfbr2^F/F;Nkx2.1-cre (cKO) mice in either normoxia or hyperoxia (Figure 3A). This approach allowed for simultaneous comparison of both injury models: normoxia versus hyperoxia in CTRL mice, and CTRL versus cKO mice in normoxia (Figure 3B). Given the reduction of PDGFRα+ fibroblasts observed by flow cytometry in both models of lung injury, we hypothesized that changes in the lung mesenchyme might play a substantial role in the pathogenesis of alveolar simplification. To ensure sufficient mesenchymal cells in our analysis, we used flow cytometric sorting to enrich for mesenchymal and epithelial cells by limiting the input of hematopoietic and endothelial cells. After processing and combining all samples, we identified 24 clusters from 26,610 cells (Figure S1A-C). We focused our analysis on the lung mesenchyme in which we used the Seurat software package, differential gene expression, and comparison to published scRNA-seq data to assign identities to seven distinct mesenchymal clusters. (Figures 3C and S2) ^47–50.

Figure 3.

Although flow cytometry had already shown a profound loss of PDGFRα+ cells in these two injury models, we could not attribute this loss to a specific fibroblast subset. Our single cell data enabled higher resolution assessment of the dynamics of fibroblast subsets (Figure 3D). In particular, we observed a significant loss of alveolar and ductal myofibroblast clusters in both models of injury (Figure 3D-E), suggesting that the reduction of PDGFRα+ cells by flow cytometry was due to the loss of alveolar and ductal myofibroblasts. To externally validate these findings, we re-analyzed two recently published scRNA-seq studies in which mice were treated with 85% oxygen from P0-P14 ^47,50. Consistent with our own data, hyperoxia-treated mice in both studies showed a significant loss of myofibroblasts within the lung mesenchyme (Figure S3).

To gain insight into the molecular mechanisms driving the loss of myofibroblasts, we used Qiagen Ingenuity Pathway Analysis (IPA) to compare differential gene expression in the myofibroblasts in both models of injury ⁵¹. After filtering for upregulated pathways, we identified 32 predicted upstream regulators that were shared in both models (Figure 3F). Of these pathways, we noted several predictions corresponding to TGFβ signaling, inhibitors of cell cycle, and Wnt signaling. These results are notable because TGFβ is a known driver of myofibroblast differentiation and profibrotic programs in fibroblasts ^21,30, while its overall contribution to alveolar simplification and BPD is unclear ^5,52. The prediction of increased cell cycle inhibition is also of interest because it suggests that the reduction of myofibroblasts in these models might be caused by their impaired proliferation.

We used the NicheNet software package to gain insight into the specific cell-cell interactions which might be impacting the myofibroblasts in these two models of lung injury ⁵³. We reasoned that epithelial-myofibroblast signals would be particularly important, since the defect in the cKO model is restricted to epithelial cells, yet the most prominent phenotype is found in myofibroblasts. To explore this further, we used NicheNet to identify shared patterns of changes in ligand-receptor signaling originating from the lung epithelium to the myofibroblasts (Figure 4A). Here we again noted an increase in TGFβ signaling to the myofibroblasts (Figure 4B). Of the decreased pathways, the most notable were Pdgfa-Pdgfra and Shh-Hhip, both of which are known to be critical for myofibroblast differentiation and function (Figure 4C) ^29,54. Additionally, the predicted reduction of Pdgfra signaling aligns with our flow cytometry data showing a loss of PDGFRα+ fibroblasts in both injury models as well as a reduction in the mean fluorescence intensity of PDGFRα antibody staining in these cells (Figure 2B-D). Together, these data suggest that epithelial dysfunction may cause myofibroblast defects through loss of supportive signals and gain of inhibitory signals.

Figure 4.

TGFβ Signaling Plays Homeostatic Role in Normal Alveolar Development and in Hyperoxia

Both IPA and NicheNet analyses predicted that TGFβ signaling is activated in myofibroblasts in both models of alveolar simplification (Figures 3 and 4). Given the conflicting literature regarding the role of TGFβ in normal lung development ^5,18,52, we sought to determine the functional significance of TGFβ signaling during postnatal alveologenesis in both normoxic and hyperoxic conditions. We hypothesized that if excessive TGFβ signaling is a pathologic response in hyperoxia, titration of a pan-TGFβ-blocking antibody (1D11) could identify a therapeutic window to protect from disease while permitting normal development under normoxic conditions. Alternatively, worsened disease with 1D11 treatment would suggest that increased TGFβ is a compensatory response to injury rather than a primary driver of disease in the hyperoxia model. To test these hypotheses, we injected wildtype mice with 0, 10, 20 or 30 mg/kg of 1D11 from P2-P10 in either normoxia or hyperoxia and then harvested at P40 for analysis by histology (Figure 5A). To our surprise, we found that while only the highest dose of 1D11 treatment caused alveolar simplification in normoxia, all doses of 1D11 treatment caused worse alveolar simplification in hyperoxia compared to PBS (Figure 5B). Consistent with these histology results, lung physiology studies of 1D11-treated mice showed a pattern of emphysematous changes with increased compliance and decreased elastance (Figure 5C). These changes are similar to what we previously observed in hyperoxia-treated mice and in Tgfbr2^F/F;Nkx2.1-cre mice and complement the histologic findings of alveolar simplification in 1D11-treated mice. Together, these results support the alternative hypothesis that TGFβ signaling plays a homeostatic role in both normal alveolar development and in hyperoxia.

Figure 5.

Because 1D11 antibody treatment is systemic, we were unable to attribute the results of these studies to a specific cell type within the lung. Our current study, consistent with earlier reports, demonstrates that TGFβ signaling to the epithelium is required for normal alveolar development as well as an adaptive response to hyperoxia ²⁶. To interrogate the role of TGFβ signaling to the lung mesenchyme, we sought to cross Tgfbr2^F/F mice to either the Pdgfra-CreERT2 or Gli1-CreERT2 allele to target Tgfbr2 in myofibroblasts and mesenchymal cells during alveologenesis. While Pdgfra-CreERT2 has been used to specifically target PDGFRα+ cells ⁵⁵, many studies have also used the broader Gli1-CreERT2 to target myofibroblasts during alveolar lung development ^29,54. Because the Pdgfra-CreERT2 and Gli1-CreERT2 alleles were generated using a knock-in/knock-out approach, mice carrying the CreERT2 allele are haploinsufficient for Pdgfra or Gli1, respectively ^55,56. We conducted hyperoxia studies in each of these two lines and found that the Gli1-CreERT2 allele alone does not disrupt alveolar development in either normoxia or with hyperoxia treatment (Figure S4A-B). In contrast, mice with the Pdgfra-CreERT2 allele undergo normal development in normoxia, but they develop worse alveolar simplification with hyperoxia treatment compared to their cre-negative littermates (Figure S4A-B). These data support use of the Gli1-CreERT2 allele to target the lung mesenchyme in either normoxia or hyperoxia while limiting the use of Pdgfra-CreERT2 mice (and other Pdgfra-haploinsufficient alleles) to normoxic conditions.

Based on these validation studies, we generated Tgfbr2^F/F;Gli1-CreERT2 mice to study the effects of TGFβ signaling to the lung mesenchyme. Interestingly, we observed a similar pattern of lung disease as 1D11 treatment when deleting Tgfbr2 in the lung mesenchyme: Tgfbr2^F/F;Gli1-CreERT2 mice developed alveolar simplification in normoxia when compared to littermate controls, and they developed worse disease than their littermates with hyperoxia exposure (Figure S5A-B). Together, these results show that TGFβ signaling is required to both the lung epithelium and mesenchyme for normal development, that it plays a homeostatic role in response to neonatal hyperoxia treatment, and that all our efforts to inhibit TGFb made hyperoxia-induced alveolar simplification worse.

To address mechanisms of TGFβ activation in normal alveolar development and in the hyperoxia model of alveolar simplification, we generated Itgb6^F/F mice and crossed them to the Nkx2.1-cre allele. The αvβ6 integrin is expressed predominantly in epithelial tissues, activates TGFβ in vivo, and has been shown to be a critical mediator of lung injury in models of acute lung injury and pulmonary fibrosis ^57–62. Itgb6^F/F;Nkx2.1-cre mice lack αvβ6 on lung epithelial cells and undergo normal alveolar development, but to our surprise, show no difference in disease severity with hyperoxia treatment (Figure S5C-D). Within the mesenchyme, αvβ1 and αvβ8 have been shown to activate TGFβ in vivo ^62–65, so we next deleted all αv-integrins in the lung mesenchyme by crossing Itgav^F/F mice to the Gli1-CreERT2 allele. These mice develop spontaneous alveolar simplification in normoxic conditions, have lung physiology parameters consistent with an emphysematous phenotype, and develop worse lung disease with hyperoxia exposure (Figure 6A-C). These data suggest that epithelial αvβ6 does not play a role in TGFβ activation during normal development or neonatal hyperoxia, while αv-integrins in the lung mesenchyme are required for normal development and play a protective role in response to hyperoxia.

Figure. 6.

Impaired Proliferation of PDGFRα+ Fibroblasts with Hyperoxia Treatment

Flow cytometry and scRNA-seq both showed a reduction of PDGFRα+ fibroblasts in hyperoxia-treated mice, which suggests these cells are exquisitely sensitive to injury. To address whether this might be due to impaired proliferation as suggested by IPA analysis, we exposed wildtype mice to hyperoxia and treated them with the nucleoside analogue 5-ethynyl-2’-deoxyuridine (EdU) to identify proliferating cells. We found a decrease of PDGFRα+ cells in hyperoxia-treated mice by P8 both by percentage and by absolute number (Figure 7A-C), which persisted to P10 and continued even with 4 days of normoxic recovery. Using EdU analysis, we found hyperoxia caused a proliferative defect in PDGFRα+ cells as early as P4, preceding the decrease in cell number. This early impairment of proliferation by PDGFRα+ cells was the most striking change observed across all subsets in our EdU analysis, and suggests that a proliferative defect contributes to the decreased cell number.

Figure 7.

After observing impaired proliferation of PDGFRα+ cells in the hyperoxia model of lung injury, we explored whether this mechanism is conserved across our other models of alveolar simplification. We used EdU treatment and flow cytometry to characterize cell subsets and proliferation in Tgfbr2^F/F;Nkx2.1-cre mice and 1D11-treated wildtype mice. Remarkably, both models showed a similar reduction of PDGFRα+ cells along with decreased PDGFRα+ cell proliferation as quantified by EdU uptake (Figure 8A-D). These results provide functional validation of our IPA analyses showing enrichment of cell cycle inhibitory pathways in injured myofibroblasts, and they suggest impaired PDGFRα+ cell proliferation might be a conserved feature of developmental lung injuries which result in alveolar simplification.

Figure 8.

Impaired Proliferation of PDGFRα+ Fibroblasts Is Sufficient to Cause Alveolar Simplification

After observing impaired PDGFRα+ cell proliferation in multiple models of alveolar simplification, we sought to determine whether inhibiting proliferation in these cells would be sufficient to cause disease. We generated mice with conditional deletion of Ect2, a protein required for cytokinesis ⁶⁶, by crossing the Ect2^F/F mice with the Pdgfra-CreERT2 allele. We found earlier that PDGFRα+ cells underwent robust proliferation through P10 in normoxic conditions (Figure 7C), so we analyzed Ect2^F/F;Pdgfra-CreERT2 mice at P14 to characterize the impact of this mouse model on the cellular composition of the lung during early alveologenesis. Using flow cytometry, we observed a significant reduction in PDGFRα+ cells within the lung by percentage and absolute number (Figure 9A-B). While this reduction in PDGFRα+ cells was expected in Ect2^F/F;Pdgfra-CreERT2 mice, it was interesting that Pdgfra expression was reduced in the remaining PDGFRα+ cells as quantified by mean fluorescence intensity (MFI). Earlier studies found a correlation between Pdgfra expression levels and fibroblast proliferation ⁶⁷, so we wondered whether this phenotype of reduced Pdgfra expression was conserved across the other injury models in this study which showed reduced PDGFRα+ cell proliferation. Indeed, quantification of PDGFRα MFI by flow cytometry confirmed that Pdgfra expression is significantly reduced in the lung mesenchyme and specifically in PDGFRα+ cells in each of the injury models we analyzed in this study (Figure S6A).

Figure 9.

Next, we aged Ect2^F/F;Pdgfra-CreERT2 mice until P40 to analyze their lungs by either morphometry or lung physiology. By histology, we observed enlarged alveolar airspaces and increased MLI’s in conditional knockouts (Ect2^F/F;Pdgfra-CreERT2) compared to controls (Figure 7C). Importantly, conditional heterozygous mice showed no difference in MLI compared to cre-negative mice, which confirms that haploinsufficiency of Pdgfra caused by using the Pdgfra^CreER allele in normoxic conditions is not responsible for this phenotype (Figures 7C and S6A-B). By lung physiology, Ect2^F/F;Pdgfra-CreERT2 mice produced the same emphysematous phenotype of increased compliance and decreased elastance which we observed in several other models of alveolar simplification in this study (Figure 7D). Together, these data demonstrate that impaired proliferation of PDGFRα+ cells is sufficient to cause alveolar simplification, even in the absence of environmental insult or modulation of TGFβ signaling.

Discussion

Using a combination of flow cytometry, lung physiology, and scRNA-seq, we identified several core features of lung injury that are conserved across multiple mouse models of alveolar simplification. We first compared neonatal hyperoxia exposure (75% hyperoxia P0-P10) to the loss of epithelial TGFβ signaling (Tgfbr2^F/F;Nkx2.1-cre). Flow cytometric analysis of both models showed a significant reduction of PDGFRα+ cells, while scRNA-seq studies attributed these changes to a loss of alveolar and ductal myofibroblasts. Qiagen IPA and NicheNet ligand-receptor analyses identified several pathways that were upregulated in injured myofibroblasts in both models: increased TGFβ signaling, enrichment of inhibitors of cell cycle, and decreased Pdgfa-Pdgfra signaling. Using a combination of pharmacologic and genetic approaches, we demonstrated that increased TGFβ signaling in the mesenchyme does not seem to be a driver of alveolar simplification, but rather TGFβ signaling is critical for both normal alveolar development and for protecting against impaired alveolar development in response to hyperoxia. In contrast, we show that PDGFRα+ cell proliferation is reduced in multiple models of alveolar simplification, and that this impaired proliferation on its own is sufficient to cause alveolar simplification.

Several recent studies have identified lung mesenchymal cells as critical mediators of alveologenesis. Li et al. established the requirement of myofibroblasts during postnatal development by specifically ablating PDGFRα+ cells at the onset of alveologenesis and showing that these mice developed alveolar simplification ³¹. Other work has identified essential functions of myofibroblasts during alveolar development, including contractility ³⁷ and mitochondrial energetics ⁶⁸. Ricetti et al used bulk RNA sequencing to show a skewing from myofibroblast to matrix fibroblast phenotypes after hyperoxia exposure along with reduced Ki-67 uptake and contractility in myofibroblasts ³⁶. Our current study builds on this existing body of work by characterizing significant changes in PDGFRα+ cell number in several models of developmental lung injury. By flow cytometry, we found reduction in the number of PDGFRα+ cells in hyperoxia-treated mice and in Tgfbr2^F/F;Nkx2.1-cre mice during the early phase of alveologenesis. Using scRNA-seq, we attributed this reduction to the alveolar and ductal myofibroblast populations. Importantly, we identified myofibroblast proliferation as a critical feature of normal alveolar development and validated this observation in several experimental models.

NicheNet analysis of our scRNA-seq data predicted decreased Pdgfa-Pdgfra signaling from the epithelium to the myofibroblasts in both hyperoxia-treated and Tgfbr2^F/F;Nkx2.1-cre mice. This observation is interesting because several studies have focused on the requirement of Pdgfa-Pdgfra interactions for myofibroblast identity and function. Deletion of Pdgfra in mesenchymal populations of the neonatal lung has been shown to disrupt alveologenesis ^54,69. Other studies have found lower Pdgfra expression in lung samples from BPD patients ^40,70. In our current study we use flow cytometry to show a loss of PDGFRα+ cells in multiple models of alveolar simplification: hyperoxia-treated mice, Tgfbr2^F/F;Nkx2.1-cre mice, 1D11-treated mice, and Ect2^F/F;Pdgfra-CreERT2 mice. In addition to the loss of PDGFRα+ cells, we also observed decreased Pdgfra expression within the remaining PDGFRα+ cells (Figure S6A). Given that EdU uptake is also reduced amongst the remaining PDGFRα+ cells, these results suggest Pdgfa-Pdgfra signaling might regulate myofibroblast proliferation in addition to their identity and function during neonatal lung development.

As part of our current studies, we conducted control experiments with the Gli1-CreERT2 and Pdgfra-CreERT2 alleles and discovered results with broad implications for others utilizing similar genetic tools for their work. Both alleles were generated using a knock-in/knock-out gene targeting approach which disrupts the native expression of Gli1 and Pdgfra, respectively ^55,56. Because Hedgehog and Pdgfa-Pdgfra signaling pathways are important for lung development ^71,72, we sought to validate the effects of haploinsufficiency of Gli1 and Pdgfra when using these alleles. Both Gli1-CreERT2 and Pdgfra-CreERT2 mice underwent normal alveolar development in normoxic conditions when compared to cre-negative littermates. In contrast, Pdgfra-CreERT2 mice developed worse alveolar simplification after hyperoxia treatment compared to their cre-negative littermates. The Gli1-CreERT2 results suggest this allele can be used for conditional knockout studies in either normoxic or hyperoxic conditions with little impact from Gli1 haploinsufficiency. The Pdgfra-CreERT2 results, however, should serve as a caution to researchers utilizing genetic tools with Pdgfra-haploinsufficiency.

Recent work has established an essential role for TGFβ signaling in normal alveolar development ⁵², but the functional role of TGFβ in models of BPD and alveolar simplification remains unclear. Our scRNA-seq data showed increased TGFβ signaling in myofibroblasts in two injury models, so we hypothesized that excess TGFβ signaling might be a driver of disease pathogenesis. To determine the functional significance of TGFβ in the hyperoxia model, we treated mice with the pan-TGFβ-blocking antibody 1D11. Because TGFβ is required for normal alveologenesis, we hypothesized a dose titration approach might identify a “sweet spot” to neutralize an excess of TGFβ in hyperoxia while preserving normal development in normoxia. Instead, we found that every dose of 1D11 treatment worsened hyperoxia-induced alveolar simplification while only the highest dose caused disease in normoxia. While the 1D11-normoxia results confirm an essential role for TGFβ in normal alveolar development, the 1D11-hyperoxia results suggest TGFβ plays a homeostatic function rather than a pathologic role in the hyperoxia model of injury. Because 1D11 treatment neutralizes TGFβ ligands systemically, we used Nkx2.1-cre and Gli1-CreERT2 alleles to conditionally delete Tgfbr2 in the lung epithelium and mesenchyme, respectively. Confirming earlier reports, we found that disrupting TGFβ signaling to either of these populations causes alveolar simplification ^26,29. Of significance, however, both Tgfbr2^F/F;Nkx2.1-cre and Tgfbr2^F/F;Gli1-CreERT2 mice developed worse disease with neonatal hyperoxia treatment. Taken together with our 1D11-hyperoxia experiments, we conclude that TGFβ signaling to both the lung epithelium and mesenchyme is required for normal alveolar development and is protective rather than pathologic in hyperoxia-perturbed alveolar development.

Little is known about mechanisms of TGFβ activation in alveolar development and neonatal lung injury. TGFβ ligands are secreted as latent, inactive complexes that are anchored to the cell surface or the extracellular matrix ⁶². These complexes prevent TGFβ from engaging its receptors, and therefore TGFβ activation is tightly regulated in tissues despite the high levels of latent TGFβ complex⁶². Integrins are cell-surface proteins that regulate interactions in the extracellular matrix, and several studies by our group have demonstrated that the epithelial integrin αvβ6 is critical for TGFβ activation in models of acute lung injury and pulmonary fibrosis ^58,60,73,74. However, TGFβ activation by epithelial αvβ6 does not play a role in alveolar development in either normoxia or hyperoxia. We then deleted αv-integrins in the lung mesenchyme since αvβ1 and αvβ8 are the only other integrins known to activate TGFβ in vivo ⁶². Interestingly, deletion of αv-integrins under normoxic conditions only produced a modest degree of alveolar simplification but clearly led to more severe alveolar simplification after hyperoxia.

We originally hypothesized that increased TGFβ signaling was a pathologic response to neonatal hyperoxia treatment and that targeting the TGFβ pathway might protect from disease pathogenesis. While our results disproved this hypothesis, they also seem at odds with published data where blocking TGFβ signaling was beneficial in various lung injury models. For example, it was unexpected that Tgfbr2^F/F;Nkx2.1-cre mice develop worse lung injury in our hyperoxia-induced model of BPD because earlier work showed these same mice to be protected from hyperoxia ²⁶. We suspect this difference in outcome can be explained by the more aggressive models of lung injury used by this earlier study. Specifically, these authors used either transgenic overexpression of TGFβ or 100% hyperoxia exposure which both resulted in neonatal mortality and severe inflammation, neither of which are observed in our model of 75% hyperoxia. Additionally, many earlier studies focused on TGFβ activation in hyperoxia-induced BPD utilized chronic 85% hyperoxia such as P0-P14 or P0-P28, both of which result in inflammation and fibrosis ^8,14–16,75. Similarly, one earlier study used transgenic expression of Il-1b in the developing airways to produce a BPD-like phenotype due to inflammation ⁷⁶. Deletion of Itgb6, and therefore loss of epithelial αvβ6 integrins, conferred partial protection from disease in this injury model⁶¹. A common theme across these earlier studies implicating TGFβ in pathology is that they utilized aggressive models of lung injury which caused severe inflammation, fibrosis, and neonatal mortality, all of which recapitulate the histopathologic features of “old BPD” in the clinical setting. In contrast, the hyperoxia model used in our current study more closely represents the “new BPD” phenotype of emphysematous changes with simplified alveolar structures in the absence of lung scarring and fibrosis ^1,2. Collectively, these data suggests that while inhibition of TGFβ may have some protective role in severe injury models associated with significant fibrosis, the major functions of TGFβ signaling during normal alveolar development and in response to moderate hyperoxia are homeostatic and protective.

There are several limitations to our current study. We used flow cytometry to phenotype population-level changes across several models of neonatal lung injury. While flow cytometry is a powerful tool to quantify population subsets and sort cells for scRNA-seq, it requires tissue digestion and generation of single-cell suspensions for downstream analysis. The cells which survive these protocols may represent a skewed population and may not be completely representative of in vivo conditions. For example, many fragile and dying/apoptotic cells do not survive this processing and would be excluded from analysis. Another drawback to flow cytometry is that our ability to define population subsets is limited by the knowledge and feasibility of using cell-surface markers to faithfully distinguish cell types within each population. Recent work by our group has identified cell-surface markers to characterize fibroblast subsets in the adult lung by flow cytometry, but we are unaware of similar protocols to address fibroblast heterogeneity in the neonatal lung ⁷⁷. We used scRNA-seq to characterize cell subsets within the lung mesenchyme, but the identification and validation of cell-surface proteins suitable for flow cytometric analysis of these same populations is beyond the scope of this current study. Our scRNA-seq data, together with emerging scRNA-seq from multiple other laboratories should provide a rich dataset for others to identify cell surface markers that allow more precise analysis and sorting of distinct mesenchymal populations in the developing lung.

Our current work highlights the essential role of myofibroblast proliferation and TGFβ signaling in normal alveologenesis. We show that myofibroblasts are depleted during neonatal lung injury and that their loss is at least partially due to impaired proliferation and expansion during this critical window of development. We confirm the results of others who similarly observed increased TGFβ signatures in the hyperoxia-injured neonatal lung and have generated interest in targeting this pathway as a therapeutic intervention for BPD. However, after conducting exhaustive studies targeting TGFβ ligands, receptors and activating integrins, we conclude that increased TGFβ signaling in myofibroblasts more likely represents a failed compensatory mechanism rather than a central driver of disease pathogenesis. Clinical BPD remains a heterogeneous disease encompassing both the severe inflammation, scarring and fibrosis of the “old BPD” phenotype as well as the increasingly prevalent “new BPD” phenotype of alveolar simplification and emphysematous changes ^1,2. Therefore, further work is necessary to contextualize our results with respect to these different clinical BPD phenotypes.

In summary, our results underscore the importance of impaired myofibroblast proliferation as a central feature of alveolar simplification in several models of murine lung injury. Further work is needed to validate these findings in human BPD. If validated, these results would suggest that efforts to prevent or reverse this process could have therapeutic value in BPD. Our sequencing data should provide useful insights to the broader community studying alveolar development and neonatal lung injury. While our current work focuses on a few of the genes and pathways highlighted by these data, we are optimistic that others will utilize this dataset to expand our understanding of the molecular mechanisms driving both normal and aberrant alveolar development.

Materials and Methods

Mice

C57BL/6 (Stock No. 000664), Nkx2.1-cre (Stock No. 008661), Pdgfra-CreERT2 (Stock No. 032770), and Gli1-CreERT2 (Stock No. 007913), and Itgav^F/F (Stock No. 032297) mouse lines were obtained from the Jackson Laboratory. Tgfbr2^F/F mice (exon 2 conditional allele) were described previously ⁷⁸. Ect2^F/F were described previously ⁶⁶, and were obtained from Dr. Alan Fields at Mayo Clinic. All lines were maintained on the C57BL/6 genetic background except for the Ect2^F/F line, which was originally generated in BALB/c and backcrossed 5 generations to C57BL/6 for these studies. All animal experiments were in accordance with protocols approved by the Institutional Animal Care and Use Committee and Laboratory Animal Resource Center.

Neonatal Hyperoxia Treatment

The hyperoxia animal chamber (BioSpherix) was attached to a medical oxygen source controlled by a ProOx single gas controller (BioSpherix) set to maintain 75% oxygen under normobaric conditions. Birth was defined as <12 hours of life, and pups assigned to hyperoxia were transferred into the chamber with lactating dams and maintained from P0-P10 and then recovered in normoxia. During hyperoxia treatment, lactating dams were rotated between hyperoxia and normoxia to prevent maternal injury and to control for nutrition amongst both cohorts. Mice in both conditions were given nestlets and trail mix for additional enrichment. Mice in both conditions remained under typical 7a-7p light cycling, and chamber was checked daily to monitor temperature, humidity, and gas controller function.

Tamoxifen, Antibody, and EdU Treatments

Tamoxifen (Sigma) was dissolved in corn oil (Sigma) at 15 mg/ml, and 150 μg (10 μl) was administered via intraperitoneal (i.p.) injections on P2 and P4. Antibody clone 1D11 was used for pan-TGFβ-blocking studies ⁷⁹, and was generated in our laboratory from a hybridoma obtained from ATCC. Antibody was diluted in PBS and administered via i.p. injections at 0, 10, 20 or 30 mg/kg on P2, P4, P6, P8 and P10. For proliferation studies, EdU (Thermo Fisher) was reconstituted in DMSO at 100 mg/ml, diluted in PBS to 5 mg/ml, and injected i.p. at 75 mg/kg 24 hours prior to harvest.

Generation of Itgb6 Flox Mice

Itgb6 flox mice (Itgb6^F/F) were made by CRISPR/Cas9-aided homology-directed repair. Loxp sequences were inserted to flank exon 4 of Itgb6. Two guide RNA target sequences were chosen in the introns upstream and downstream of exon 4. crRNAs were obtained from IDT with input sequences CAGCTTATCATCCATCTAAA (upstream) and ACCTTCCTCTGACGCACTTT (downstream). Two 200bp donor DNAs were obtained from IDT. EcoRV sites (gatatc) were inserted following loxp sequences for screening purposes. The sequences of the donor DNAs are as follows:

Upstream donor DNA: TAATCTCTCCTTTATTTGGCTCACCTTTTCTGCAACCACACACCAAGAAAGGGCAGCTTATCATCCA TCTAAAATAACTTCGTATAGCATACATTATACGAAGTTATgatatcTGGATGCTACTTCTCCCTAGGAG ATATAAAATATCCCAACATACACCTCCTTCTGTCCTTCAATCCTCAC

Downstream donor DNA: TAACCTACATTTTTTTCTCTGAGTTTTTCTATCAAAATAACAATTTTTGCACCTTCCTCTGACGCACT TTATAACTTCGTATAGCATACATTATACGAAGTTATgatatcGGGAAATGTGGCTTTCACTCATTGCTG AGAGCAGCAGCCTTCATTGCAATTAAAGTCAAGAGGAAATGGG.

CRISPR/Cas9 complex (Cas9, crRNA, trRNA) and donor DNA were injected into C57BL/6 fertilized zygotes, which were then implanted into the oviducts of pseudopregnant female mice. 28 pups were born, and 6 of them had at least one allele with desired loxp insertion with two of them being homozygous for recombinant alleles. We picked one founder to expand the colony. Genotyping was performed with a forward primer 5’- CTGCAACCACACACCAAGAA-3’ and a reverse primer 5’- GCGTGACCTTATTAAGCTGCA-3’, which provide 196bp bands for wild type alleles and 236bp bands for flox alleles.

Histology

For morphometry studies, lungs were inflated with ice-cold 4% paraformaldehyde (PFA) under constant pressure of 25 cm H₂O for 5 minutes. Lungs were carefully dissected of attached structures and transferred into vial containing cold 4% PFA. Samples were rocked in 4°C overnight, washed 3x with PBS, and dehydrated in a series of ethanol (30%, 50% and 70%). Tissues were submitted in 70% ethanol to the UCSF Gladstone Histology and Light Microscopy core for further processing for paraffin embedding and tissue blocks were sectioned for Hematoxylin and Eosin (H&E) staining. For quantification of mean linear intercept (MLI), 6-8 sections of 5 μm thickness were sampled at 20 μm levels through each set of lungs. Images were acquired on a Nikon Ti Inverted Microscope using a 10x objective and DS-Ri2 color camera. For each sample, approximately 30 images were acquired from which 12 were randomly selected for MLI quantification using ImageJ software as described previously ⁸⁰. All samples were blinded during imaging and MLI quantification.

Lung Physiology

Pulmonary compliance and elastance were analyzed using the flexiVent system (SCIREQ) as previously described ⁴⁵. Mice were anesthetized with ketamine (100 mg/kg), xylazine (10 mg/kg), and acepromazine (3 mg/kg) before a tracheostomy was performed to cannulate the trachea with a 20-guage catheter. Mice were then paralyzed with pancuronium (0.1 mg/kg) and analyzed using the flexiVent rodent ventilator. All mice were analyzed in a blinded fashion.

Tissue Dissociation

Mouse lungs were harvested after perfusion through the right ventricle with PBS. Lungs were dissected into individual lobes, minced with razor blades, and suspended in protease solution [0.25 % Collagenase A (Millipore Sigma), 1 U/ml Dispase II (Millipore Sigma), 2000 U/ml Dnase I (Millipore Sigma) in Dulbecco’s Modified Eagle Medium (Thermo Fisher) containing 10 mM HEPES (Millipore Sigma) and 2% FBS (Millipore Sigma). The suspension was incubated in a 37°C water bath in a 15 ml conical for 25 minutes with aggressive trituration by glass Pasteur pipette every 8 minutes. The digestion was then quenched for 5 minutes in ice-cold PBS containing 5% FBS and 2 mM EDTA before passing cells through a 70 μm cell strainer. Cells were pelleted and resuspended in Red Blood Cell Lysing Buffer HybriMax (Sigma) for 5 minutes. RBC lysis was quenched dropwise with 5% FBS/PBS and cells were passed through a 40 μm cell strainer, washed, and resuspended in 5% FBS/PBS containing antibodies for FACS purification.

Flow Cytometry and Fluorescence-activated Cell Sorting (FACS)

After tissue dissociation, cells were resuspended in 5% FBS/PBS with Fc-blocking antibody (TruStain FcX; BioLegend) at 0.01 mg/ml for 5 minutes at room temperature. All antibodies were resuspended in Brilliant Stain Buffer (BD Biosciences) and added to Fc-blocked cells for final antibody concentration of 1:200. Cells were stained for 30 minutes on ice for surface antibody staining, washed, and then passed through another 40 μm before analysis and/or cell sorting. DRAQ7 (BioLegend) was used 1:1000 to identify dead cells. Flow cytometric cell counting was performed using CountBright Plus Absolute Counting Beads (Thermo Fisher). The following antibodies were used in this study: anti-CD9 (clone MZ3, PE; BioLegend), anti-CD31 (clone 390, A488, BV605; BioLegend, BD Biosciences), anti-CD45 (clone 30F-11, BV786, APC/Cy7, PE/Cy7; BioLegend, BD Biosciences), anti-Epcam (clone G8.8, FITC, PE, PE/Cy7; BioLegend), anti-I-A/I-E (clone M5/114.15.2, APC/Cy7, Spark UV 387; BioLegend), anti-Mcam (clone ME-9F1, A488, BUV496; BioLegend, BD Biosciences), anti-Pdgfra (clone APA5, APC, BV421, PE/Cy7; BioLegend). For proliferation studies, cells were stained for surface markers as described above and then washed, fixed, permeabilized, and processed for Click-iT EdU detection according to manufacturer’s protocol (Click-iT Plus EdU Alexa Fluor 647 Flow Cytometry Assay Kit; Thermo Fisher). All samples were analyzed and sorted using an Aria Fusion (Becton Dickinson) with 85 μm nozzle except for scRNA-seq samples, which were collected using a FACS Aria III (Becton Dickinson) with 100 μm nozzle. Flow cytometry data were analyzed using FlowJo v10.8 (Becton Dickinson).

Single-cell RNA-seq Library Preparation and Sequencing

Lung tissues were harvested from Tgfbr2^F/F and Tgfbr2^F/F;Nkx2.1-cre mice at P7 and P14 after hyperoxia or normoxia treatment. Two pups were harvested for each genotype, exposure, and timepoint for total of 8 mice on P7 and 8 mice on P14. Single cell suspension was obtained as described above. To enrich for epithelial and mesenchymal populations for sequencing input, 1 x 10⁵ CD45+ cells, 1 x 10⁵ CD31+ cells, and 1 x 10⁶ CD45−/CD31-cells were sorted for each sample and collected in 10% FBS/PBS. Sorted cells were then counted and labeled with oligonucleotide tags for multiplexing using 10x Genomics 3’ CellPlex Kit Set A. All 8 biologic samples for the P7 timepoint were pooled and 60,000 cells / lane were loaded onto 2 lanes of a Chromium Next GEM Chip (10x Genomics). The same workflow was used for P14 samples with all 8 biologic samples for P14 pooled together and loaded onto 2 lanes as well. Lanes 1 and 2 were therefore technical replicates for each biologic sample and similarly lanes 3 and 4 were also technical replicates. Chromium Single Cell 3’ v3.1 (10x Genomics) reagents were used for library preparation according to the manufacturer’s protocol. The libraries were sequenced on an Illumina NovaSeq 6000 S4 flow cell.

Sequencing Data Processing

Fastq files were uploaded to the 10x Genomics Cloud (https://www.10xgenomics.com/products/cloud-analysis) and reads were aligned to the mouse reference genome mm10. The data were demultiplexed, and cells with multiple oligonucleotide tags were identified as multiplets and removed by the 10x Genomics cloud analysis function with default parameters. Raw count matrices were imported to the R package Seurat v4.1.1 ⁸¹, and cells with fewer than 200 detected genes, more than 7500 detected genes, or more than 15% mitochondria genes were excluded. We used the DoubletFinder package⁸² for each sample using an estimated multiple rate of 1% remove doublets that were not detected upon alignment. We then merged all the sample objects, identified the top variable genes using Seurat’s FindVariableGenes function, and integrated the samples using the RunFastMNN function of the SeuratWrappers R package ⁸³. For visualization, Seurat’s RunUMAP function was performed using MNN dimensional reduction.

Twenty-three clusters were initially identified from a total of 27,839 cells using Seurat’s FindNeighbors and FindClusters functions with resolution = 0.3. The expression of canonical linage markers (Epcam, Col1a1, Pecam1, Ptprc, Msln) was used to define major cell types of the lung (epithelium, mesenchyme, hematopoietic, endothelium, mesothelium). Cluster 9 (729 cells) was identified as a contaminant and excluded from further analysis due to cells showing up in multiple locations across the UMAP embedding. While re-clustering and annotating epithelial, hematopoietic, and endothelial subpopulations we identified an additional 500 cells that clustered independently from known cell types within each subcluster, had less than 1000 detected genes, and were enriched for the expression of multiple canonical cell types, suggesting these cells were also contaminants. After removing these cells and cluster 9 cells, the remaining 26,610 cells were re-clustered with FindVariableGenes, RunFastMNN, RunUMAP, FindNeighbors, and FindClusters functions with clustering resolution = 0.8. Thirty clusters were identified at this stage and differentially expressed genes for each cluster were identified using FindAllMarkers focusing on genes expressed by more than 20% of cells (either within or outside of a cluster) and with a log fold change greater than 0.2. Using publicly available data, we were able to merge and annotate clusters to obtain the 25 clusters depicted in this study. Mesenchymal cells were re-clustered using the same workflow outlined above with clustering resolution = 0.3. Clusters were annotated based on a combination of three publicly available data sets for the neonatal lung. For Qiagen Ingenuity Pathway Analysis (Qiagen, https://digitalinsights.qiagen.com/IPA) ⁵¹, differentially expressed genes in myofibroblasts were identified by comparing either CTRL RA vs O2 cells or RA CTRL vs cKO cells using Seurat’s FindMarkers function. These lists were uploaded to Qiagen IPA to identify predicted upstream regulators for each comparison. To identify shared predicted upstream regulators, lists were filtered for z-score >1.5 (upregulated) or z-score <-1.5 (downregulated), and the filtered results were compared by Venn diagram between the two comparisons.

NicheNet⁵³ was used to compare signaling between conditions, as described in the “Differential NicheNet analysis between conditions of interest” vignette on the GitHub repository. Standard NicheNet statistical thresholds described in the vignette were used in the analysis. Alveolar epithelial subtypes (AT2, AT2 Lyz1, AT2 Activated, AT1/AT2, and Prolif. AT2) were used as a combined sender cell population, and alveolar and ductal myofibroblasts were used as combined receiver cell populations. Differential signaling from sender-to-receiver cells was compared across two conditions, CTRL RA vs CTRL O2 and CTRL RA vs cKO RA, and the results were examined for pathways present in both comparisons.

Data Analysis

Mean linear intercept morphometry was quantified using ImageJ software as described previously ⁸⁰. scRNA-seq data analysis was performed in R version 4.1.1. Statistical tests were performed in GraphPad Prism version 9.4.0.

Data Availability

The data that support this study are available upon reasonable request. Sequencing data and Seurat objects were deposited in the Gene Expression Omnibus Series GSE243129.

Grants

I.S.K was supported by the Pediatric Scientist Development Program, the Nina Ireland Program for Lung Health, and the UCSF Division of Neonatology. This work was supported by HD000850 (NICDH, I.S.K.); UCSF Division of Neonatology (I.S.K.); and HL145037 and HL142568 (NHLBI, D.S.).

Acknowledgements

We thank members of the Sheppard lab for helpful advice and discussion. We thank the following UCSF core facilities for their technical support: Laboratory for Cell Analysis (supported by P30CA082103), Center for Advanced Light Microscopy (supported by UCSF PBBR), Gladstone Genomics Core, Gladstone Histology and Light Microscopy Core, Gladstone Transgenic Gene Targeting Core, and the Center for Advanced Technology (supported by UCSF PBBR, RRP IMIA, and 1S10OD028511-01).

Disclosures

D.S. is a founder of Pliant Therapeutics, a member of the Genentech Scientific Review Board, a member of the Amgen Immunology Scientific Advisory Board, and a member of the Scientific Review Board for Lila Biologics. No funding or reagents from any of these companies were used for this project. None of the other authors has any conflicts of interest, financial or otherwise, to disclose.

Author Contributions

I.S.K and D.S. conceived the study, interpreted the data, and wrote the manuscript. I.S.K., C.M., X.R., V.A., and M.C. performed the experiments and/or analyzed the data. T.T. and A.A. provided critical resources for this study. D.S. supervised the study.

Abbreviations

• α-SMA: alpha-smooth muscle actin

• BPD: bronchopulmonary dysplasia

• cHet: conditional heterozygous

• cKO: conditional knockout

• EdU: 5-ethynyl-2’-deoxyuridine

• PDGFRα: platelet-derived growth factor receptor alpha

• P: postnatal day

• scRNA-seq: single-cell RNA-sequencing

• TGFβ: transforming growth factor-beta

• TGFβR2: transforming growth factor-beta receptor 2

• UMAP: Uniform Manifold Approximation and Projection.

Supplementary Figure Legend

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.