Temporal, spatial, and phenotypical changes of PDGFRα expressing fibroblasts during late lung development^☆

2.1. PDGFRα⁺ fibroblasts switch from a proximal bronchiolar to distal alveolar location

To assess temporal-spatial localization of PDGFRα fibroblasts and identify their potential interactions with other cell subsets, we performed immunofluorescence microscopy on lungs of Pdgfrα^GFP knock in mice from E16.5 to PN28. Distribution of PDGFRα-GFP⁺ fibroblasts along the proximal-distal axis of the lung was determined by immunofluorescence using the epithelial marker E-Cadherin, confocal microscopy with 3D reconstruction. Differentiation of PDGFRα-GFP⁺ fibroblasts into smooth muscle cells and myofibroblasts was assessed by colocalization of the nuclear GFP with αSMA (Fig. 1).

During the canalicular stage (E16. 5–17.5) a dense cell layer of PDGFRα-GFP⁺ fibroblasts are located along bronchi and terminal bronchioles with the cells closest to the epithelium co-expressing αSMA. PDGFRα-GFP⁺ fibroblasts line the dilating acinar tubules but are absent from the more distal acinar buds (Fig. 1A, E). During the saccular stage (E18.5 - PN4) peribronchiolar PDGFRα-GFP⁺ fibroblasts are sparse, and predominantly negative for αSMA. The majority of PDGFRα-GFP⁺ fibroblasts associate with the terminal saccular structures, which will give rise postnatally to the alveolar ducts and alveoli (Fig. 1B, F). During the alveolar phase (PN5-PN28), PDGFRα-GFP⁺ fibroblasts become more restricted to the alveolar structures. At PN7, a significant number of interstitial PDGFRα-GFP⁺ fibroblasts activate expression of αSMA and differentiate into alveolar myofibroblasts (Fig. 1C, G). By PN28, individual PDGFRα-GFP⁺ fibroblasts are scattered throughout the alveolar interstitium; the relative number is reduced compared to other time points, and co-expression of αSMA is largely lacking (Fig. 1D, G, H). Assessment of Pdgfra^GFP lung single cell suspensions confirmed that during distal lung development, as a percentage of total lung cells, PDGFRα-GFP⁺ fibroblasts significantly diminish (Fig. 1H). These spatiotemporal relationships are consistent with previous findings (Branchfield et al., 2016; McGowan et al., 2008; Ntokou et al., 2015). These data demonstrate that PDGFRα fibroblasts form a layer along the conducting airways directly beneath the epithelium and contribute to peribronchiolar αSMA⁺ smooth muscle cells during terminal branching and sacculation. These peribronchiolar smooth muscle cells form a concentric pattern around mature airways, which is distinct from the perivascular mesh pattern found in mature arteries (Data in Brief Fig. 1). During saccular maturation and alveolarization, individual PDGFRα fibroblasts are localized to terminal bronchioles, terminal saccules and alveolar structures. In both compartments, and at different developmental time points, bronchiolar and alveolar PDGFRα⁺ fibroblasts can differentiate into contractile cells expressing αSMA.

2.2. Spatial relation of PDGFRα-GFP⁺ fibroblasts to distinct epithelial subpopulations

Apposing layers of epithelial and mesenchymal cells, separated by a basement membrane and the capillary network, form epithelial-mesenchymal trophic units (Brewster et al., 1990; Evans et al., 1999). In the lung, two distinct epithelial-mesenchymal trophic units have been identified: 1) proximal conducting airways consisting of bronchiolar cells and 2) distal alveolar structures containing alveolar epithelial cells. The cells that comprise the mesenchymal layer of these distinct units have not been defined. To further assess temporal and spatial proximity of PDGFRα-GFP⁺ fibroblasts with distinct epithelial subpopulations, we performed immunohistochemistry on Pdgfra^GFP lungs for epithelial-specific transcription factors SOX2, HOPX and SOX9 (Que et al., 2009; Barkauskas et al., 2013; Ng et al., 1997; Perl et al., 2005; Rockich et al., 2013; Figs. 2–5). SOX2 is expressed in proximal bronchiolar cells and located in the nucleus. HOPX, is expressed in bipotential pneumocytes located in the acinar tubules during the canalicular stages and becomes a more restricted marker specific for type 1 alveolar epithelial cells (AECI) in adult lungs (Jain et al., 2015). HOPX can be found in the cytoplasm and nucleus. SOX9 is restricted to distal acinar buds and is essential for early branching morphogenesis, distal epithelial differentiation, and extracellular matrix organization (Perl et al., 2005; Rockich et al., 2013; Ustiyan et al., 2016). The SOX9 protein is localized to the nucleus and is also rarely expressed in some fibroblasts adjacent to the mature peribronchiolar smooth muscle layer.

2.3. Canalicular lung at E16.5

SOX2 expression in the respiratory epithelium gradually decreases from proximal to distal and correlates with expression of αSMA in the underlying mesenchyme. Proximal bronchiolar epithelial cells that express SOX2 protein are in close apposition to a discrete layer of peribronchiolar smooth muscle cells consisting of double-positive PDGFRα-GFP⁺/αSMA⁺ fibroblasts (Fig. 2A). The region of the bronchiolar tree where αSMA expression in the underlying mesenchyme tapers off and gradually vanishes, identifies a bronchiolar transition zone that is continuous with the acinar tubules. In the late canalicular stage, the epithelium in this bronchiolar transition zone expresses HOPX, and the underlying fibroblasts form a PDGFRα-GFP⁺ mesenchymal layer in which αSMA is greatly diminished and absent more distally (Fig. 2B). Acinar buds are the most distal structures and that ultimately will mature into alveolar ducts and alveoli. During the canalicular stage, the epithelial cells of the acinar buds express SOX9 and the underlying mesenchyme contains almost no PDGFRα-GFP⁺ fibroblasts (Fig. 2C). The lipofibroblast phenotype was assessed by immunohistochemistry for perilipin (ADFP). Perilipin staining was not co-localized with nuclear PDGFRα-GFP, suggesting the absence of lipofibroblasts in the canalicular PDGFRα⁺ fibroblast population (Fig. 2D). Proximity of PDGFRα⁺ fibroblasts to capillaries and vessels was assessed by immunohistochemistry for endomucin. At E16.5, PDGFRα-GFP⁺ fibroblasts form a clear layer between bronchiolar epithelial cells and stromal cells adjacent to endomucin positive cells of the vasculature (Fig. 2E).

2.4. Saccular lung at E18.5

During the saccular stage (E18.5), SOX2-expressing bronchiolar epithelial cells are in juxtaposition to an αSMA⁺ smooth muscle layer that contains only a few isolated PDGFRα-GFP⁺ fibroblasts (Fig. 3A). In the terminal saccules, which will give rise to the alveolar ducts and alveoli, HOPX⁺ epithelial cells and PDGFRα-GFP⁺ fibroblasts are in close proximity (Fig. 3B). Moreover, PDGFRα-GFP⁺ fibroblasts in the terminal saccules are negative for αSMA. Compared to E16.5, SOX9 expression is significantly more restricted to individual epithelial cells in the developing saccules, which are rarely adjacent to PDGFRα-GFP⁺ fibroblasts (Fig. 3C). Perilipin and PDGFRα-GFP⁺ did not co-localize (Fig. 3D). In contrast to E16.5, PDGFRα-GFP⁺ fibroblasts are now adjacent to endomucin positive cells (Fig. 3E).

2.5. Alveolar lung at PN7

At PN7, SOX2-expressing epithelial cells in the bronchiolar alveolar duct junction are not associated with PDGFRα-GFP⁺ fibroblasts (Fig. 4A). In the alveolar epithelium, HOPX is expressed in AECIs that comprise the alveolar wall (Fig. 4B), while SOX9 expression is restricted to single cells that are predominantly located in corners that are shared amongst multiple alveoli (Fig. 4C). PDGFRα-GFP⁺ fibroblasts are located in close proximity to AECI but not SOX9⁺ cells. The majority of PDGFRα-GFP⁺ fibroblasts reside within the bases and tips of secondary crests, where endothelial cells (Fig. 4E), epithelial cells, PDGFRα-GFP⁺/αSMA⁺ myofibroblasts, and PDGFRα-GFP⁺/ADFP⁺ lipofibroblasts are positioned intimately (Fig. 4D–F). Closer inspection of developing secondary crests by electron microscopy revealed fibroblasts containing lipid vesicles at the base of the emerging septum (presumed lipofibroblasts), and elastin deposition in fibroblasts residing at the septal tip (presumed matrix fibroblasts). (Fig. 4F, note no myofibrils detected in this image, graphic overlay is for illustration purposes only).

2.6. End of alveolarization PN28

Toward the end of the alveolar stage at PN28, SOX2 is expressed in the bronchiolar epithelium, including the epithelium of the bronchiolar-alveolar duct junction (Fig. 5A). Some PDGFRα-GFP⁺ fibroblasts are found adjacent to but predominately distinct from the αSMA⁺ smooth muscle layer (Fig. 1D, insert). This spatial proximity of PDGFRα-GFP⁺ fibroblasts to HOPX⁺ and SOX9⁺ alveolar epithelial cells is comparable to PN7 (Fig. 5B, C). Contrary to the earlier period of alveolar development, where there was robust expression of αSMA, now there are only sparse interstitial PDGFRα-GFP⁺/αSMA⁺ myofibroblasts but there are clearly identifiable interstitial PDGFRα-GFP⁺/ADFP⁺ lipofibroblasts (Fig. 5B–D).

2.7. Dynamic changes in gene expression profiles in PDGFRα positive fibroblasts from the canalicular to the alveolar stages

PDGFRα-fibroblasts compartmentalize with proximal bronchiolar epithelial cells in the canalicular stage and later with distal alveolar epithelial cells during saccular and alveolar development and exhibit different functions co-expressing αSMA and perilipin. We therefore compared their gene expression profiles throughout the canalicular, saccular and alveolar stages. PDGFRα⁺ fibroblasts were isolated from mouse lungs at E16.5, E18.5, PN7 and PN28, and RNA obtained from the resulting enriched PDGFRα population was subjected to RNA-Seq analysis. Bayesian Analysis of Time Series (BATS) was performed and 1908 genes were identified that were differentially expressed across these developmental stages (Bayesian Factor < 0.01) (Angelini et al., 2008; Fig. 6A). Hierarchical clustering and principle component analyses demonstrated clear separation of gene expression profiles in PDGFRα⁺ fibroblasts isolated at different stages of distal lung development (Fig. 6B, C). Differentially expressed genes were subjected to pattern recognition using Short Time series Expression Miner, STEM (Ernst and Bar-Joseph, 2006). STEM clustered 1607 genes into 7 major gene expression profiles (Fig. 6D–J). Unique functional enrichments of genes within each profile were identified by comparing all profiles to each other using Toppcluster (See Tables 1 and 2) (http://toppcluster.cchmc.org/) (Kaimal et al., 2010).

2.8. Gene expression and functional enrichment at E16.5

Profile #10 consists of genes that are gradually downregulated throughout distal lung development and it inversely correlates with profiles #39 and #23 which consist of genes that are gradually upregulated (Fig. 6F, J). To account for molecular agonists and antagonists of biological processes we looked for shared functional enrichment categories of either up or downregulated genes. With this approach, we identified gene enrichment for negative regulation of cell migration, regulation of response to wounding, and negative regulation of cell death (See Table 2). These profiles suggest that cell migration decreases over time; that developmental processes are also used in responses to wounding; and that cell death is negatively regulated upon completion of distal lung development.

2.9. Gene expression and functional enrichment at E18.5

During the saccular stage at E18.5, the majority of profiles indicate that gene expression is relatively downregulated (Profile #20 and #13). Only 66 of the 1607 genes are significantly upregulated at E18.5 when compared to all other time points (Data in Brief Table 3). Unique gene enrichment categories for these genes include “protein kinase inhibitor activity” and “delayed wound healing.” Co-expression analysis revealed high correlation with “genes upregulated in cultured stromal stem cells from adipose tissue.” These gene expression data suggest that during saccular lung development, cells are less mitotic and are shifting from the smoothmuscle phenotype of the canalicular bronchioles to the interstitial matrix- and lipofibroblasts.

2.10. Gene expression and functional enrichment at PN7

Profiles #1, #18 and #39 comprise genes that are upregulated specifically during active alveolar septation. Additionally, genes in #39 increase until PN28. Analysis revealed that expression is enriched for signal transduction, actin cytoskeleton, ECM organization, cell migration, angiogenesis and epithelial development (Table 1). These enrichment categories identify myo-, matrix-, and lipo-fibroblast functions within the PDGFRα⁺ population, and imply that PDGFRα⁺ interstitial fibroblasts play an important role in regulating both angiogenesis and epithelial differentiation during alveolar septation.

2.11. Gene expression and functional enrichment at PN28

Profiles #39, #13, and #23 are comprised of genes that have high expression in the late phases of distal lung development. Expression profile #39 contains genes that steadily increase throughout lung development. It consists of 271 genes, 18 of which encode transcription factors. The major unique enriched gene categories include “cell motility”, “cell matrix adhesion”, “tube development” and “renal development”. Profiles #13 and #23 include genes that are increasing from E18.5 onwards. These gene ontology categories contain functional terms that suggest differentiation of lipid and matrix fibroblasts (“regulation of lipid metabolic process”, “vesicle mediated transport”, and “Golgi membrane”) (Eyden, 2008; Tsupykov et al., 2016). The predominance of PDGFRα-GFP⁺/Perilipin lipofibroblasts that we detected at PN28 corroborates these observations (Fig. 5D).

2.12. Smooth muscle versus myofibroblasts

Despite the use of single or combinations of expression markers, including PDGFRα and αSMA, the ability to identify the distinct developmental functions of the heterogeneous population of lung fibroblasts remains challenging. Previously, single cell RNA-Seq analysis performed by the LungMAP consortium of E16.5 and E18.5 lungs classified individual fibroblast populations based on unique gene expression signatures and named them based on functional gene enrichment categories. Searching this single cell RNA-Seq dataset for Pdgfra and Acta2 expressing cells identified these smooth muscle and myofibroblast progenitors in 4 different fibroblast populations with distinct functional phenotypes, Myofibroblast/SmoothMuscle, Proliferative Mesenchymal Progenitor cells, Intermediate Fibroblasts and MatrixFBs (Du et al., 2015; Guo et al., 2016). Our present work expands upon this data by examining gene expression patterns of PDGFRα expressing fibroblasts throughout the spectrum of distal lung development and comparing/contrasting patterns of gene expression at each distinct stage. Additionally, we employed immunohistochemistry and confocal microscopy to map the unique locations of these cells in the developing lung to better understand their individual microenvironments and infer how cell-cell interactions might specify each population’s distinct function. Our data demonstrate that at E16.5, PDGFRα-GFP⁺/αSMA⁺ smooth muscle cells are located proximally and adjacent to the SOX2⁺ bronchiolar epithelium (insert Fig. 1 A and Fig. 2). Starting at E18.5, expression of PDGFRα-GFP can be detected in the developing alveolar compartment. Previously, it was postulated that interstitial septal myofibroblasts were derived from the migration of peribronchiolar PDGFRα-GFP⁺/αSMA⁺ smooth muscle cells to the distally developing alveolar compartment (Lindahl et al., 1997). While migration of proximal PDGFRα⁺ cells cannot be excluded, recent linage tracing using a constitutive Pdgfra^Cre demonstrated that proximal smooth muscles cells do not contribute to PDGFRα-GFP⁺/αSMA⁺ alveolar interstitial myofibroblasts (Ntokou et al., 2015). Together with our RNA-Seq profiling data, the notion that E16.5 PDGFRα-GFP⁺/αSMA⁺ fibroblasts do not migrate from their proximal position to the alveolar compartment and rather fully differentiate into proximal PDGFRα^neg smooth muscle cells is probable. While further linage tracing experiments and or co-expression of markers will be necessary to fully answer the question, we speculate that distal interstitial PDGFRα-GFP⁺/αSMA⁺ myofibroblasts are derived from the intermediate fibroblast population identified by LungMAP single cell RNA-Seq and that these cells undergo de novo activation of PDGFRα expression after E16.5 (Lindahl et al., 1997; Ntokou et al., 2015; Guo et al., 2016, 2015).

Contrary to the proximal association with bronchiolar epithelium at E16.5, PDGFRα-GFP⁺/αSMA⁺ fibroblasts associate with alveolar epithelial cells distally at PN7 (insert Fig. 1C, Fig. 4). These fibroblasts transiently differentiate into contractile αSMA positive myofibroblasts during active septation, but once septation is complete, αSMA expression does not persist. To identify differences between proximal E16.5 smooth muscle fibroblasts (E16.5) and distal PDGFRα-GFP⁺/αSMA⁺ myofibroblasts (PN7), we interrogated gene profiles consisting of high gene expression at the time of rapid alveolarization (PN7) and contrasted them to high expression at the time of bronchiolar smooth muscle development (E16.5). Profiles #1 and #18 identify genes with high expression at E16.5 and PN7 and predict a contractile cellular function. A list of 676 actin-related genes was obtained by Ingenuity Pathway Analysis (IPA). Venn diagram comparison of these actin-related genes with the 1908 differentially expressed genes identified 147 genes associated with muscle contractile function that are differentially expressed at PN7 compared to E16.5 and, thus, distinguish the intermediate myofibroblast program of the interstitial PDGFRα⁺ fibroblast at PN7 from the E16.5 peribronchiolar smooth muscle program. (Fig. 6K, Data in Brief Table 4). Persistence and overabundance of myofibroblasts in the alveolar interstitium, because of premature birth or perinatal hyperoxia exposure, is found in association with arrested alveolar septation in bronchopulmonary dysplasia (BPD) (Ahlfeld et al., 2013; Bozyk et al., 2012; Toti et al., 1997). Further analysis of Profile #18, which contains 66 genes that are uniquely upregulated on PN7, in combination with gene expression studies from BPD samples and other septation defects will provide valuable insight into the regulatory networks that result in fibroblast differentiation to smooth muscle cells or transient myofibroblasts.

Profile #1 contains 459 genes, including 75 transcription factors, which are up-regulated at E16.5 and PN7 (Fig. 6D). The up-regulated genes comprise unique enriched functions in cell cycle, cell projection and muscle cell differentiation (Table 1). One hundred fourteen out of the 1908 differentially expressed genes were associated with cell cycle and their expression was up at E16.5 and PN7 (Data in Brief Table 2). Profile #18 contains 66 genes, including 6 transcription factors, which are uniquely up-regulated at PN7. The unique gene enrichment categories of profile #18 are cell differentiation, tube morphogenesis, response to TGFβ signaling, and muscle contraction (Fig. 6E).

We compiled a list of 288 PDGFR signaling pathway genes by combining GO categories for “Signaling by PDGF”, “PDGFR interactions”, “PDGFR signaling path” and “PDGFR binding”. Out of these 288 genes 48 were significantly differentially expressed at one of the time points (Fig. 6L).

2.13. Transcription factor promoter enrichment analysis

To predict potential transcriptional regulators of the gene profiles that we identified, we used two complementary computational approaches. For each profile, we identified all genes encoding transcription factors (TFs), and assessed the enrichment of their (1) predicted binding sites or (2) ChIP-Seq peaks, in the promoter regions of all genes contained within that profile (see Methods). The first approach identified 9 enriched TF binding site motifs in profiles #1 and #13, 2 enriched motifs in profile #10 and 1 motif in #18, #23 and #39 (Data in Brief Table 5). Several of these TFs are involved in adipogenesis, fibrosis, tumor progression, smooth muscle differentiation, or cell-matrix adhesion. However, a role in lung development has only been demonstrated for NFIA and NFIX (Hsu et al., 2011).

The second approach identified three TFs whose ChIP-Seq peaks significantly overlap the promoters of the genes contained in a particular profile: RE1-silencing transcription factor (REST/NRSF) in profile 1, CCCTC-binding factor (CTCF) in profile 1, and myc-associated factor X (MAX) in profile 13 (Data in Brief Table 6 and 7). Of note, CTCF and MAX were also implicated in the first computational approach.

The NRSF ChIP-Seq dataset that overlaps 207 of the 465 gene promoters in profile 1 has been performed in murine C2C12 myoblasts, that were differentiated into myocytes in vitro (GEO:GSM915175). Profile #1 has the two peak expression levels at E16.5 and PN7, which are the time points of smooth muscle and interstitial myofibroblast differentiation. While PDGFRα⁺ fibroblasts are not C2C12 myoblasts, our analysis infers that both cell types use a similar gene expression program for muscle and myofibroblast differentiation.

The CTCF ChIP-Seq dataset that overlaps with 286 of the 465 gene promoters in profile 1 has been performed in murine erythroleukemia cells (MEL) (GEO:GSM918744). MEL have been extensively studied for their ability to specifically bind fibronectin, and lose adhesion during cellular differentiation (Fraser and Gordon, 1994; Patel and Lodish, 1986). PDGFRα⁺ fibroblasts produce and bind fibronectin and change fibroblast phenotypes dependent on matrix adherence. An overlap with 286 genes of the CTCF ChIP-Seq dataset in MEL cells with the expression of genes in PDGFRα⁺ fibroblasts imply a cellular function for these genes in fibronectin interaction and cell adhesion.

Among the 207 NRSF and 286 CTCF genes, 144 genes were common (Data in Brief Table 7). Gene ontology analysis using Toppcluster identified “cell cycle” (32 genes), “non-canonical Wnt signaling” (10 genes), and “cell junctions” (29 genes) as significantly enriched biological processes and cellular compartments among these shared genes. In the 142 CTCF-specific genes, “mitotic cell cycle” (25 genes) were identified, and in the 63 NRSF-specific genes “regulation of cell morphogenesis” (11 genes) and “regulation of anatomical structure morphogenesis” (14 genes) were the most strongly enriched. “Interaction with fibronectin” was another strongly enriched term in both the CTCF-bound gene list (21 genes) and NRSF-bound gene list (21 genes).

The MAX ChIP-seq experiment has been performed in murine B-cell lymphoma cells (GEO:GSM912908). 192 of the 257 profile #13 gene promoters overlap with this dataset (Data in Brief Table 7). Gene ontology analysis using ToppFun identified the ER as the major cellular component (59 genes); “ER to Golgi vesicle-mediated transport” as one of the most significant biological process (18 out of 178 genes in the annotation), and “protein processing in ER” (18/167), “fatty acid elongation” (6/23) and “lysosome” (10/122) as major KEGG pathways.

2.14. Temporal changes in the composition of PDGFRα fibroblast population

To complement conclusions from gene expression data about functional differences and proliferation in PDGFRα⁺ cell populations at different stages throughout lung development, we performed flow cytometry of lung single-cell suspensions at E16.5, E18.5, PN7 and PN28 (Figs. 7 and 8). Contrary to previous approaches that have utilized differential adherence to enrich for lung fibroblasts prior to flow cytometry analysis (Chen et al., 2012; Kimani et al., 2009; McGowan and McCoy, 2011), we utilized a whole-lung single-cell suspension and gated for CD45⁻, CD326⁻, CD31⁻, PDGFRα⁺(CD140a⁺) fibroblasts (Fig. 7A). Compared to a whole-lung approach, isolating fibroblasts by differential adherence results in a 10% loss of all PDGFRα-GFP⁺ fibroblasts and reduces by 13% the relative contribution of PDGFRα-GFP with high GFP expression (Data in Brief Fig. 2).

Temporal changes in the overall composition of the single cell digest were assessed in PDGFRα⁺(CD140a⁺) and CD140a^neg stromal cells, leukocytes (CD45⁺), epithelial cells (CD326⁺), and endothelial cells (CD31⁺) (Fig. 7B). Consistent with our PDGFRα-GFP immunofluorescence data, the contribution of PDGFRα⁺(CD140a⁺) positive cells to total lung cells was more than 50% at E16.5, declined substantially to less than 10% by E18.5 and PN7, and contributed ~5% of total lung cells at PN21 and PN28 (Fig. 7B, E–H). Flow cytometry analysis of analysis of single cell suspensions from Pdgfra^GFP lungs corroborated these findings (Fig. 1H). CD140^neg stromal cells peaked at E18.5, at which time they represented ~30% of total lung cells, but were around 10% at all other time points. Leukocytes and epithelial cells comprised 30–50% of all total lung cells from E18.5 to PN28. Endothelial cells comprised ~10% of the total lung cell population throughout lung development (Fig. 7A and B).

The PDGFRα⁺ fibroblast population was further analyzed for composition of myo- and lipofibroblasts by using flow cytometry for CD140a⁺, αSMA and lipidTOX (Fig. 7C, D). At PN7 over 75% of PDGFRα⁺ fibroblasts were myofibroblasts (CD140a⁺/αSMA⁺), which is in accordance with immunofluorescence staining and the enrichment for expression of genes with contractile function (Figs. 1C, 4, 6K). Although not discernible at other time points, at PN7 PDGFRα-GFP⁺ fibroblasts could be differentiated in PDGFRα-GFP^dim and PDGFRα-GFP^bright subpopulations. Notably, the majority of lipidTOX⁺ lipofibroblasts were PDGFRα-GFP^dim and αSMA⁺ myofibroblasts were PDGFRα-GFP^bright (Data in Brief Fig. 3).

PDGFRα⁺(CD140a⁺/lipidTOX⁺) lipofibroblasts were most abundant at PN7 but were only a minority fibroblast population when compared to the size of the myofibroblast population, which might explain why a lipofibroblast signature at PN7 was not readily detectable (Fig. 7D, G). To identify lipofibroblast distribution, we performed cytohistochemistry for Oil-red-O to assess lipid droplets in PDGFRα-GFP⁺ cells and immunohistochemistry for ADFP staining (Figs. 4D, 5D, 7E–H). Based on flow cytometry not all myo- and lipofibroblasts were PDGFRα⁺(CD140a⁺), as αSMA and lipidTOX were also found in CD140^neg stromal cells (Data in Brief Fig. 4E). These findings are in accordance with the gene expression data, which indicated postnatal lipofibroblast differentiation.

To validate the prediction of increased cell-cycle activity at E16.5 and PN7, we assessed cell proliferation by flow cytometry using Ki67 expression. From E16.5 until PN7, PDGFRα⁺ fibroblasts demonstrated high rates of proliferative activity ( > 50%), which is subsequently reduced dramatically by PN21 and PN28 (Fig. 7D, Data in Brief Fig. 4D). These findings are in accordance with previously published data on fibroblasts isolated by differential adherence (Kimani et al., 2009).

Proliferation was also assessed in non-stromal cell populations. Peak proliferation of CD45⁺ cells occurred between E16.5 and E18.5, while peak proliferation in CD326⁺ epithelial and CD31⁺ endothelial cells occurred at E16.5 (Data in Brief Fig. 4A–C). These data demonstrate that different cell compartments proliferate at different developmental stages, which adds an additional layer of complexity to gene expression analysis from whole lungs versus preparations enriched for specific cell populations.

Expression of CD90 (Thy-1), CD73 (Nt5e), CD105 (Endoglin), CD44, and CD146 (Mcam), and have previously been described in the PDGFRa-GFP expressing fibroblast lineage (Ntokou et al., 2015). Relative mRNA expression of these potential progenitor markers was extracted from our RNA-Seq expression analyses and demonstrate a dynamically changing expression pattern as previously described. Relative mRNA expression of CD34 and CD29 (Itgb1) was two to four-fold higher than expression of these previously-described fibroblast markers and was significantly different at the peak of interstitial myofibroblast differentiation at PN7 (Fig. 7I, J).

2.15. Immunophenotyping PDGFRα⁺ fibroblasts using CD34 and CD29

Although PDGFRα⁺ myo- and lipo-fibroblasts can be identified by staining with αSMA and lipidTOX, respectively, there is a population of interstitial PDGFRα⁺ fibroblasts that are neither myo- nor lipofibroblasts. In a re-alveolarization model following partial pneumonectomy, we previously identified αSMA/lipidTOX-negative matrix fibroblasts that express CD34. In contrast, expression of CD29 (β1 integrin) identified proliferative myofibroblasts. Our previously published realveolarization data suggested that PDGFRα⁺/CD29⁺/CD34⁺ fibroblasts contribute to the transient re-emergence of αSMA⁺ myofibroblasts as they lose CD34 expression (Green et al., 2016). We performed flow cytometry to determine whether similar functional differences within PDGFRα⁺(CD140a⁺) fibroblasts are associated with CD29 and CD34 during lung development. Our comparison of CD140a⁺ fibroblasts isolated by lung digest versus differential plate-down revealed a loss of CD29+ cells due to trypsinization after in vitro culture (Brown et al., 2007) (Data in Brief Fig. 4).

At E16.5 and E18.5, the majority of PDGFRα⁺(CD140a⁺) fibroblasts are CD29⁺(Fig. 8A). Starting at PN7 there is a shift toward CD29⁺/34⁺ co-expression followed by a loss of both CD29 and CD34 by PN28 (Fig. 8A). To determine whether CD29 or CD34 are associated with proliferation, αSMA or lipid we analyzed CD140a⁺/Ki67⁺, CD140a⁺/αSMA⁺ and CD140a⁺/lipidTOX⁺ cells for the presence of CD29 and CD34 (Fig. 8B–D). Prenatally, the majority of CD140a⁺ fibroblasts express CD29. At PN7 and PN21 the percentage of cells that express both CD29 and CD34 increases to ~20%. The percentage of cells that express neither CD29 nor CD34 also increases from PN7 to PN28. At both E16.5 and PN7, the majority of proliferative CD140a⁺ cells express CD29, but neither CD29 nor CD34 is expressed by PN21 and PN28 (Fig. 8B). The majority of CD140a⁺/αSMA⁺ cells express CD29 at all time points (Fig. 8C). CD140a⁺/lipidTOX⁺ cells were only found at the postnatal time points (Fig. 7G, H). At PN7, CD140a⁺/lipidTOX⁺ cells express CD29, at PN21 they co-express both CD29 and CD34, and at PN28 they express neither (Fig. 8D). These data suggest that both peribronchiolar and alveolar myofibroblasts are proliferative and express CD29. Although there was no evidence that myo- and lipofibroblasts express CD34 alone, at PN21 there is a significant peak of CD29/CD34 co-expression in lipofibroblasts.

4.1. Materials and methods

Additional detail on the methods is provided in the Data in Brief.

4.2. Confocal microscopy

Lung tissues were harvested, fixed with 4% PFA in PBS and frozen, or paraffin embedded, sectioned and further processed. Tissue preparations were blocked with donkey or goat serum and incubated with anti-αSMA (Sigma-Aldrich, St. Louis, MO), SOX2, SOX9, TTF1, Endomucin, E-cadherin, HOPX (Santa Cruz, CA, USA), Pro-SPC and ADFP (Abcam, Cambridge, MA) followed by incubation with secondary antibodies conjugated to Alexa Fluor-488, 595 or 647 (Invitrogen, Carlsbad, CA) and counterstained with DAPI-containing mounting media. GFP signals were detected using a chicken polyclonal anti-GFP antibody (Abcam, Cambridge, MA). Data was analyzed by Imaris software, version 7.6.

4.3. Transgenic mice

B6.129S4-PDGFRα^{tm11(EGFP)Sor}/J herein designated Pdgfra^GFP (Hamilton et al., 2003): The endogenous PDGFRα promoter drives expression of the H2B-eGFP fusion gene.

4.4. PDGFRα⁺ fibroblast sorting, RNA extraction and RNA-Seq

Lung tissue was harvested and digested into single cell suspension. PDGFRα⁺ fibroblasts were isolated, using microbeads (Miltenyi Biotec MACS technology, Gladbach, Germany), by negative selection for CD45, CD326, and CD31 (to deplete leukocytes, epithelia, and endothelial cells, respectively) followed by positive selection for CD140a. Cells of 6–10 pubs per time point were pooled. RNA was extracted from MACS-sorted PDGFRα⁺(CD140α) cells using RNeasy Mini Kit (No. 74104) (QIAGEN, Germantown, MD). RNA was quality checked, transcribed into complementary DNA using the Verso complementary DNA synthesis kit (AB-1453; Life Technologies, Carlsbad, CA), sheared, amplified, adaptor ligated and sequenced for single-end and pair-end RNA-Seq reads (Illumina Inc. San Diego, CA, USA).

4.5. Bioinformatics data analysis

RNA-Seq data was quantitated using TopHat and Cufflinks (Trapnell et al., 2012), normalized using trimmed means and genes were pre-filtered out if the expression level (FPKM) was less than 1 in all samples. Normalized data were subject to Bayesian Analysis of Time Series (BATS) using a normal error distribution model with a cutoff of Bayesian Factor < 0.01 for differentially expressed genes (Angelini et al., 2008). Differentially expressed genes were subject to pattern recognition using STEM (Ernst and Bar-Joseph, 2006). Gene set enrichment analysis was carried out using Toppgene (http://toppgene.cchmc.org) (Chen et al., 2009). Unique functional enrichments within each profile were identified and all profiles were compared to each other using Toppcluster (http://toppcluster.cchmc.org/) (Kaimal et al., 2010). Bayesian-based analytic tools that measure time-dependent differential gene expression require at least 5 time points - therefore three additional time points (E17, PN0 and PN18) were interpolated using the R script smooth.spline function (Bioconductor 3.3).

4.6. Transcription factor promoter enrichment analysis

We used two complementary computational methods to identify transcription factors (TFs) that might regulate the genes contained in each gene expression profile. The first method performs TF binding site motif enrichment analysis on the promoters of each profile. Promoters were defined as −200 to +100 relative to the transcription start site. To this end, we applied the HOMER motif enrichment algorithm (Heinz et al., 2010) and a large library of mouse position weight matrix binding models obtained from the CisBP database (Weirauch et al., 2014).

The second method compares the promoter regions of a profile with a large library of experimentally determined TF-genome interactions. We created the library by compiling 498 mouse ChIP-seq datasets from a variety of sources, including ENCODE (Consortium, 2012), Cistrome (Liu et al., 2011), and PAZAR (Portales-Casamar et al., 2009). As input, our method takes a set of genomic regions of interest (in this case, the promoter regions of the genes in a given profile) and systematically overlaps them with each ChIP-Seq dataset. The observed overlap of the input set with each ChIP-Seq dataset is calculated by counting the number of input regions it overlaps by at least one base. Next, a P-value describing the significance of this overlap is estimated using a simulation-based procedure. A distribution of expected overlap values is created from 1000 iterations of randomly choosing RefSeq gene promoters with the same length as the input set (e.g., if 50 promoters of length 300 bp are used as input, then 50 randomly chosen promoters of length 300 bp will be used in each simulation). The distribution of the expected overlap values from the randomized data resembles a normal distribution, and is used to generate a Z-score and P-value estimating the significance of the observed number of input regions that overlap each ChIP-Seq dataset. The resulting data thus provide ranked lists of datasets, based on enrichment for experimentally determined ChIP “peaks” located in the promoters of each gene set. From the results of our analysis all of these ChIP experiments were performed in mouse cell lines.