Lethal lung hypoplasia and vascular defects in mice with conditional Foxf1 overexpression

ABSTRACT FOXF1 heterozygous point mutations and genomic deletions have been reported in newborns with the neonatally lethal lung developmental disorder, alveolar capillary dysplasia with misalignment of pulmonary veins (ACDMPV). However, no gain-of-function mutations in FOXF1 have been identified yet in any human disease conditions. To study the effects of FOXF1 overexpression in lung development, we generated a Foxf1 overexpression mouse model by knocking-in a Cre-inducible Foxf1 allele into the ROSA26 (R26) locus. The mice were phenotyped using micro-computed tomography (micro-CT), head-out plethysmography, ChIP-seq and transcriptome analyses, immunohistochemistry, and lung histopathology. Thirty-five percent of heterozygous R26-Lox-Stop-Lox (LSL)-Foxf1 embryonic day (E)15.5 embryos exhibit subcutaneous edema, hemorrhages and die perinatally when bred to Tie2-cre mice, which targets Foxf1 overexpression to endothelial and hematopoietic cells. Histopathological and micro-CT evaluations revealed that R26Foxf1; Tie2-cre embryos have immature lungs with a diminished vascular network. Neonates exhibited respiratory deficits verified by detailed plethysmography studies. ChIP-seq and transcriptome analyses in E18.5 lungs identified Sox11, Ghr, Ednrb, and Slit2 as potential downstream targets of FOXF1. Our study shows that overexpression of the highly dosage-sensitive Foxf1 impairs lung development and causes vascular abnormalities. This has important clinical implications when considering potential gene therapy approaches to treat disorders of FOXF1 abnormal dosage, such as ACDMPV.


INTRODUCTION
The FOXF1 (Forkhead box F1) gene located on chromosome 16q24.1 encodes a member of the FOX family of transcription factors characterized by a distinct forkhead DNA binding domain, and plays an important role in epithelium-mesenchyme signaling as a downstream target of Sonic hedgehog (Mahlapuu et al., 2001;Dharmadhikari et al., 2015). FOXF1 is expressed in fetal and adult lungs, placenta, and prostate (Hellqvist et al., 1996;Bozyk et al., 2011;van der Heul-Nieuwenhuijsen et al., 2009). In the mouse embryonic lungs, Foxf1 expression is restricted to mesenchymederived cells such as alveolar endothelial cells and peribronchiolar smooth muscle cells (Kalinichenko et al., 2001a). Lung mesenchymal cells include numerous subtypes such as airway smooth muscle cells, fibroblasts, pericytes, vascular smooth muscle cells, and alveolar endothelial cells. Alveolar endothelial cells are differentiated vascular cells of mesenchymal origin, and line blood vessels that are in close proximity to air spaces in the lung.
Heterozygous point mutations and genomic deletions involving FOXF1 or its upstream enhancer have been reported in newborns with a lethal lung developmental disorder alveolar capillary dysplasia with misalignment of pulmonary veins (ACDMPV, OMIM 265380) with or without defects involving heart, gastrointestinal, or genitourinary systems (Stankiewicz et al., 2009;Bishop et al., 2011;Sen et al., 2013). A majority of Foxf1 +/− mice die perinatally, exhibiting defects in lung vasculature, similar to those in patients with ACDMPV (Kalinichenko et al., 2001b;Stankiewicz et al., 2009).
FOXF1 has also been reported to be epigenetically inactivated in breast and colorectal cancers (Lo et al., 2010;Mitchell et al., 2014) and overexpressed in Patched-associated tumors, including basal cell carcinoma, medulloblastoma, and rhabdomyosarcoma (Wendling et al., 2008a;Armeanu-Ebinger et al., 2011). FOXF1 overexpression was found in lung fibroblasts from patients with idiopathic pulmonary fibrosis (Melboucy-Belkhir et al., 2014). Isolated gastrointestinal abnormalities such as pyloric stenosis, mesenterium commune, and aplasia of the appendix were found to be associated with 16q24.1 duplications involving FOXF1; however, RNA or protein lung studies could not be performed . Hence, while effects of the loss of function of FOXF1 in lung disease and cancer are well documented, effects of the overexpression of FOXF1, specifically in the context of lung development, are not currently known.
To study the consequences of Foxf1 overexpression, we developed a Cre-inducible Foxf1 allele by knocking-in Foxf1 into the ROSA26 (R26) locus after a lox-STOP-lox (LSL) cassette (ROSA26-lox-STOP-lox-Foxf1). Here, we show that Tie2-creinduced heterozygous R26Foxf1 mice manifest a wide range of phenotypes that cause death from embryonic to perinatal stages due to lung and vascular defects. Determining dosage sensitivity of Foxf1 is important to inform future gene therapy approaches to potentially treat patients with ACDMPV and other disease conditions due to genetic alterations in FOXF1.

Tie2-cre-mediated overexpression of Foxf1 results in embryonic vascular and perinatal lung defects
We generated a floxed Foxf1 construct to target the ROSA26 locus (R26-LSL-Foxf1; Fig. S1) as previously described (Hohenstein et al., 2008;Carofino and Justice, 2015). First, we tested overexpression of Foxf1 in all tissues by crossing the R26-LSL-Foxf1 and the CMV-cre lines. We found that it led to embryonic lethality of R26Foxf1; CMV-cre embryos, displaying hemorrhage at E12.5 (Fig. S2). Foxf1 is endogenously expressed in mesenchymederived cell types such as endothelial and smooth muscle cells in the lungs (Kalinichenko et al., 2001a). To determine the effects of overexpression of Foxf1 in a similar context, we mated R26-LSL-Foxf1 mice to Tie2-cre mice, overexpressing Foxf1 in the endothelial and hematopoietic lineages.
Lungs of R26Foxf1: Tie2-cre pups are hypoplastic At P0.5, the wet lung weight (LW) to body weight (BW) ratios of R26Foxf1; Tie2-cre pups were significantly lower (approximately 3.0%) compared to that of R26-LSL-Foxf1 control littermate pups (approximately 4.1%) (Fig. 3A). Computational analyses of micro-CT imaging of E15.5 and E18.5 R26Foxf1; Tie2-cre embryos confirmed that their lungs are hypoplastic compared to R26-LSL-Foxf1 control embryos ( Fig. 3B-E), and did not show mispatterning of major pulmonary vessels (Fig. S7). Furthermore, the LV:BW ratios were significantly lower (approx. 50%) in E18.5 R26Foxf1; Tie2-cre embryos compared to controls (Fig. 3F). There was also a trend towards a lower ratio of lung surface area to body weight (approx. 30% lower) (Fig. 3G). Other organs, including the heart and liver, did not show any hematomas or evidence of blood pooling. R26Foxf1: Tie2-cre pups exhibit respiratory defects Respiratory function in P0.5 pups was tested under room air, hypoxic, and hypercapnic conditions using head-out plethysmography (Fig. 4). Under room-air conditions, the R26Foxf1; Tie2-cre pups showed significantly higher respiratory rates (V f ) (Fig. 4A) while tidal volume per breath (V T ) was significantly lower than in sibling controls (Fig. 4B). Total respiratory output per minute (V E ) matched that of sibling controls (Fig. 4C). Examination of R26Foxf1; Tie2-cre respiratory traces showed a significant decrease in the number of apneas, but the average length of apneas was not significantly different ( Fig. 4D,E). Respiratory traces also showed that R26Foxf1; Tie2-cre animals had more regular breathing patterns, as they had significant reductions in measures of variability in interbreath interval (IBI) and breath volume ( When challenged with hypoxia, R26Foxf1; Tie2-cre animals showed a significant decrease in the relative change in respiratory rate as compared to sibling controls (Fig. 4H). As overall tidal volume did not significantly change as compared to room air values (Fig. 4I), R26Foxf1; Tie2-cre animals showed a trend of reduced relative increase in overall minute ventilation as well as reduced absolute V T and V E (Fig. 4J, Fig. S9B,C). R26Foxf1; Tie2-cre animals also showed a trend toward an increase in apneas but it was not significant (Fig. 4K). There was no difference in IBI or volume variability ( Tie2-cre animals maintained a steady breathing pattern upon hypoxic exposure, three out of the 11 animals ceased breathing (see Fig. 5E for example trace) but resumed upon re-exposure to room air.
Upon a mild 5% hypercapnic challenge, R26Foxf1; Tie2-cre animals showed a blunted response with a reduction in the relative increase in V f as compared to sibling controls (Fig. 4O), but no difference in relative change of V T (Fig. 4P) or V E (Fig. 4Q), resulting in a trend of reduced absolute V T and V E (Fig. S9E,F). Under the increased drive mediated by hypercapnia, there were no significant differences in the number or length of apneas (Fig. 4R,S) and no significant differences were found in IBI and volume variability (Fig. 4T,U).  (Table S1). Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Huang et al., 2009) analyses identified gene ontology (GO) terms related to protein transport, protein localization, cell adhesion, and blood vessel morphogenesis to be associated with the deregulated genes (Fig. 6B). The microarray data was verified for the genes Igfbp3, Pparg, and Rcan1 using qRT-PCR (Fig. S10). Using a threshold of 1.2-foldchange, comparison of the R26Foxf1;Tie2-cre microarray dataset with our previously published Foxf1 knock-out P0.5 lung dataset (Sen et al., 2014) revealed 215 genes to be commonly deregulated in both datasets (P-value for overlap: 6.5×10 −13 ; Fig. 6C). Blood circulation, blood vessel development, and lung development were biological processes of interest found to be associated with the Breeding data from ROSA26-LSL-Foxf1 homozygote female mice×+/Tie2-cre male mice suggest perinatal lethality due to a significant decrease in number of normal R26Foxf1; Tie2-cre mice at P0.5. commonly deregulated genes in both the datasets (Fig. 6D). 165 out of 215 genes exhibited reciprocal gene expression trends, indicating that they could be direct targets of FOXF1 (P-value for reciprocal gene trends: 7.6×10 −16 ; Fig. 6E; Table S3).

E18.5 FOXF1 ChIP-seq analysis reveals targets involved in vasculature development
Using a threshold of 5×10 −7 for irreproducible discovery rate (IDR), ChIP-seq analysis identified 697 significant peaks for FOXF1 binding sites in E18.5 wild-type lungs. 1033 genes were found to be associated with these peaks (Table S4). Distribution of the peaks with respect to transcription start sites (TSS) is shown in Fig. 7A, with most peaks located 50 to 500 kb upstream or downstream to the TSS. Representative peaks in both biological replicates compared to the input control near the gene Cdh5 are shown in Fig. 7B. Functional analysis using the Genomic Regions Enrichment of Annotations Tool (GREAT) (McLean et al., 2010) revealed enrichment of biological processes related to vasculature development, heart development, and embryonic development (Fig. 7C). Similar analysis for abnormal mouse phenotypes showed enrichment of phenotypes such as abnormal vascular, cardiovascular and embryonic development ( Fig. 7D). Overlaying the FOXF1 ChIP-seq data with the R26Foxf1;Tie2-cre and Foxf1 +/− lung microarray datasets identified 11 common genes: Arhgap18, Sox11, Zswim6, Tnfrsf19, Ednrb, Ghr, 2510009E07Rik, Ostf1, Smarca2, Slit2 and Nup210 (Fig. 7E). This layered and multi-step analysis further narrows down the list of potential direct targets of FOXF1 in the developing lung.

DISCUSSION
In contrast to Foxf1 haploinsufficiency, the phenotypic effects of the increased dosage of Foxf1 during murine embryonic development are largely unknown. Previous studies include Foxf1-enforced expression in vitro, shown to repress hematopoiesis (Fleury et al., 2015), and overexpression of Foxf1 in murine skin under the control of the basal cell-specific promotor Keratin 5 (Krt5), resulting in severe hair loss, increase in size and number of sebaceous glands, and growth failure of these animals (Wendling et al., 2008b). The R26-LSL-Foxf1 line described here is a novel mouse model to overexpress Foxf1 in a time-and tissue-specific manner.
Tie2-cre-mediated overexpression of Foxf1 beginning at E7.5 in endothelial and hematopoietic lineages resulted in a combination of pulmonary, vascular, lymphatic and platelet defects. At E15.5, 35% of the heterozygous R26Foxf1; Tie2-cre embryos showed localized subcutaneous hemorrhages and edema, indicating a vascular fragility defect. Moreover, whole-body overexpression of Foxf1 using CMV-cre also led to embryonic lethality, with R26Foxf1; CMV-cre embryos exhibiting hemorrhages. We elected to specifically focus on the effects of Foxf1 overexpression in the lungs because of the pulmonary phenotypes seen in patients with ACDMPV. The majority of R26Foxf1; Tie2-cre pups died within 24 h of birth, due to respiratory failure.
The lungs of the heterozygous R26Foxf1; Tie2-cre pups were smaller compared to those of their littermate controls. qRT-PCR showed a 1.7-fold overexpression of Foxf1 in E18.5 R26Foxf1; Tie2-cre lungs and an increase in the amount of the FOXF1 protein was also observed by IHC in E18.5 lungs. The lungs were found to be immature by histopathological evaluation and were associated with decreased expression of the vascular marker FLK1, but not the epithelial marker proSPC. In mice, a wet LW:BW ratio less than 4% defines pulmonary hypoplasia (Seegmiller et al., 1986), indicating that the R26Foxf1; Tie2-cre lungs with a wet LW:BW ratio of approximately 3% are hypoplastic. In support of this notion, micro-CT imaging showed that the E18.5 average lung volume and lung surface to body weight ratios of R26Foxf1; Tie2-cre embryos were approximately half and one-third smaller, respectively, compared to control lungs.
Examination of respiratory function revealed multiple deficits under baseline (room air), hypoxic and hypercapnic conditions. Under baseline conditions, P0.5 R26Foxf1; Tie2-cre mice had a reduced V T , likely due in part to the reduced volume of the hypoplastic lungs. However, an increase in V f was compensatory and resulted in equivalent minute ventilation V E . The increase in respiratory drive is also the most probable reason behind the reduced number of apneas and reduced periodic and volume instability in comparison to the typical irregular breathing pattern of sibling controls (Hilaire and Duron, 1999), with no additional change to apnea length. Despite the equivalent ventilatory output, R26Foxf1; Tie2-cre pups still exhibit lethality, suggesting  (n=11) show a reduced relative increase in respiratory rate compared to baseline room air respiratory rate (H, P=0.03), similar relative tidal volume increase (I), and a trend towards a reduced relative increase in overall minute ventilation (J, P=0.07) as compared to sibling controls (n=10). R26Foxf1; Tie2-cre neonates show similar numbers of apneas (K), and matched average apnea length (L) and levels of instability in both inter-breath interval (IBI) and tidal volume (M,N). (O-U) Under hypercapnia (5% CO 2 /21% O 2 /74% N 2 ), R26Foxf1; Tie2-cre neonates (n=9) show a reduced relative increase in respiratory rate (O, P=0.04), matched increase in tidal volume (P), and a statistically insignificant reduced relative increase in minute ventilation (Q) as compared to sibling controls (n=10). R26Foxf1; Tie2-cre pups also show matched number (R) and length of apneas (S) and levels of instability in both IBI and volume (T,U). Data are shown as mean±s.e.m.
that changes in respiratory rate are insufficient to overcome the observed pulmonary and vascular defects, an observation supported by our hypoxic studies. Whereas the Tie2-cre is expressed in the neurovasculature, central-or peripheral-neuralor glia-specific expression has not been reported or observed in our hands. We hypothesize that increasing pulmonary-vascular dysfunction results in failure of adaptive regulation through greater levels of oxygen desaturation that ultimately leads to cardiopulmonary failure.
Microarray analyses of the E18.5 R26Foxf1; Tie2-cre lung transcriptome showed enrichment of genes related to protein transport, protein localization, blood vessel development, and cell adhesion. Some of the downregulated genes included Igfbp3, Pparg, Rcan1 and Prkcdbp. IGFBP3 and PRKCDBP are also inversely upregulated in ACDMPV lungs (Sen et al., 2014), suggesting that these genes might be relevant to the role of FOXF1 in the pathology of ACDMPV. Downregulation of Igfbp3 is associated with nitrofen-induced pulmonary hypoplasia (Ruttenstock et al., 2010). Interestingly, Pparg agonists enhance lung maturation and attenuate hypoxiainduced inhibition of lung development Nicola et al., 2011), whereas Rcan1 is implicated to play a role in angiogenesis as a regulator of calcineurin (Riper et al., 2008), and is activated by Vegf (Holmes et al., 2010). Some of the upregulated genes included Fgfr2, Egfl7 and Robo4. Fgfr2 upregulation has been described to be associated with nitrofen-induced pulmonary hypoplasia (Friedmacher et al., 2012). Egfl7 overexpression in mice has been reported to cause embryonic lethality due to impaired angiogenesis (Nichol et al., 2010), whereas Robo4 is a vascular specific receptor known to inhibit endothelial migration (Park et al., 2003). Furthermore, comparison of the genes deregulated in the R26Foxf1; Tie2-cre lung dataset with the genes differentially expressed in the dataset associated with the heterozygous loss of Foxf1 in P0.5 lungs (Sen et al., 2014) showed involvement of genes associated with a variety of biological processes, including blood circulation, vasculature development and lung development. Of note, 165 genes showed reciprocal trends in gene expression, when compared between the two datasets. These genes are potential targets of FOXF1 in the lung, as their expression changes reciprocally with the loss or gain of Foxf1.
ChIP-seq analysis in E18.5 wild-type lungs identified binding of FOXF1 in the proximity of genes involved in biological processes such as blood vessel, cardiovascular and embryonic development. Among the genes associated with multiple FOXF1 binding sites, Nrp1 signaling has previously been shown to be essential for fetal pulmonary development (Joza et al., 2013). Interestingly, mice homozygous for a knock-out allele of Sdpr (Cav2) exhibit abnormal caveolae formation in the lung endothelium affecting endothelial cell function (Hansen et al., 2013). Similarly Cdh5 and Itgb1 are genes involved in endothelial cell development. Recently, it was shown that Itgb1 controls Cdh5 localization and blood vessel stability (Yamamoto et al., 2015). Another forkhead gene, Foxa2, show variability in both interbreath interval (IBI) and tidal volume with high numbers of apneas. In contrast, R26Foxf1; Tie2-cre neonates (B) show increased respiratory rate, reduced tidal volume, and more regular breathing with fewer apneas. (C-E) Under hypoxic conditions, apart from changes in amplitude, sibling controls (C) and R26Foxf1; Tie2-cre neonates (D) have similar respiratory traces, however, three out of eleven R26Foxf1; Tie2-cre animals ceased breathing (see E for example) and resumed once they were returned to room air conditions. (F,G) Under hypercapnic conditions, apart from differences in amplitude, R26Foxf1; Tie2-cre neonates and their sibling controls show similar patterning in breathing. showing gene expression trends between the commonly regulated genes in the two datasets. The light blue bar represents genes upregulated in the Foxf1 knock-out dataset and downregulated in the Foxf1 knock-in dataset; dark blue bar represents genes down regulated in the knock-out dataset and upregulated in the knock-in dataset; orange bar represents genes upregulated in both datasets, and the red bar represents genes downregulated in both the datasets.  Biology Open (2016) 5, 1595-1606doi:10.1242 has been described to be required for the transition to breathing at birth (Wan et al., 2004). When the FOXF1 ChIP-seq data were compared to the Foxf1 knock-in and knock-out microarray datasets, 11 genes were found to have binding sites for FOXF1 and were reciprocally deregulated in the microarray datasets. These included the genes Arhgap18, Sox11, Zswim6, Tnfrsf19, Ednrb, Ghr, 2510009E07Rik, Ostf1, Smarca2, Slit2 and Nup210. Interestingly, Sox11 −/− mice exhibit lung hypoplasia and die at birth (Sock et al., 2004). Ghr signaling is involved in early lung growth, oxidative protection, and lipid metabolism in the developing lung (Beyea et al., 2006). These genes could potentially be direct targets of FOXF1 in the embryonic and early postnatal lung.

RESEARCH ARTICLE
In conclusion, the R26-LSL-Foxf1 mice develop lung and vascular defects when crossed to the Tie2-cre line that activates recombination in endothelial and hematopoietic lineages. The defects manifest as clear functional deficiencies in neonate respiratory function. Our study shows that Foxf1 is highly dosage sensitive, with both loss-and gain-of-function of Foxf1 having implications in development and disease conditions. Additionally, the dosage sensitivity of Foxf1 suggests that conventional gene therapy approaches to treat ACDMPV and other FOXF1 related disorders may not be successful. Instead, alternative approaches that manipulate targets downstream of FOXF1 would need to be investigated.

Animal care
All mouse experiments were carried out under the approval of the Institutional Animal Care and Use Committee (IACUC) at Baylor College of Medicine (BCM). Mice were housed in the Transgenic Mouse Facility (barrier level 3) under the care of the Center for Comparative Medicine (CCM), which is accredited by the Association for Assessment and Accreditation of Laboratory Animal Care International (AAALAC).

Generation of the ROSA26 targeting construct
The ROSA26 locus is ubiquitously expressed; disruption of the endogenous transcript has no apparent phenotype (Zambrowicz et al., 1997) and is widely used as a reporter locus to determine expression from tissue-specific cre drivers (Srinivas et al., 2001) and to model targeted overexpression of oncogenes (Carofino et al., 2013).

Chimeric mouse generation and breeding
Targeting vector DNA was electroporated into the C57BL/6N embryonic stem cell (ESC) line JM8A3 (Pettitt et al., 2009) by the Mouse Embryonic Stem Cell Core at BCM. Correctly targeted clones were microinjected into C57BL/ 6-Tyr c-Brd (albino) blastocysts and transplanted into pseudopregnant foster mothers by the Genetically Engineered Mouse (GEM) Core at BCM. Chimeric male offspring were crossed to C57BL/6-Tyr c-Brd females to test for germline transmission of the targeted ROSA26 allele, referred to here as R26-LSL-Foxf1. R26-LSL-Foxf1 animals were crossed to the B6.Cg-Tg(Tek-cre)1Ywa/ J (Tie2-cre) line (Kisanuki et al., 2001), obtained from the laboratory of Dr. Daryl Scott (BCM) and the B6.C-Tg(CMV-cre)1Cgn/J line (Schwenk et al., 1995). The Tie2-cre line is expressed in endothelial and hematopoietic cell lineages beginning at E7.5 and the X-linked CMV-cre line is likely activated before implantation during early embryogenesis and is expressed in all tissues, including germ cells.

Endothelial cell flow sorting
E18.5 lungs (R26Foxf1; Tie2-cre and control +; Tie2-cre) were dissected in 3 cm dishes and minced finely with forceps and scissors. Then 2 ml of the digestion medium (HBSS supplemented with 0.05% Trypsin, 0.1% Collagenase and 25 µM HEPES) was added to the dish and the tissue was incubated for 1 h at 37°C. The reaction was stopped by addition of cold MEM supplemented with 10% FBS. The tissue was then disaggregated by pressing through a 40 µm nylon mesh. Cells were centrifuged at 290 g (1200 rpm) for 10 min and the pellet was re-suspended in 250 µl PBS supplemented with 1 mM EDTA, 25 mM HEPES, and 1% FBS to ∼10 6 cells/ml. Endothelial cells were stained with anti-CD31 (FITC conjugated) antibodies (Abcam, Cambridge, MA, USA) at a 1:50 dilution for 1 h at 4°C in the dark. Washing was done three times in the same buffer that was used for staining. Sorting of CD31-positive cells was done on FACSAriaII cell sorter (BD Biosciences, San Jose, CA, USA) using a 70 µm nozzle. Cells were collected directly into QIAzol and processed for RNA isolation using miRNeasy Mini Kit (Qiagen, Hilden, Germany).

Hematology measurements
For hematology measurements trunk blood was obtained from P1.5 pups by decapitation. CBCs were performed by the Comparative Pathology Laboratory at BCM on an Advia 120 automated hematology system (Siemens) on samples pooled from 2-3 pups.

Micro-CT imaging
E15.5 and E18.5 embryos were imaged by micro-CT at the Optical Imaging and Vital microscopy (OIVM) core at BCM. Embryos were prepared by a method called STABILITY as described previously (Wong et al., 2013). In short, embryos were fixed in 4% PFA; hydrogel stabilized using acrylamide and stained with 0.1N iodine. 3D datasets for each embryo were acquired using SKYSCAN 1272 micro-CT scanner (Bruker). Images were obtained using a 0.5 mm aluminum filter with a rotation step of 0.2 at a resolution of 11 μm. The acquired datasets were reconstructed using the NRecon software (Bruker) and visualized using the CTVox software (Bruker).

Lung volume quantification
Individual lung images were extracted using the region of interest (ROI) function of the CTAn software (Bruker) from reconstructed NRecon datasets to create volume of interest (VOI) datasets for each embryo. The VOI lung datasets were visualized using the Imaris software (BitPlane) and lung volumes calculated using the surface rendering function in Imaris.
Neonatal head-out pneumotachography P0.5 respiration was measured in a custom built head-out maskpneumotachograph system as described previously (Cummings and Frappell, 2009), that was engineered and machined for a minimum of dead space to increase sensitivity. Additional facemask ports were engineered for gas flow-through and calibration. For calibration and experiments, room air or mixed gasses were drawn through the mask-pneumotachograph by a vacuum pump attached to the gas flow-through port. Mixed gasses were supplied by a magnetically coupled bell housing that rested over the end of the pneumotachograph while remaining partially open to the room to avoid pressure fluctuations from the inflowing gas. All measurements were done between 7 am and 3 pm on the day of birth.
Prior to an experiment, the facemask was sealed with a piece of nitrile rubber. Ventilation was calibrated as a series of 20 µl pipetman injections into an empty facemask at a rate of 3 Hz. The rate of gas flow-through was continuously controlled via two rotameters in series.
For experimental assays, a small opening was made in the nitrile rubber to fit the snout (nose and mouth) of a P0.5 mouse. The mouse was affixed to the facemask with Impregum F, Polyether Impression material (Patterson Dental, St. Paul, MN, USA). The mouse rested on a platform attached to the facemask that fit inside a temperature controlled chamber to maintain the mouse pup at 36°C.
Pneumotachograph pressure changes and chamber temperature were recorded with LabChartPro (AD Instruments, Colorado Springs, CO, USA) in real time. The pneumotachograph trace was integrated to produce a respiratory waveform. Waveforms were analyzed offline to determine respiratory rate (V f ), tidal volume (V T ), minute ventilation (V E ) and pattern analysis.
After attachment to the face mask, mice were allowed to acclimate for 10 min in room air. Data was recorded for another 20 min under room air conditions and then switched for an additional 20 min to either a mix of 10% O 2 /90% N 2 for hypoxia or a mix of 5% CO 2 /21% O 2 /74% N 2 for hypercapnia before returning to room air for 20 min.
Respiratory waveforms were collected when the neonate was at rest and readings were free from movement artifacts. A minimum of 1 min cumulative data compiled from traces at least 10 s long from the last 5 min of a given experimental condition were analyzed. Apart from integration, no filtering, smoothing or other manipulations were applied to the pressure waveform. Tidal volume (V T ) was determined by comparing peak (mV) height to calibration injections (mV/µl).
Poincaré plots and apnea measurements were determined using 10 min of movement-free traces from each breathing condition. Apneas were defined as an interbreath interval (IBI) that was longer than 1.5 s. The coefficient of variation (CV) of the IBI and amplitude was also calculated from the same 10-minute trace compilation of each breathing condition (standard error IBI or amplitude/mean IBI or amplitude).

Gene expression arrays
Illumina mouse WG-6 v2.0 expression Beadchip analysis was performed on three R26Foxf1;Tie2-cre E18.5 lung RNA samples and three control R26-LSL-Foxf1 E18.5 lung RNA samples. 500 ng of total RNA was labeled using Illumina TotalPrep RNA Amplification Kit (Thermo Fisher Scientific, Waltham, MA, USA) and hybridized as per manufacturer's instructions. The array data were analyzed using the lumi bioconductor package (Du et al., 2008), normalized by robust spline normalization and transformed using variance stabilization transformation (VST) as previously described (Sen et al., 2014). A two-sample t-test was applied to determine differentially expressed genes between R26Foxf1; Tie2-cre and the R26-LSL-Foxf1 lung groups. Differential expression P-values were adjusted for false discovery rates (FDR). Fold changes were calculated using reverse VST. Although our sample size was relatively small for parametric tests, we justified the use of t-tests in this study by large effect sizes of our analysis (DeWinter, 2013). Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Huang et al., 2009) was used for gene ontology and pathway analyses.

Quantitative RT-PCR
Total lung RNA was reverse transcribed to cDNA using the SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific, Waltham, MA, USA) and amplified using the Power SYBR Green PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA). Primers used for qRT-PCR validation of the mouse microarray are listed in Table S2. Amplification and data analysis were conducted on an ABI 7900HT Fast real time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). Relative gene expression was calculated using the ΔΔC T method (Pfaffl, 2001). To determine Foxf1 overexpression, TaqMan probes (Thermo Fisher Scientific, Waltham, MA, USA) were used for Foxf1 and Gapdh (internal control). Amplification was performed using the TaqMan Universal PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA).

ChIP-seq assay and analysis
Two biological replicates of pooled E18.5 lungs (n=3) were crosslinked with 37% formaldehyde to a final concentration of 1%. ChIP was performed using the SimpleChIP Enzymatic Chromatin IP Kit (Cell Signaling, 9005S). Micrococcal nuclease-digested chromatin was further sonicated using Diagenode Bioruptor with 20 pulses of 15 s on and 15 s off at high power to yield sheared chromatin. Ten micrograms of chromatin was used per immunoprecipitation with 7 µg of FOXF1 antibody  (AF4798: R&D Biosystems, Lot#B1508). ChIP-grade normal rabbit IgG #2729 and Histone H3 (D2B12) XP ® Rabbit mAb (ChIP formulated) from Cell Signaling Technology were used as negative and positive controls, respectively; 2% input was used as a control. The Genomic and RNA Profiling (GARP) Core at BCM conducted sample quality checks using the NanoDrop spectrophotometer, Invitrogen Qubit 2.0 DNA quantitation assay and Agilent Bioanalyzer. The Rubicon ThruPlex DNA-Seq library preparation system was used to prepare ChIP-Seq libraries for sequencing on the Illumina HiSeq sequencing system. Sequence reads were mapped on to the mm10 genome using Bow-tie2 (Langmead and Salzberg 2012). Percentage of uniquely mapped reads were 77.85, 73.10 and 73.10%, respectively, for the two biological replicates and input control samples, acceptable according to Bailey et al. (2013). Peak calling procedure was performed using the Model-based Analysis of ChIP-Seq (MACS2) tool (Zhang et al., 2008). The assessment of peak calling was done using IDR-score (Li et al., 2011). The consensus FOXF1 binding motif was identified from the DNA sequences underlying the peaks, using the WebLogo tool (http://weblogo.berkeley.edu/logo.cgi). Functional analysis of selected peak regions was performed using the GREAT tool (McLean et al., 2010) and classification of genes common to the ChIP-seq and microarray datasets was done using LungGENS (Du et al., 2015).

Data analysis
Graphical and statistical analyses were conducted using Prism (GraphPad, La Jolla, CA). A χ 2 test was used for analysis of timed matings and an unpaired t-test was used for all other analyses. Standard error of the mean is shown on all figures. Statistical significance was considered at P≤0.05.