Temporal dynamics of the developing lung transcriptome in three common inbred strains of laboratory mice reveals multiple stages of postnatal alveolar development

Tissue and RNA

Embryos and mice for this study were collected from timed pregnant mice for C57BL/6J (B6; stock # 000664), A/J (AJ; stock # 000646), and C3H/HeJ (C3H; stock # 000659) inbred lines from The Jackson Laboratory. Whole lungs were surgically dissected at embryonic (E) days E11.5, E12.5, E13.5, E14.5, E15.5, E16.5, E17.5, E18.5, E19.5 (A/J only) and postnatal (P) days P0, P2, P3, P4, P5, P7, P9, P11, P12, P13, P14, P18, P21, P24, P30, and P56 (Fig. 1). Dissected tissues were placed immediately in RNALater (Ambion) and stored at -80°C. Due to the difficulty of manually dissecting lung buds at E9.5, whole embryos were harvested for RNA extraction at this time point. For time points E11.5 and E12.5, tissues were pooled from multiple embryos (four and three, respectively) to obtain sufficient RNA for array analysis. For all subsequent time points, tissues from at least three animals were collected and analyzed separately. As the A/J strain has an extra day of gestation relative to the other inbred strains used in this study, embryos were collected at E19.5 only for this strain. Only male mice were used in this study. The sex of every embryo and pup less than 21 days of age was determined by PCR of DNA extracted from tail fragments using primers designed to amplify the male-specific Sry gene (MGI:98660) (Data S1). All animal work was performed in accordance with Jackson Laboratory Animal Care and Use Committee (ACUC) protocol 101011 (CJB).

Total RNA was extracted using mirVana RNA isolation kits (Ambion) following manufacturer’s directions. An additional DNase I (Promega) digestion step was added before organic extraction. The quality of each RNA sample was assessed using RNA 6000 Nano LabChip assays and a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). RNA purity (OD260/280 and OD 260/230 ratios) was assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific); only samples with ratios above 2.0 were included in this study.

Data processing

Gene expression was assayed using the Affymetrix Mouse Gene 1.0 ST array platform (Affymetrix). The .CEL files from the arrays were deposited in the Gene Expression Omnibus (GEO) archive at NCBI (Barrett et al., 2013) (Accession GSE74243). The array data were processed using the DeSNPer pipeline developed at the Center for Genome Dynamics (CGD) at The Jackson Laboratory (https://github.com/jaxcs/DeSNP). The DeSNPer algorithm removed probes with single nucleotide polymorphisms (SNPs) (Yalcin et al., 2012) among the three strains from further analysis to ensure that gene expression differences detected among the strains were not due to hybridization biases resulting from strain-specific sequence variation. Following probe removal, the data were log2-transformed, then probe-level hybridization intensities were quantile normalized and median polished to generate summarized gene-level expression values (Irizarry et al., 2003). A Z-score transformation was then applied, standardizing each gene to an average of zero and standard deviation of one across all samples. Finally, to focus data interpretation on the dominant patterns of gene expression changes across the developmental time series, the dataset was filtered to retain the top 30% genes with the largest pre-Z-transform standard deviation (referred to as the “filtered dataset”; Data S2). The temporal expression patterns for all of the genes measured on the arrays can be visualized using the web interface implemented for this project (http://lungdevelopment.jax.org/).

Principal Component Analysis (PCA)

The filtered dataset described above (Data S2) was analyzed using Principal Component Analysis (PCA), a data transformation method to reduce the dimensionality of a dataset to those orthogonal components that account for the maximum variance in the data with the fewest number of observations. PCA was performed using JMP^® v11.1.1 (SAS Institute) using the principal component toolset within the multivariate library. All default settings were used and no weighted bias was applied. PCA was first used to identify any systematic biases between samples, specifically among batches of arrays that were processed together (i.e., batch effects). This process identified a subset of 42 arrays that were strongly correlated with each other independent of strain or time point. The batches containing these arrays corresponded with changes in reagent chemistry used for array hybridization and were removed from further analysis. PCA was then performed on the final set of 216 arrays to identify global patterns of expression variation between samples. This analysis generated two data matrices: gene loading values (Data S3), which represent the correlation between each gene and the expression patterns captured by each principal component (PC), and sample scores (Data S4), which reflect the relationship between each sample and each principal component. The gene loading values are either positive or negative, with positive (negative) values reflecting correlation (anti-correlation) between the gene expression and the sample score plot (Fig. S1). Characteristic gene sets for each component were created using the most highly correlated genes (positively and negatively) from each PC.

Regression analysis

Two standard least squares regression analyses were performed in JMP^® (v11.1.1) to model the strain and developmental time point effects associated with the gene expression patterns identified by PCA. In the first analysis, the PCA score was modeled with independent additive effects from developmental time point and mouse strain: P_ij = a_i + b_j + e_ij, where P_ij is PCA score, a_i is the effect of time point i, b_j is the effect of strain j, and e_ij is the residual for strain j at time i. There were not sufficient data to model a strain-by-time point interaction term. Multiple hypothesis correction was implemented by controlling for FDR < 0.1 (Benjamini & Hochberg, 1995). The computed time point effects were plotted by component and used to identify contiguous bins of time points with similar PCA scores. These bins, in turn, were used to define gene expression-based stages of lung development. In cases of overlapping blocks between different block endpoints, a visual inspection was used to manually select the best grouping (e.g., postnatal time point 3 was grouped with P0 and P2, but could alternatively have been included with P4 and P5).

In the second analysis, PCA score was modeled as function of strain, stage, and the interaction of strain and stage using the stages identified in the first regression analysis: P_ij = d_i + b_j + g_ij + e_ij, where P_ij is PCA score, d_i is the effect of developmental stage i, b_j is the effect of strain j, e_ij is the interaction term between stage i and strain j, and e_ij is the residual for strain j at time i. Effect significance was corrected for by multiple hypothesis testing at FDR < 0.1, and pairwise comparisons between effects (stages, strains, and stage*strain) were carried out with the Tukey-Kramer Honestly Significant Difference (HSD) test (Kramer, 1956). The second regression analysis was performed to estimate effects for each strain (presumed to be constant across all time points), effects for each stage (presumed to be constant across all strains), and effects for each strain at each stage, allowing for the identification of stage-specific differences among strains for each PCA component.

Genes from the filtered dataset (Data S2) were binned based upon their patterns of expression as a function of stage and/or strain to identify strain-dependent differential expression. Standard linear pairwise correlation analysis (Pearson correlation) was also used to quantify relationships between strains from each PCA score (Table 1).

Differential gene expression analysis: Significance of Microarrays (SAM)

Differential gene expression between stages of development was determined using Significance Analysis of Microarrays (Tusher, Tibshirani & Chu, 2001) with a median FDR = 0 using the TM4 package (Saeed et al., 2003). Binning of individual time points was based on the developmental stages defined by the regression analysis described above.

Quantitative Real-Time PCR (qPCR) validation of strain-specific gene expression

We used quantitative Real-Time PCR (qPCR) to validate strain-specific expression of selected genes. cDNA was generated from nine samples of total RNA originally used for microarray analysis, using RT² First Strand Kit (Qiagen) following manufacturer’s directions. cDNA was generated from three biological replicates for each strain at postnatal day 13 (P13). This time point was chosen, in part, because the strain-specific expression differences of the genes selected for validation (see below) were evident at this time and there was sufficient quantities of the same high quality RNA used for the microarray analysis to use for qPCR. Gene expression was assayed using RT² SYBR Green ROX qPCR MasterMix Kit (Qiagen) following manufacturer’s directions. Three genes with low variation in expression by time point or strain (Actb, Rpl10, Rpl13a) were used to normalize expression values of the target genes. A geometric mean of the control genes was determined for each biological triplicate to generate a normalization value. Following normalization, relative expression values were calculated for each target gene (Fggy, Saa3, Wnt11, Wif1) within all biological groups. The genes for validation were selected because of their relevance to fibrosis-related biological processes that are known to differ among inbred strains of mice (Walkin et al., 2013). The significance of differences between strains in relative expression arbitrary units (A.U.) was assessed by Tukey Multiple testing corrected ANOVA using Prism v7.0a (GraphPad Software). Primer sequences used for qPCR were as follows: Actb, 5′-GTGACGTTGACATCCGTAAAGA-3′, 5′-GCCGGACTCATCGTACTCC-3′; Rpl10, 5′-CGTGGTGTCCCTGATGCTAAG-3′, 5′-GTTGGCACAAATACGGGCAG-3′; Rpl13a, 5′-AGCCTACCAGAAAGTTTGCTTAC-3′, 5′-GCTTCTTCTTCCGATAGTGCATC-3′; Fggy, 5′-CGGCTCTAGTGGACCAGAGA-3′, 5′-ATGTGCATCAATCCCTTGAACA-3′; Saa3, 5′-AGAGAGGCTGTTCAGAAGTTCA-3′, 5′-AGCAGGTCGGAAGTGGTTG-3′; Wnt11, 5′- ATGCGTCTACACAACAGTGAAG-3′, 5′-GTAGCGGGTCTTGAGGTCAG-3′; Wif1, 5′-CCTGCCGAAATGGAGGTAAA-3′, 5′-GCTGGCTCCATACCTCTTATTG-3′.

Annotation term enrichment analysis

Annotation enrichment analysis was performed to guide the interpretation of gene sets associated with specific developmental stages and for those contributing to individual principal components (PCs). Probe set identifiers from the Affymetrix arrays were first mapped to mouse gene identifiers using Mouse Genome Informatics (MGI) database Batch Query tool (Bult et al., 2010). Sets of mouse gene symbols or gene accession identifiers were submitted for Gene Ontology (2015) and Mammalian Phenotype Ontology (Smith & Eppig, 2012) term enrichment analysis using VisuaL Annotation Display (VLAD) v1.5.1 (Richardson & Bult, 2015). Mouse genes that were converted to human homologs were also analyzed using the Reactome Pathway Browser v52 (Croft et al., 2014) to identify enrichment of homologous pathway annotations.

Mouse and human gene homology

Gene homology between mouse and human was used to compare transcriptional profile results from this study to those from a previous study of human lung development (Kho et al., 2010). A gene homology resource containing 10,266 genes was generated manually using homology assertions from Panther (Mi, Muruganujan & Thomas, 2013) and Homologene (Sayers et al., 2012). If one gene mapped to multiple homologs, each homolog was represented in the final set. For example, the mouse alcohol dehydrogenase 1 (Adh1) gene is homologous to ADH1A, ADH1B, and ADH1C in human; the carboxylesterase (CES) gene in human is homologous to Ces1d, Ces1e, and Ces1g in mouse. Many-to-one homolog relationships were identified in 18 instances. Overlap between lists of homologous genes was assessed with the GeneVenn tool (Pirooznia, Nagarajan & Deng, 2007) using the NCBI GeneID for the human gene as the primary gene accession identifier.

Principal component and regression analysis reveals nine stages of murine lung development

The first three principal components (PCs) accounted for 66.2% of the variation among samples; the first ten PCs accounted for 80.9% of the sample variation (Table 1). The patterns of variation in gene expression represented by PCs 1–3 represent common patterns of expression that do not very significantly between the three mouse strains (Fig. 2). In contrast, PCs 4–10 demonstrated strain-specific patterns of expression variation (Fig. S2). Regression modeling supported the PCA results (Table 1); no significance was detected among the strain or strain by stage effects for PCs 1–3 whereas PC 4–10 all were found to have differences between one or more of the strains for some of the developmental stages (Fig. 3). To identify possible temporal shifts in gene expression patterns between strains, correlations across all strain by PC combinations were performed. No significant correlations from this analysis were observed.

Regression analyses of the PCA results support the grouping of sampled time points into nine stages of lung development (Fig. 4). The four prenatal stages, embryonic (EMB, E9.5–E12.5), pseudoglandular (PSG, E13.5–E15.5), canalicular (CAN, E16.5–E17.5), and saccular (SAC, E18.5–E19.5) are concordant with those defined previously by histology and morphology. We identified four molecularly distinct stages of alveolar development between P0–P18 (ALV1-4) that are defined by the expression patterns and functional properties of differentially expressed genes. Finally, the time points following alveolarization were grouped under the common heading of mature lung (MAT, P21–P56).

Strain-independent principal components 1–3 define a murine developing lung characteristic subtranscriptome (mDLCS)

The first PC (55.1% of the sample variation) was significantly correlated (P < 0.0001) with developmental time point, capturing the patterns of gene expression across the entire developmental timeline. Over 50% of the genes in our filtered dataset (Data S2) had relatively high (positive) or low (negative) loading values on PC1. GO term enrichment analysis of genes contributing to the prenatal signal (PC1^pos) revealed enrichment of genes associated with nucleic acid metabolic process (GO:0090304) and RNA processing (GO:0006396). Genes previously associated with lung cell differentiation were among the top 10% of contributors to PC1 (Fig. S6); a 3.2-fold enrichment (Fisher exact test; P < 1.7 × 10⁻³). Annotation enrichment analysis of genes contributing to the postnatal signal (PC1^neg) identified enrichment of immune system processes (GO:0002376), cellular communication (GO:0010646), and localization (GO:0051179). Specifically, we observed postnatal induction of genes associated with RAS protein superfamily (RAS), Ras-related protein 1 (Rap1), phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt), and mitogen-activated protein kinase (MAPK) signaling as well as chemokine signaling and phagocytosis/endocytosis, including seven toll-like receptors, four lymphocyte antigens, six chemokine receptor/ligands, and eighteen interleukins (Table S1). Many surfactants (4 of 4), mucins (2 of 4), and extracellular matrix genes (6 of 10 laminins; 7 of 21 collagens; 10 of 18 integrins) also contribute significantly to PC1^neg.

The pattern of sample variation captured by PC2 revealed major changes in gene expression between the EMB/PSG and ALV3/ALV4 stages, depicted as shifts in PC2 sample coordinates between stages (Fig. 2). The shift between EMB/PSG stages represents a molecular transition from early embryonic organogenesis to prenatal epithelial/vascular development; genes that were significantly down regulated between these stages are primarily associated with embryonic morphogenesis whereas up regulated genes are associated with cell migration and angiogenesis. The shift at ALV3/ALV4 represents a transition from active cellular growth to homeostasis; genes that were down regulated between these stages are associated with cell cycle whereas genes that were up regulated are associated with metabolism (lipid, oxidation–reduction) and immune response.

Shifts in gene expression captured by PC3 occur between the PSG/CAN, ALV1/ALV2, and ALV4/MAT stages (Fig. 2). Down regulated genes between PSG/CAN are primarily associated with cell cycle and up regulated genes are associated with cellular differentiation and localization, suggesting that processes related to differentiation and morphogenesis dominate late prenatal development. The shift between ALV1/ALV2 also represents a transition from metabolism/homeostasis to induction of cell cycle; down regulated genes between these stages are associated with lipid metabolism, chemical homeostasis, and localization while up regulated genes are highly enriched for cell cycle. All differentially expressed genes between molecular stages can be found in Data S5.

Strain-specific patterns of gene expression during lung development

As reported above, PC 4–10 lacked significant correlation among the strains and regression modeling revealed significant strain effects (Table 1). Of 20 calculated strain terms, 11 were significant (FDR < 0.1), including at least one term for each of PC 4–10. In contrast, only four PCs (5–6, 9–10) had a significant interaction (strain*stage) effect for one or more developmental stages; of 180 calculated interaction terms, only 7 were judged significant (FDR < 0.1). Thus, the dominant strain-specific expression patterns captured by PC 4–10 are stage-independent, with only smaller stage-specific variations observed.

PCs 4–10 were evaluated for patterns of significant strain effects, specifically how the patterns differentiate each strain from the others. Although nineteen possible patterns of strain-dependent expression are possible with three strains, the results from regression modeling demonstrated that the following four classes account for nearly all strain-dependent effects (Fig. 3):

• C3H different from indistinguishable pair AJ/B6 (PC4)

• B6 different from indistinguishable pair AJ/C3H (PC6, PC7)

• AJ different from indistinguishable pair B6/C3H (PC8)

• All strains different (PC5, PC9, PC10)

Regression modeling of gene expression as a function of strain, and subsequent Tukey HSD testing (Q < 0.05), quantified the genes that significantly contributes to each class (Table 2). Gene expression associated with lung development is most dissimilar between B6 and AJ or C3H mice; a 40/60 ratio of genes that are up versus down regulated in B6 relative to AJ/C3H was also observed. Differences between C3H vs AJ/B6 and AJ vs C3H/B6 exhibited 40/60 and 60/40 ratios of up versus down regulated genes, respectively. Differences between all strains were less common; only 57 strain-dependent genes (7.3%) were differentially expressed between all three strains. Strain-dependent gene expression was validated by qPCR (Fig. S11) in four target genes (Fggy, Saa3, Wnt11, Wif1). The qPCR results confirmed the trends observed in the array data despite the variation observed among biological replicates (Table S2). All strain effects for Saa3 and Wif1 were confirmed by qPCR; that is, Saa3 expression is greatest in C3H and indistinguishable between AJ and B6 whereas Wif1 expression is greatest in AJ and indistinguishable between B6 and C3H. The qPCR results for Fggy were consistent with the trend of AJ < B6 < C3H; however, the expression differences were not statistically significant between AJ and B6. Similarly Wnt11 expression was greatest in B6 but that B6 was not statistically significant different from the AJ strain. These results suggest that the general trends of strain-dependent expression captured by the arrays are robust but that expression of strain-differences for individual genes should be validated by qPCR or single cell sequencing prior to experimental follow up to investigate the biological significance of these differences.

Comparative analysis of the mouse and human embryonic developing lung characteristic subtranscriptomes

To investigate the genomic elements of lung development that are conserved between human and mouse we performed PCA on a reduced set of murine samples (E12.5–E16.5) that corresponded to the time points assayed in a previously published human Developing Lung Characteristic Subtranscriptome dataset containing 3,209 genes (Kho et al., 2010). We generated a prenatal mDLCS using the union of rank-ordered genes from PC1-3 that resulted in a set of 3,077 genes. Most of the genes (2,302) in the prenatal mDLCS overlapped the mDLCS that included all 26 developmental time points.

Homologous human genes could not be identified for 464 of the prenatal mDLCS genes and homologous mouse genes could not be identified for 201 of the hDLCS genes, leaving 2,634 genes in the prenatal mDLCS for comparison to the 2997 hDLCS gene set. The genes from the prenatal mDLCS genes for which no human homolog could be found included predicted genes, microRNAs, long non-coding RNAs (lncRNAs), and species-specific immune-related proteins (e.g., major histocompatibility complexes). There were 771 genes shared in common (RF = 2.3, P < 2.4 × 10⁻¹³⁴) between the two datasets following the conversion of the mouse genes to human homologs; 2,226 genes were unique to the hDLCS and 1,861 genes were unique to the prenatal mDLCS (Fig. 5).

There are two primary reasons for the lack of overlap of specific genes between the mouse and human subtranscriptomes. First is the difference in the genes represented on the two gene expression Affymetrix array platforms (Mouse 1.0ST and the Human 133 Plus 2.0). Mouse homologs for 68 genes included in the hDLCS were not assayed by the Mouse 1.0ST array. Second, mouse homologs for 1,216 of the genes in the hDLCS lacked variance in expression during embryonic development and were removed at the variance-filtering step prior to PCA. These genes therefore, were not included in the prenatal mouse DLCS.

Although the overlap of individual genes between the hDLCS and mDLCS was small, the biological processes and pathways represented by the DLCS gene lists were similar. Both DLCS gene sets were enriched in genes involved in broad developmental processes associated with lung development including cell cycle, lung-specific metabolism, signal transduction, and a wide range of immune functions (Table 3).

The 2,226 genes unique to the hDLCS set were significantly enriched for cell cycle and DNA repair processes, which may reflect differences in tissue quality, harvesting, and/or processing between studies. The 1,861 genes represented only in the mDLCS were enriched for high-level biological processes (cell cycle, DNA replication) as well as ECM organization, specifically non-integrin membrane-ECM interactions (LAMA3/4, LAMB1, COL4A1/4A2, COL2A1), ECM proteoglycans, degradation of the ECM, and integrin cell surface interactions. The set of genes associated with enrichment of ECM organization includes alpha/beta integrins (ITGA6/8, ITGB6), laminins (LAMA3/4, LAMB1), collagens (COL1A2, COL2A1, COL4A1/2, COL6A1/2/3, COL14A1, COL26A1), adhesion molecules (ICAM1/2, JAM2/3, VCAM1), cathepsins (CTSB, CTSS, CTSV), and metallopeptidases (ADAMTS5/9, MMP14/19). These results suggest that underlying differences in ECM remodeling and/or composition may exist between the mouse and human prenatal developing lung microenvironment.

The previous human lung transcriptome study by Kho et al. (2010) revealed evidence of a novel pseudoglandular substage between the 13th–17th weeks of human lung development (corresponding to E15.5 in mouse). Our mouse embryonic transcriptome data did not recapitulate this finding. Although the plots of PC sample scores for PC1 and PC2 in the B6 mouse strain show some similarity to the plots for the human data, strain-dependent variance of the PCA sample scores and the variance in gene expression for the B6 E15.5 time points complicate the comparisons. Only three of the more than 50 genes reported to be differentially expressed between the early and late pseudoglandular (PSG) time points in humans are differentially expressed between E14.5 and E15.5 in mouse: Sulf1, Muc1, and Sftpc. The lack of concordance between the mouse and human data is likely related to the sampling of lung development time points in mouse. The PSG sub stage defined in humans lasts approximately five weeks whereas the comparable developmental time in mice occurs over several hours. To definitively address the presence of a novel pseudoglandular stage in mouse would require sampling lung development at a much finer temporal scale (hourly) than was performed for the current study (one sample per day). An alternative explanation to the lack of a PSG substage in our murine dataset may lie in the fact that only male mice were used in this study; recent work in humans has revealed gender-based differences in transcriptional modulation surrounding this substage (Kho et al., 2015).

Principal components 1–3 suggest that lung development utilizes superimposed periodic patterns of transcriptional control

The plots of PC sample scores for the first three principal components across the sampled developmental time points revealed distinct periodic patterns of gene expression during lung development (Fig. 2) similar to those reported recently for developmental gene expression in the nematode, C. elegans (Hendriks et al., 2014). The temporal expression pattern for PC1 divides the developmental timeline in two segments: (1) embryonic, pseudoglandular, and canalicular stages versus (2) saccular, alveolar, and mature stages. The global trend is partially reversed in the second and third alveolar stages. The singular transition point from embryonic to post-natal pattern that occurs on PC1 between the canalicular and saccular stages supports the previous assertion that the saccular stage of lung development is a critical period of preparation for the switch to breathing air in mice (Kho et al., 2009). PC2 captures a developmental pattern of gene expression that has two major changes, occurring between EMB and PSG stages and then reversing again between ALV3 and ALV4 stages. As previously observed in Kho et al. (2009), PC2 effectively distinguishes stages temporally nearer to birth (PSG, CAN, SAC, ALV1, ALV2, and ALV3) from those more distant (EMB, ALV4, and MAT). Finally, PC3 captures a pattern with four distinct phases, with transitions found between PSG and CAN, SAC and ALV1, and ALV4 and MAT stages of development. In effect, the PC3 expression profile from early embryo to birth is repeated from birth to maturity. Somewhat strikingly, where PC2’s expression profile is inverted in mature lungs compared to late embryological development, PC3’s profile matches that of late development. The regulatory control and consequences of these periodic patterns of increasing frequency are topics for future inquiry.

Alveolarization in mouse has four distinct transcriptional stages suggesting waves of vascularization and innervation

A regression analysis of the principal components identified nine distinct stages of lung development represented by the gene expression data. Five of these stages were consistent with previous definitions of lung development (embryonic, pseudoglandular, canalicular, saccular, and maturation following alveolarization). In contrast to previous studies that proposed alveolar development as an extension of the saccular stage (Bonner, Lemon & You, 2003) or as two stages (Mund, Stampanoni & Schittny, 2008; Schittny, Mund & Stampanoni, 2008) our data suggest subdivision of postnatal alveolarization into four distinct stages (Fig. 4).

Differential gene expression analysis between successive stages of development revealed two distinct waves of angiogenesis during alveolarization (Fig. 6). Multiple angiogenic regulators, including Vegfa, Adora2b, Adrb2, Lgals3, Tgfb1 and Ang, were significantly induced in the first alveolar stage (ALV1), relative to the saccular stage (SAC), and again in the fourth alveolar stage (ALV4) relative to the third (ALV3). Further support of the two-phase pattern of angiogenesis is found in the stage-specific (ALV2/ALV3) induction of genes associated with negative regulation of angiogenesis (Thbs2, Agt) and/or vascular stabilization (Angpt1, Angpt2, Serpine1, Igf1). Seven of these vascular stabilization genes were also expressed in the mature lung (MAT) suggesting that there are periods of angiogenic “rest” during postnatal maturation and adult homeostasis (Fig. S3). A similar dual-phase pattern of induction (ALV1 and ALV4-MAT) was detected among genes known to be involved in alveolarization and neurogenesis (Fig. S4 and S5), suggesting extensive genomic coordination between vascular, neurogenic, and alveolar processes in the postnatal lung development.

PCs 1–3 define a “Developing Lung Characteristic Subtranscriptome” in mouse

We combined the genes from the principal components that had no significant strain dependent effects into a murine Developing Lung Characteristic Subtranscriptome (mDLCS). The mDLCS represents the core set of genes expressed in specific temporal patterns during lung development in the mouse, independent of strain-specific genetic effects. When we compared our mDLCS to lists of genes known previously to be involved in murine lung development we determined that the majority of these genes were represented in the mDLCS. There were two primary factors that excluded previous lung development genes from our mDLCS. First, some genes did not pass the threshold we used to exclude genes from the PCA/regression analysis that defined the mDLCS. For example, 7/20 GO-annotated genes and 67/115 MP-annotated genes that were not included in the mDLCS were near the threshold criterion and would have been included with only minor changes in the threshold requirements. Second, some genes were not included because they displayed strain specific expression patterns and therefore were captured in PC4–10 instead of PC1-3.

Strain specific gene expression during lung development

The inclusion of multiple inbred strains of mice in our study gave us the opportunity to identify strain-dependent patterns of gene expression. Regression modeling revealed the majority of these patterns represent global differences between strains that are independent of developmental stage. The four primary patterns of strain-dependent expression that we uncovered were largely characterized by expression differences in genes associated with ECM composition, lung development, neurogenesis, immune system chemotaxis and function, and profibrotic pathway signaling. Differences in leukocytic composition between strains may explain the strain-dependent differential expression of chemotactic factors; however further investigation is necessary to shed light on putative strain-dependent differences in pulmonary ECM composition, alveolar versus bronchiolar development, and pulmonary innervation.

The expression profiles of several immune-related or profibrotic genes, including those that influence inflammatory response, ECM remodeling through MMPs, and epithelial-mesenchymal crosstalk, were among the genes showing strain differences and may help explain the known differences in susceptibility for asthma and pulmonary fibrosis between B6 (susceptible) and AJ (resistant) or C3H (very resistant). Further investigation of the profibrotic Wnt signaling pathway also identified gene expression differences between B6, AJ, and C3H mice that may explain phenotypic differences in fibrotic susceptibility between strains. We observed increased expression of mast-cell secretory factors (Cpa3, Mcpt4) in AJ and C3H mice compared to B6 immediately after birth, which may reflect differences in pulmonary inflammatory response between strains and is consistent with experimental data showing elevated mast cell numbers and exacerbated alveolitis in AJ and C3H relative to B6 following radiation (Haston, 2012).

Comparison of mouse and human embryonic lung development

Our comparison of a murine embryonic DLCS to the previously published human embryonic DLCS revealed only moderate overlap at the gene level. Among the genes we identified in common between the mDLCS and hDLCS were key cell differentiation factors (AGER, ABCA3, GPRC5A, EGFL6), primary surfactant-/mucin-associated proteins (ADGRF5, CTSH, MUC1, SRFP1/2, SFTPA1/2, SFTPB, SFTPC), or major structural components, including adhesive junction-associated proteins (CDH1, CADM1, MATR3, MLLT4, SCAPER) and tight junction-associated claudins (CLDN5, CLDN18) (Table S2). The lack of substantial overlap of the DLCS gene lists at the level of individual genes may be due, in part, to the differences in the array platforms used in the two studies. Mouse lung expression was assayed using Affymetrix 1.0 ST gene arrays where probesets are largely targeted to protein coding regions. In sharp contrast, the human lung expression dataset used U133 Plus 2.0 arrays, which has probesets targeted primarily to 3′untranslated regions (UTRs). Prior studies have frequently identified developmental and cell-specific transcript isoform variation specifically in the 3′-UTR. The two array designs will have different sensitivity to such variations, potentially resulting in very different measurements, even with common transcript expression profiles, and subsequently changing the identification of genes with enough variation for inclusion into downstream analysis.

In contrast to the modest overlap of individual genes, the enriched annotations associated with the genes in the mDLCS and hDLCS revealed many shared functional properties. Both subtranscriptomes were enriched in genes involved in axon guidance, specifically through L1CAM/NCAM interactions and signaling through semaphorins/Robo, suggesting pulmonary neurogenesis events during prenatal development of the lung are highly conserved between mouse and human. Several signaling transduction pathways were also enriched among conserved elements of the mDLCS and hDLCS, including signaling through insulin receptors, protein kinases, receptor tyrosine kinases, Wnt, and Rho GTPases. Developmentally associated growth factor mediated signaling was also conserved including epidermal growth factor, fibroblast growth factor receptor, platelet-derived growth factor, nerve growth factor, and stem cell factor.

Both hDLCS and prenatal mDLCS gene sets were also enriched in general immune system cytokine signaling, including interferon- and interleukin-signaling, as well as adaptive immune system, including TCR signaling, MHC-II antigen presentation, and PD-1 signaling. This was surprising since both subtranscriptomes are based entirely on pre-birth development time points and suggests an unexpected role for adaptive immune system elements in prenatal lung development. Notable differences in annotation enrichment results between the mouse and human subtranscriptomes were also detected. For example the mouse embryonic DLCS gene set was enriched for genes involved in ECM organization whereas the human lung DLCS (but not murine) showed enrichment of genes involved in protein biosynthesis and protein modification.

The data reported here expand the current knowledge of genome wide temporal changes in gene expression during the complex process of murine lung development, revealing novel stages of postnatal alveolar development. The data generated for this study complement existing transcriptomics data sets for murine lung development such as those that are available from the LungMAP resource (http://www.lungmap.net/). The normal mouse lung gene expression time course data available from LungMAP includes only two embryonic and two post-natal time points for a single strain (C57BL/6). Our data support a much finer-grained assessment of the temporal dynamics of gene expression across lung development and the ability to identify possible strain-specific differences in expression. In contrast to our study, the LungMAP resource includes some data for cell-specific expression (from sorted endothelial cell, myeloid cell, and type II pneumocytes) at different development time points. These cell-specific expression data can be integrated with our data to develop hypotheses about which cell types may be contributing to the tissue level gene expression patterns observed in our data. Although a full treatment of combining the data set reported here with LungMAP data is beyond the scope of this report, the expression data for the angiogenin (Ang) gene serves as an example of the complementary nature of these data resources. The general pattern of increased expression of Ang over development is observed in both data resources (see www.lungmap.net and lungdevelopment.jax.org). However our data reveal that the relative level of expression of Ang is much higher in C57BL/6J than in A/J or C3H/HeJ. This strain-specific difference is not possible to detect in the existing LungMAP data. The cell-specific expression data from LungMAP suggests that type II pneumocytes are contributing most to the increased expression of angiogenin between postnatal day 7 and day 28. This testable hypothesis would not be possible to formulate using our data set alone.

Strain-specific differences in gene expression patterns were identified that serve as the basis for future investigations into the genetic basis for well-known differences in adult lung physiology and respiratory disease susceptibility among inbred mouse strains. Comparison of developmental gene expression data between mouse and human highlighted conserved biological processes as well as organismal differences that may be important for using the mouse as a model for human respiratory biology and disease. The raw data files and processed data for this study can be accessed at NCBI’s Gene Expression Omnibus (GEO) data archive (GSE74243). Expression patterns across stages (and time points) of lung development for all genes measured by the Mouse Affymetrix Gene 1.0ST array can be queried, visualized, and downloaded at http://lungdevelopment.jax.org/.