Study Design

All animal experiments were performed in accordance with protocols approved by the University of Colorado Denver (UCD) Animal Care and Usage Committee. Our analysis focused upon Theiler stages (TS) 17–20.3, which are equivalent to E10.5–E12.5 of mouse embryonic development. This time-span is the most relevant to the development of human orofacial clefts, and also marks the period in which initial differentiation and morphogenesis of the craniofacial skeleton occurs [3], [10], [11]. The earliest time point (E10.5) corresponds to when all three facial prominences can first be distinguished alongside the olfactory pits (TS17). The final time point (E12.5) is when the individual prominences are about to initiate fusion to form the primary palate after which they cannot be distinguished as discrete entities (TS20.3). Facial growth and morphogenesis is rapid during this period so samples were isolated at 0.5-day intervals to capture dynamic changes in gene expression.

The inbred C57BL/6J strain (Jackson Labs) was chosen for analysis to reduce genetic variation among samples and select changes in gene expression that result more from temporal- and spatial-specific differences. Moreover, the most extensive staging classification for development of the mouse face is based on this strain [10]. Following mating and visualization of a vaginal plug, pregnant females were sacrificed between E10.5 and E12.5. Mice employed for the E11, and E12 time points were housed in a reverse light cycle room to facilitate dissections at midday for all samples. Embryos were dissected from the uterus in ice-cold phosphate buffered saline solution (PBS) and then individually placed into drops of PBS in a petri dish. Embryos were staged by examination of somite numbers and craniofacial features (Table 1) and embryos conforming to the criteria for the 0.5-day interval under analysis were pooled. Pooling was necessary to obtain sufficient RNA for screening of the microarrays and also smoothes out variations between individual samples. The alternative, using fewer embryos but then incorporating a PCR amplification step, would potentially introduce considerable skewing and bias into the analysis. Table 1 indicates the number of embryos required at the various time points to obtain sufficient RNA for each prominence to be used for one microarray screening experiment. A schematic overview of the dissection process is shown in Figure 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Isolation of the mouse facial prominences.

Graphic representation of the embryonic mouse head illustrating the facial tissues collected for analysis.

https://doi.org/10.1371/journal.pone.0008066.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Staging criteria for assignment of embryos to particular timepoints and corresponding RNA yields. TS  =  Theiler stage.

https://doi.org/10.1371/journal.pone.0008066.t001

After staging, embryos were bisected with forceps at the level of the heart and the caudal portion was discarded. Tungsten needles were employed to isolate the combined maxillary and mandibular prominences and these were then separated into their individual components. Next, the medial and lateral frontonasal prominences were removed from the remainder of the head. Briefly, an incision was made in the ectoderm overlying the prominences to produce a flap of ectoderm that was peeled back with the loosely aggregated cells of the mesenchyme still attached. Note that if the cut was made too deep and extended into the forebrain region, loosely packed mesenchymal cells were no longer visible. Samples that were potentially contaminated with forebrain tissue were discarded. Suitably pure maxillary (MxP), mandibular (MdP), and frontonasal (FNP) samples were placed in separate tubes containing RNAlater (Ambion) and stored at −20°C for subsequent pooling and processing.

RNA preparation and quality assessment

Tissue pools were dissolved in 1 ml Trizol (Invitrogen) and total RNA was extracted according to the manufacturer's instructions and resuspended in a final volume of 100 µl DEPC-treated water. Subsequently, the sample was further purified using the RNeasy Mini Kit (Qiagen, Inc.). The concentration of each total RNA sample was determined from absorbance at 260 nm (A260) and the quality of each sample was determined based on the ratio of A260 to A280. Samples with a ratio between 1.9–2.1 were considered adequately pure for cDNA synthesis. Prior to microarray analysis, individual samples were assayed using RT-PCR to ensure RNA integrity, purity, and expression of transcripts appropriate for the particular prominences. 1–2 µg total RNA was used for cDNA synthesis using random primers and Superscript III reverse transcriptase (Invitrogen) according to the manufacturer's instructions. Subsequently, 2 µl of the cDNA mixture was used for PCR in a 25 µl reaction with Accuprime Taq DNA polymerase (Invitrogen) using the following conditions: 1 cycle at 94°C for 3 min; 35 cycles at 94°C for 45 sec, 58°C for 45 sec, 68°C for 1 min; and 1 cycle at 68°C for 10 min. Specific primers (Table S1) were directed against transcripts diagnostic for the facial prominences. Representative RT-PCR results demonstrated the presence of transcripts corresponding to expression of the marker genes Tcfap2a (GeneID: 21418), Bmp2 (GeneID: 12156), Dlx2 (GeneID: 13392), and Gsc (GeneID: 14836) in the FNP and MxP prominences, as expected (Figure S1). We also analyzed the prominence tissues for transcripts that would indicate significant contamination from surrounding tissues, particularly the presence of adjacent CNS tissue within our FNP samples e.g. Zic3 (GeneID: 22773) for forebrain tissue (Figure S1). Contaminated or degraded RNA samples were discarded.

Microarray hybridization and data capture

Microarray analyses were carried out by the UCD Gene Expression Core Facility. Total RNA (2–5 µg) was converted to double-stranded cDNA (ds-cDNA) using the Superscript Choice System (Life Technologies, Inc.) and an oligo-dT primer containing a T7 RNA polymerase promoter (Genset Corp.). Subsequently, the reaction mixture was extracted with phenol:chloroform:isoamyl alcohol and the ds-cDNA was recovered by ethanol precipitation. In vitro transcription (IVT) was performed to generate biotin-labeled cRNA using an RNA Transcript Labeling Kit (Enzo. Inc.) and 3.3 µL ds-cDNA template in the presence of a mixture of biotin-labeled ribonucleotides. Biotin-labeled cRNA was then purified using an RNeasy affinity column (Qiagen) and the cRNA was fragmented to ensure optimal hybridization to the oligonucleotide array. Fragments between 35–200 bp in length were generated by incubating the cRNA at 94°C for 35 min in a fragmentation buffer. The quality of total RNA, labeled cRNA and fragmented cRNA were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Only samples that passed all tests of quantity and quality were used for hybridization to microarrays. The sample was then added to 200 µl of hybridization solution containing 100 mM MES, 1 M NaCl, and 20 mM EDTA in the presence of 0.01% Tween 20 to give a 0.05 µg/µL final concentration of fragmented cRNA. Next, the cRNA was hybridized to Affymetrix GeneChip® Mouse430 2.0 microarrays, which contain probes corresponding to 39,000 mouse transcripts, at 45°C for 16 hr in a GeneChip® Hybridization Oven 640 (Affymetrix, Santa Clara, CA.). The microarrays were washed and stained using an Affymetrix Fluidics Station 450 according to the manufacturer's standard protocol and scanned using the Affymetrix GeneChip Scanner 3000 (Affymetrix). Raw image data from array scans was processed in the GeneChip® Operating Software (Affymetrix) using the option in which a scaling factor was applied to bring the average intensity for all probes on the array to the same target intensity value (TGT) of 500. This normalization allowed samples to be compared across arrays. All data from our analysis are MIAME compliant, and the raw data (non-log₂ transformed) are available via GEO (http://www.ncbi.nlm.nih.gov/geo) with the accession number GSE7759.

In situ hybridization

We amplified genes chosen for verification analysis from total RNA by RT-PCR using primer pairs (Table S1) as described above. Subsequently, the products were cloned using Zero Blunt TOPO vector methodology (Invitrogen) and sequenced to confirm identity. Digoxigenin-UTP labeled probes were generated from linearized plasmids with either T7 or SP6 RNA polymerase (Roche). Whole-mount in situ hybridization was then performed essentially as described using staged C57BL/6J embryos [12]. Briefly, mouse embryos were fixed overnight in 4% PFA (paraformaldehyde) at 4°C, dehydrated in a graded series of methanol buffered with PBT (PBS containing 0.1% Tween 20) and stored at −70°C in 100% methanol. After rehydration, embryos were incubated overnight at 70°C in 1 ml hybridization buffer (50% formamide, 5X SSC, 50 µg/ml yeast tRNA, 1% SDS, 50 µg/ml heparin and 1 µg digoxigenin-labeled sense or anti-sense probe). After washing to remove unhybridized RNA probe, AP-conjugated (alkaline phosphatase) anti-digoxigenin antibody (Roche) was used to detect the digoxigenin-labeled RNA probe. BM purple AP substrate (Roche) was used for colorimetric detection. Images were taken using a Spot II digital camera (Diagnostic Instruments) and processed using Adobe Photoshop.

Data filtering

To improve statistical power in light of necessary corrections for multiple testing, two filters were first applied to the dataset to discard probe sets that did not vary across any of the conditions. The first filter, the median filter [13], [14], discarded probe sets whose standard deviation across all 105 samples - using the GeneChip® Operating Software (Affymetrix) Signal values - was less than, or not significantly different than, the median standard deviation, leaving 20,868 of the original 45,101 probe sets. Selection of the median variance as a cutoff is based on the assumption that fewer than half of the genes in the genome are involved in orofacial development. The second filter removed a further 634 probe sets with a Detection call of ‘Absent’ across all 105 samples. The remaining 20,234 probe sets formed the input dataset for further statistical testing.

Statistical differences of expression

Statistically significant differential expression between any two samples was tested using the limma package [15] of Bioconductor [16] on the log₂-transformed Signal data. The limma package uses an empirical Bayes approach to create a moderated t-statistic by shrinking the estimated sample variances towards a pooled estimate thereby allowing more stable inference than an ordinary t-statistic when the sample size is small. P-values were calculated for each time point among different prominences (12 in total) as well as among adjacent pairs of time points for each prominence (30 in total). For all analyses, the false discovery rate was set to 1% with the additional requirement of at least a 2-fold change to be deemed significant (setting fdr = 0.01 and lfc = 1 in the decideTests function).

Differential expression across the entire time series was characterized by the trajectory clustering algorithm [14]. In lieu of ANOVA or Kruskal-Wallis p-values used by the original algorithm, the moderated t-statistic p-values from limma were used with a threshold of 0.05 (retaining the requirement for at least 2-fold change).

Functional category analysis

Probe sets were mapped to MGI identifiers using information from the Affymetrix netaffx web interface (downloaded 18^th Dec 2007). Each MGI was then mapped to terms from: the Gene Ontology (GO) Biological Process and Molecular Function ontologies [17] using the association file available from the GO website (downloaded 11 May 2007); InterPro categories [18] using the MRK_InterPro.rpt file from the Jackson Laboratory website (downloaded 22 Aug 2006); the Mammalian Phenotype Ontology (MP) terms using the MGI_PhenoGenoMP.rpt file from the Jackson Laboratory website (downloaded 10 Oct 2006) and; the Kyoto Encyclopedia of Genes and Genomes (KEGG) [19] pathways using the mapping files from the KEGG website (downloaded 20 June 2007). Over-representation of each term was tested using the binomial distribution (an efficiently computable and close approximation to the hypergeometric distribution when the number of items to choose among is large [20]). P-values were adjusted using the false discovery rate (fdr) multiple comparison correction [21] and setting the threshold to 0.05.

The Facial Transcriptome

A comprehensive, high-quality gene expression dataset was derived for mouse orofacial development by sampling the mandibular (MdP), maxillary (MxP) and frontonasal (FNP) facial prominences at 0.5-day intervals from E10.5-E12.5, after which point the prominences fused and were no longer separate entities. Seven independent biological replicates were used for each sample, which provided very strong statistical power - a false reading should only be obtained at the rate of 0.06%. In the remainder of the results, we first assess the numerical and biological quality of the dataset. We then provide an overview of gene expression profiles for genes involved in key biological processes, such as cell cycle regulation, and transcription. Lastly, we describe the identification of several genes which either: 1) had not been linked with orofacial development previously, but had been studied in the context of other biological systems; or 2) had little or no prior functional information available, but for which the dataset now suggested a role in orofacial development.

Overview of the Expression Data and Data Quality Assessment

Determination of the number of genes expressed in the dataset was based upon the Detection calls made by the Affymetrix analysis software and any “non-Absent” calls were considered in the expressed category. Using this criterion, more than 20,000 probes were expressed during the period of face formation from E10.5–E12.5. Of 45,101 total probes on the Affymetrix Mouse 430 2.0 microarray, 30% had detectable expression in all 105 samples while 47% were undetectable in any sample. In terms of spatial distribution, 80–82% of the expressed sequences were detected in all three regions of the developing face (Figure 2). A more limited number of probes were expressed in only one prominence or were shared between two of the three prominences (Figure 2). In addition to expression defined by Detection calls, a stringent assessment of prominence-specific expression was obtained by applying the limma package [15] of BioConductor [16] to the Signal values at each time point. Of 20,234 probes tested, 1,506 displayed at least one statistically significant difference between prominences for a given time point. The minimum observed was 171 differences between the MdP and MxP at E10.5, with a maximum of 646 differences in sequence intensity between MdP and FNP at E12.5 (Figure S2A). In general, an increasing number of differences occurred between individual prominences at each successive time point consistent with a divergence of function accompanied by greater complexity of gene expression (Figure S2A). The greatest number of statistically significant prominence-specific gene expression differences occurred between the MdP and FNP, and the least between the MdP and MxP, at all time points. With respect to temporal changes in gene expression, 1,328 probe sets displayed at least one significant increase or decrease between adjacent time points within a particular prominence. This ranged from a minimum of 44 significant changes for the FNP between E11.5 and E12.0, to 582 differences for the MxP between E11.0 and E11.5 (Figure S2B). The data potentially illustrate both an earlier onset of significant gene expression changes in the MdP and circadian aspects to the changes that are occurring in the three prominences. A more in-depth characterization of genes differentially expressed over the entire time series is given later in this section.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Gene expression during mouse facial development.

Venn diagrams showing the number of genes expressed (called Present or Marginal) in the three facial prominences at the five developmental time points analyzed. Values represent the number of probes on the Affymetrix Mouse 430 2.0 array detected by the software in all seven replicates, i.e., called Present or Marginal using the GeneChip Operating Software.

https://doi.org/10.1371/journal.pone.0008066.g002

Reproducibility of the qualitative Detection calls, namely Absent, Marginal and Present, were assessed using the kappa statistic which measures agreement among categorical variables [22] with excellent reproducibility denoted by kappa values >0.75. The median kappa value among replicates was 0.7938 (min = 0.65, max = 0.84). Relatively high kappa values also occurred among non-replicates (median = 0.7649, min = 0.58, max = 0.83) indicative of considerable similarity in global gene expression across the three prominences at all time points. This latter finding suggested that there were unlikely to be widespread radical differences in gene expression between the prominences over the time course of the analyses. Instead, the data was more consistent with subtle alterations in the expression of multiple genes potentially coupled with extensive differences in a limited set of transcripts.

Reproducibility of the quantitative Signal values was assessed using the coefficient of variation (CV), defined as the ratio of the standard deviation to the mean [22]. The CV is zero for perfect reproducibility while higher values indicate poorer reproducibility. If the CV is close to 1, the variability of measurement is on the same scale as the signal and no meaningful analysis is possible [23]. In our dataset the median CV was 0.30 among replicates (min = 0.01, max = 2.59), calculated over all probes in all prominence- and time-specific replicate groups (Figure 3). Much of the high variability was attributable to probes called Absent in all 105 samples (Figure 3, inset), an important distinction since data for all-Absent probes were removed in any further analyses. Disregarding those probes called Absent in all samples, the median CV dropped from 0.30 to 0.24, versus CV = 0.59 among all-Absent probes. The dramatic effect seen when distinguishing among probes based on Absent calls highlights not only the benefits of incorporating the Detection information, but also the importance of utilizing the kappa statistic to assess the quality of that Detection information. High quality of replication in terms of the mean and standard deviation among replicates for each sample was also observed (Figure S3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Coefficient of variation (CV) analysis demonstrates high quality of the dataset.

The main graph illustrates the CV distribution for all probes tested. In the inset graph, the data are represented by the brown line. Removal of the “probes all absent” calls (black line), which show greater signal variation, indicates that the “Present” and “Marginal” probes have even greater reproducibility, illustrated by the shift in the CV curve.

https://doi.org/10.1371/journal.pone.0008066.g003

Biological Verification of the Dataset

For initial biological verification we focused on genes with established prominence-specific expression patterns such as Hand2 (GeneID: 15111), Lhx8 (GeneID: 16875), Elavl4 (GeneID: 15572) and members of the Dlx family [24], [25], [26], [27]. There was excellent agreement between our microarray data (Figure 4A and B) and the expression patterns of these genes as previously determined by whole mount in situ hybridization analysis (WMISH). Thus, the microarray data demonstrated that Dlx1 (GeneID: 13390), Dlx2 and Lhx8 (also known as Lhx7) were preferentially expressed in the MdP and MxP (Figure 4A and B) in agreement with published WMISH results [25], [26]. Previous data had also documented that Dlx3 (GeneID: 13393), Dlx4 (GeneID: 13394), Dlx5 (GeneID: 13395) and Dlx6 (GeneID: 13396) were expressed most highly in the MdP, and this was again reflected at the earlier time points in our microarray data (Figure 4A and B). Of note, signal levels were also detectible for these latter four Dlx genes in the FNP and MxP, a finding that was consistent with in situ hybridization data ([26] and data not shown). A similar concordance was found between our data and the WMISH expression patterns for Hand2 (also known as dHand) a gene that was expressed strongly in the MdP, but not in the MxP and FNP (Figure 4A and B). We chose Elavl4 (Hud) as a marker that should be enriched in the FNP as this gene is associated with neurogenesis and is expressed in the developing olfactory epithelium of the nasal pit [27]. The microarray data again reflected appropriate FNP-specific expression for Elavl4 (Figure 4A and B). Elavl4 is also strongly expressed in the developing peripheral nervous system, including the trigeminal ganglia, but we did not detect significant levels of Elavl4 expression in our MdP and MxP samples until E12. This finding provided additional evidence that our prominence samples were not heavily contaminated with surrounding craniofacial tissue. Thus, our dissections produced prominence specific material, and there was an excellent correspondence between our microarray findings and established differential gene expression patterns.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Biological verification of the dataset.

(A). Raw data for the prominence-specific genes and associated probe sets (left) indicated at the five time points in the three prominences. (B). Corresponding heatmap for the prominence-specific genes and associated probe sets. The heatmaps show the expression data for each probe scaled such that the vector of log₂ expression values for a probe (averaged among replicate samples per time point) has a mean of zero and a magnitude of one. Red and blue indicate high and low expression, respectively.

https://doi.org/10.1371/journal.pone.0008066.g004

Functional Category Analysis of the Temporal Expression Differences

Analysis of over-represented terms among probes with temporal expression differences supported a model in which processes such as cell metabolism and cell growth were declining while differentiation was increasing. Alterations in expression between adjacent time points across the entire time series were characterized using the trajectory clustering algorithm [14]. This approach generates a set of direction change labels between adjacent time points, namely increased, flat, or decreased expression, that can then be extended into a more informative matrix by including all time points. Data collection at five time points generated a total of 81 possible trajectories for each prominence (Table S2). The majority of probes, 15,510 out of 20,234 tested, exhibited consistently flat profiles (no changes) across all five time points in the three prominences. The remaining set of 4,724 probes showing differential expression over time in at least one prominence was much larger than the set of 1,323 probes found using limma alone on adjacent time points, since trajectory clustering captures changes that unfold over longer periods of time.

Our initial analysis concentrated on trajectories with predominantly upward or downward trends. Examination of the trajectories showing at least two decreases and no increases over the time course revealed a broad and statistically significant down-regulation of genes involved in cell division and general metabolism in all three prominences between E10.5–E12.5 (Figure 5 and Tables S3A-J). Seven functional categories associated with cell motility, amino acid metabolism, and other metabolic processes showed a consistent 2-fold reduction in expression during the time course while genes involved with cell cycle, splicing, and translation showed analogous declines in the expression, although less than 2-fold (Figure 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Functional analysis of trajectories: selected categories showing at least two decreases and no increases over the time course.

Over-represented terms, reaching significance at fdr  = 0.5 for the binomial distribution, were based on probes achieving limma statistical significance. Those shown in bold and labeled “(2FC)” reach the additional criterion of a 2-fold change in gene expression.

https://doi.org/10.1371/journal.pone.0008066.g005

Across the prominences, ∼28% of the down-regulated genes were involved in protein metabolic processes, 16% in biosynthetic processes, 18% in transcription, 10% in RNA processing, 8% in translation, 11% in the cell cycle, and 3% in nucleocytoplasmic transport. With respect to individual categories, 25–35% of the genes assigned to glucose metabolism, nucleotide metabolism, DNA replication, DNA repair, RNA processing, amino acid activation, and nucleocytoplasmic transport were found to decrease in expression in the three prominences. Certain terms, including DNA replication initiation, did show a wider range of the percentage of affected genes (Table S3E). The variation in the number of down-regulated genes within such categories may indicate subtle differences in developmental programs among the prominences. Similarly, the assignment of particular terms to individual prominences - or to two of the three prominences - highlights further pathways that may contribute to differential development (Figure 5). In general, though, the relatively small magnitude of the observed gene expression decreases, and the broad range of metabolic functions involved likely reflect an equivalent slowing of growth across all three prominences rather than a drastic change in basic cellular processes.

Conversely, for genes showing an upward trend in expression for at least two of the four intervals and no decreases, the over-represented terms reflected differentiation (Figure 6 and Tables S3K-U). These increases were more robust, and consistently exceeded the criterion of a 2–fold change in the level of expression. Prominent categories displaying expression increases in all three prominences included cell adhesion, cell migration, and cell morphogenesis. Categories for transcription factors, cell signaling, protein interaction domains, calcium-dependent processes, and the extracellular matrix were also highlighted in the various prominences. Notably, ∼50% of the genes associated with phosphate transport showed increasing expression in all prominences (Table S3Q), many by more than 2-fold (Table S3K).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Functional analysis of trajectories: selected categories showing at least two increases and no decreases over the time course.

Over-represented terms, reaching significance at fdr  = 0.5 for the binomial distribution, were based on probes achieving limma statistical significance. All reach the additional criterion of a 2-fold change in gene expression.

https://doi.org/10.1371/journal.pone.0008066.g006

Genes showing an upward trend in expression reflected tissue specialization into teeth, cartilage, and nervous system components with prominence-specific categories more evident than in the gene sets showing decreased expression (Figure 6 compared with 5). In particular, the MdP exhibited increased expression of genes associated with various muscle differentiation categories, reflecting tongue development. The MdP was also over-represented for terms associated with tube morphogenesis and epithelial cell differentiation, presumably reflecting development of the taste buds, teeth, and salivary glands (Figure 6). In contrast, genes associated with the peripheral nervous system and with cyclic nucleotide dependent protein kinase signaling were apparent in the MxP. The FNP was enriched for genes associated with sulfotransferase activity and skeletal development. The MdP, MxP, and FNP were also distinguished by the different protein motifs and protein:protein interaction domains that the trajectory clustering analysis highlighted in the three prominences (Figure 6). We present a further analysis of prominence specific differences in the next section.

In addition to the generally upward or downward gene expression profiles, the trajectory clustering analysis also highlighted several other terms that were over-represented in specific trajectories within a given prominence (Table S2). Two broad themes can be identified from these additional profiles. First, categories in GO Molecular Function, and Biological Process ontologies associated with oxidative phosphorylation showed fluctuating trajectories (up and down) in all three prominences. This finding suggests that energy generation during this period of embryogenesis displays periodicity, consistent with a circadian rhythm that may reflect the nocturnal activity of the mother. Second, ontological terms containing the words “histones” or “chromatin” also had fluctuating trajectories, potentially indicating dynamic changes in transcription and DNA packaging.

Functional Category Analysis of Prominence-Specific Expression

While trajectory clustering analysis highlights dynamic changes in gene expression, this method is not well suited to visualizing instances in which expression may be higher in one of the prominences, but is not increasing or decreasing over the time course. Therefore, to obtain a prominence-specific view of expression changes we studied the set of 1,506 probe sets showing at least one significant (2-fold) difference among prominences at a given time point. This set was further divided into MdP (705 probes), MxP (495 probes), and FNP (760 probes) groups based upon differences in expression from at least one other prominence for the same time point.

These probe set groups were then analyzed with respect to GO, InterPro, MGI Phenotype, and KEGG terms to determine over-represented categories (Figure 7 and Tables S4A-E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Selected categories showing increased expression in the facial prominences.

Prominence-specific groups were created among probes achieving limma statistical significance, with the additional criterion of a 2-fold change in gene expression, by assigning a probe to a group if expression was increased in the group over at least one other prominence at a specific time point.

https://doi.org/10.1371/journal.pone.0008066.g007

Consistent with the data presented in Figure 2, in which we observed significant overlap of detectable expression among the three prominences, the functional characterization again revealed that the prominences had very similar gene expression dynamics since many of the categories observed were seen in the MdP, MxP and FNP. The data also supported our conclusions from the trajectory clustering analysis that significant changes were occurring in growth and differentiation. Alterations in RNA expression levels of transcription factors, signaling molecules, cell adhesion components, and structural proteins were apparent (Figure 7). Cell fate determination, morphogenesis, motility, migration, and adhesion were also prominent categories with ∼30% of the up-regulated genes assigned to cell differentiation, 30% to cell communication, and 10% to cell motility in each of the prominences. There were also notable prominence-specific differences in the distribution of expressed genes within particular categories that may help distinguish the ultimate fate of these separate facial regions. In this regard, 36% of all genes increasing by >2-fold within the MxP were associated with transcriptional regulation, whereas this figure was 29% and 25% for the MdP and FNP, respectively. Similarly, 16% of genes in the FNP were categorized as cell adhesion, 13% in the MxP, and only 10% in the MdP.

Progressive specialization of tissues in the prominences was highlighted by the distribution of over-represented categories associated with development of muscle, bone, cartilage, dentition, vasculature, glands, epithelia, and the nervous system between the MdP, MxP, and FNP (Figure 7). Among probes in the MdP-specific group, the over-represented terms were associated with signal transduction via MAPK and JNK second messenger pathways, as well as with Hedgehog, VEGF, BMP, chemokine, netrin, and neurotransmitter signaling. The MxP group was enriched for terms associated with development of teeth, nerves, and the endocrine system, correlating with formation of the dentition, the trigeminal ganglion, and part of the pituitary system from Rathke's pouch in the oral epithelium. Also prominent in the MxP were genes involved in Wnt signaling and transcriptional control, particularly the negative regulation of transcription. The FNP group was over-represented for categories associated with fibroblast growth factor receptor signaling, sulfate assimilation, ion transport, cell adhesion, cartilage formation, and nerve development, some of which likely reflect development of the specialized olfactory epithelium.

Over-represented categories shared between two of the three prominence groups included muscle, neural crest cell, and peripheral nervous system development in the MdP and MxP. With respect to muscle development, 6.5% and 9.5% of genes were assigned to this category in the MdP and MxP, respectively, but only 2.5% in the FNP group. TGFβ signaling, chemotaxis, angiogenesis, the extracellular matrix, and skeletal formation were highlighted in the MdP and FNP groups. In these two prominences, 3% of genes were in the ECM organization and biogenesis category, whereas the equivalent MxP dataset contained less than 1% of such genes. The MxP and FNP groups shared the categories of transcription factors, retinoic acid, cell adhesion, and claudins. Overall, these findings reflect the widespread importance of signal transduction, transcriptional regulation, and cell:extracellular matrix contacts in growth and morphogenesis of the developing face. With respect to pattern formation, it was noticeable that gene categories associated with different axes of polarity were highlighted in the three prominences. The apparent segregation of radial (MdP), anterior/posterior (MxP), and dorsal/ventral (FNP) patterning in these prominences could represent one important developmental mechanism whereby their individual identities are determined. A further point of interest is that maternally imprinted genes were clearly over-represented in the MdP compared with the MxP and FNP (Figure 7 and Table S4E). Finally, as expected, many genes we identified as expressed in our dataset were linked with the MGI Phenotype categories “cleft palate”, “midline facial cleft” “abnormal cartilage morphology” and “abnormal myogenesis” (Table S4E).

Detailed Examination of Facial Prominence Expression by Functional Category

The assignment of prominence expression changes to particular GO, InterPro, KEGG or MGI Phenotype categories identified some tissue-specific differences, but these were relatively limited compared with the commonalities in cell metabolism, cell growth, and differentiation. Therefore it was necessary to study gene expression changes within a particular category, such as transcription and cell adhesion, in greater detail to identify critical prominence-specific differences. Accordingly, we examined individual genes involved in cell cycle progression, general metabolism, the cytoskeleton and cell adhesion, the extracellular matrix, transcription, and chromatin dynamics. In each category, genes were divided into those that were coordinately regulated in all three prominences (level, increasing, or decreasing) and those that were differentially regulated. With the exception of genes in the “level” group that were coordinately regulated, all the examples presented were identified as differentially expressed by our statistical analyses. Moreover, we concentrated on genes that displayed at least a 2-fold expression difference between the first and last time point. The exceptions, genes and probe sets that are classified as increasing, decreasing, or prominence-specific, but which failed to fulfill this additional 2-fold expression difference criterion, are marked in Figures 8–19 with an asterisk.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 8. Expression profiles of cell cycle regulatory genes.

Each panel shows the expression profile for the gene or genes indicated. The x-axis represents the five time points for each prominence (MdP, MxP, and FNP). The y-axis shows the relative expression level on a log₂ scale – i.e. every integer represents a doubling in the level of expression from the preceding number. The dashed black line indicates the average expression level for all probes. The lines of other colors illustrate the data obtained for the indicated gene shown in the key at the right. If only one gene is shown in a panel the various colors indicate the data obtained from different probes corresponding to that gene. Genes marked by an asterisk in the increasing, decreasing, or prominence-specific categories have statistically significant differences in expression over the time course or between prominences but fail to satisfy the additional criterion of a 2-fold change in gene expression.

https://doi.org/10.1371/journal.pone.0008066.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 9. Expression profiles of genes involved in general cellular metabolism.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g009

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 10. Expression profiles of keratin genes.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g010

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 11. Expression profiles of collagen genes.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g011

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 12. Expression profiles of genes encoding collagen processing enzymes and chaperones.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g012

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 13. Expression profiles of structural genes of the extracellular matrix.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g013

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 14. Expression profiles of genes encoding extracellular matrix modifying components.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g014

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 15. Expression profiles of genes encoding junctional complex proteins.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g015

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 16. Expression profiles of genes encoding cell adhesion molecules.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g016

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 17. Expression profiles of genes encoding transcription factors with coordinately regulated (A–C) or complex (D) patterns.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g017

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 18. Expression profiles of genes encoding transcription factors with coordinate regulation in two of the three prominences.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g018

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 19. Expression profiles of genes encoding chromatin binding and remodeling factors.

Legend as in Figure 8.

https://doi.org/10.1371/journal.pone.0008066.g019

Genes involved in Cell Cycle Progression and General Metabolism

Analysis of genes involved in replication and cell division is shown in Figure 8A–D. In agreement with the Gene Ontology analysis there was an ∼2-fold reduction in the transcript levels of multiple genes involved in promoting cell cycle progression in all prominences (Figure 8A). Such genes included many cyclins and cyclin dependent kinases as well as genes required for nucleotide metabolism, DNA replication, DNA repair, and cell division. Concurrently, transcript levels for genes involved in the negative regulation of cell cycle progression showed a tendency to increase (Figure 8C). Transcripts encoding Rb1 (GeneID: 19645), the G2 checkpoint kinase Plk2 (GeneID: 20620), and several cyclin dependent kinase inhibitors (cdki/cdkn), rose 2 to 4-fold in the prominences. A general increase in expression was also noted for the DNA damage response gene Gadd45b (GeneID: 17873), for cdc42ep3 (GeneID: 260409), a negative regulator of cdc42, (GeneID: 12540) and for Cdkn1c (GeneID: 12577), which encodes p57 cdki. There were some modest prominence-specific differences in Cdkn1c expression, though, with the FNP displaying ∼2-fold lower transcript levels at E10.5 and the MdP having 3 to 4-fold higher levels by E12.5. Nevertheless, the overall similarity of the expression profiles for multiple cell cycle associated genes in the three prominences is striking and suggests coordinate regulation. We note that there was a small subset of such genes that displayed prominence-specific gene expression (Figure 8D). In the FNP, transcripts for DDit4l (GeneID: 73284), Pla2g16 (GeneID: 225845), Rprm (GeneID: 67874), Gadd45a (GeneID: 13197), and Gadd45g (GeneID: 23882) were found at 2 to 4-fold higher levels than in the other prominences. Moreover, expression of Pmaip1 (GeneID: 58801), a gene involved in the induction of apoptosis that is also termed Noxa [28], was 4 to 8-fold lower in the MdP and remained relatively constant throughout the analysis, whereas in the MxP and FNP the Pmaip1 expression profiles were more dynamic, peaking and then decreasing.

Concurrent with decreased expression of genes involved with replication, we also observed ∼2-fold reductions in the expression of multiple genes associated with basic cellular processes, including nuclear import/export, transcription, RNA processing, and RNA degradation, in all three prominences between E10.5 and E12.5 (Figure 9A). Similar transcript reductions were observed for several genes required for ribosome biogenesis, protein synthesis, protein folding, protein degradation, and aspects of mitochondrial function (Figure 9A). Transcripts for many other genes involved in these basic processes remain essentially level during this developmental period (Figure 9B). Only a few genes displayed a slight increase in transcript levels, including several involved in ubiquitination, as well as the catabolism of nucleotides, fatty acids and cholesterol (Figure 9C). There was also an MdP-specific increase in transcript levels for enzymes associated with skeletal muscle function, including Car3 (GeneID: 12350) and Gatm (GeneID: 67092) (Figure 9D), reflecting the development of the tongue and masticatory muscles from this prominence [29]. Taken together, though, these observations indicate that there is a simultaneous overall reduction in the rate of growth and in many aspects of metabolism in all three facial prominences from E10.5–E12.5.

Keratin Genes

We next examined whether expression of genes involved in differentiation were activated to counter-balance the decreased expression of those associated with replication. Several classes of genes were chosen to study the process of differentiation, including intermediate filament proteins (Figure 10), extracellular matrix molecules (Figures 11–14), cell adhesion molecules (Figures 15, 16), and transcription factors (Figures 17, 18). As described in detail below, these studies demonstrated that the three sets of prominences exhibited striking commonalities in the onset and progression of differentiation. We focused initially on keratin genes (Figure 10A, B) because their expression provides an excellent readout of development of the surface ectoderm from a simple to multilayered epithelium [30]. Over the time course, there was a 3 to 4-fold reduction in Krt8 (GeneID: 16691) and Krt18 (GeneID: 16668) transcripts, which are associated with simple embryonic epithelia, and the levels of Krt7 (GeneID: 110310) and Krt19 (GeneID: 16669) transcripts also dropped initially before recovering (Figure 10A). Concurrently, there was a rise in the transcript levels for keratins associated with more mature stratified epithelia (Figure 10B). Transcripts derived from Krt5 (GeneID: 110308), Krt14 (GeneID: 16664), and Krt15 (GeneID: 16665), which are expressed in the basal epithelial layers, rose exponentially from E10.5 onwards. Krt1 (GeneID: 16678) and Krt10 (GeneID: 16661) transcripts, associated with the suprabasal epidermal layers, were increasing more slowly during the same period. Lastly, Krt6a (GeneID: 16687) and Krt17 (GeneID: 16667), which are expressed in the suprabasal layers of epithelia of the oral cavity, including the tongue [30], began to rise at the later time points.

The switch from keratin genes expressed in simple epithelia to those present in stratified epithelia occurred simultaneously in the three prominences. Nevertheless, evident differences in the three prominences were superimposed on this overall pattern (Figure 10A, B). In the FNP the switch between keratins associated with simple (Krt8 and Krt18) versus complex (Krt15) epithelia progressed at a much slower rate. The MxP also displayed differences from the other prominences, especially in the levels of Krt5 and Krt14 expression, which were higher, and the expression profiles of Krt6a and Krt7, which were less dynamic. MdP specific differences were less pronounced. Overall, though, each prominence exhibited a unique molecular fingerprint of keratin gene expression. These prominence-specific differences were also reflected in the expression of Trp63 (GeneID: 22061), encoding the p63 transcription factor that promotes the commitment to epithelial stratification [31]. Transcripts encoding p63 were lowest in the FNP reflecting the slower rate of change between keratins characteristic of simple versus stratified epithelia in this prominence (Figure 10C).

Genes Encoding Components of the Extracellular Matrix

Genes Encoding Cell:Cell Communication and Cell Adhesion Molecules

The Gene Ontology and InterPro analyses shown in Figures 6 and 7 and Tables S3 and S4 highlighted gene categories associated with cell:cell adhesion as increasing in all prominences, and also suggested that claudin family members were differentially expressed in the FNP. Examination of individual genes associated with gap junctions, tight junctions, and desmosomes confirmed these dynamic and tissue-specific expression changes (Figure 15A–D). Genes encoding catenin proteins, which link the cytoskeleton with adherens junctions [34], showed coordinate decreasing (Ctnnal1, GeneID: 54366; Fig 15A), level (Ctnna1, GeneID: 12385; Ctnnb1, GeneID: 12387; and Ctnnd1, GeneID: 12388: Fig 15B), or increasing (Ctnna2, GeneID: 12386; Fig. 15C) expression profiles in the three prominences. Similarly, expression of the gap junction protein gene Gjb3 (GeneID: 14620, connexin 31) was decreasing in all prominences, while there was a concomitant rise in Gjb2 (GeneID: 14619, connexin 26) and Gjb6 (GeneID: 14623, connexin 30) expression. Genes encoding desmosomal proteins, including Dsc2 (GeneID: 13506), Dsp (GeneID: 109620), and Perp (GeneID: 64058), were also increasing in expression within in the MdP, MxP, and FNP between E10.5 and E12.5. These findings reflect increased adhesion, cell communication, and differentiation within epithelial sheets shared by the three facial prominences. Several genes associated with tight junctions, including Tjp1 (GeneID: 21872), Tjp2 (GeneID: 21873), and Tjp3 (GeneID: 27375), had slightly decreasing or essentially level expression in the facial prominences. In contrast, the expression profiles for many other genes encoding adherens and tight junction proteins [35], [36] displayed markedly different dynamics in the FNP from the other two prominences (Figure 15D). With respect to adherens junctions, Pkp2 (GeneID: 67451), Ssx2ip (GeneID: 99167) and Shroom3 (GeneID: 27428) all showed higher expression in the FNP, as did the tight junction genes Cgn (GeneID: 70737) and Ocln (GeneID: 18260). Higher expression of the tight junction gene Cldn1 (GeneID: 12737) was characteristic of the MdP and MxP, whereas Cldn3 (GeneID: 12739), Cldn6 (GeneID: 54419), Cldn7 (GeneID: 53624), Cldn9 (GeneID: 56863) were expressed at greater levels in the FNP throughout this developmental window. Cldn11 (GeneID: 18417) and Cldn23 (GeneID: 71908) expression also rose more rapidly in the FNP than in the other two prominences. Overall, these findings support a different composition and/or function of the junctional complexes in the FNP compared with the MdP and MxP. In this context, the regulatory molecule protein kinase C zeta (Prkcz, GeneID: 18762), which associates with tight junctions [35], [36], is also expressed at a higher level in the FNP. The MdP had higher levels of Amot (GeneID: 27494) expression, presumably associated with accelerated vasculature development in this prominence, as well as the cell:matrix adherens junction gene Nexn (GeneID: 68810).

The dynamic changes in the expression of genes encoding cell:cell junction components prompted a more detailed examination of genes involved in intercellular communication and cell adhesion (Figure 16A–D). A subset of these genes was coordinately down-regulated as development progressed (Figure 16A), especially Igdcc3 (previously known as Punc, GeneID: 19289) whose expression decreased by 15 to 30-fold during our analysis. In agreement with our findings, previous reports showed that Igdcc3, which encodes a neural cell adhesion molecule, has widespread expression in the mesenchyme of the facial prominences at E10.5, but soon after becomes limited to discrete regions of the CNS [37], [38]. Several other genes had relatively stable expression levels (Figure 16B), including Pvrl1 (GeneID: 58235), a gene that when mutated results in human CL/P [39], [40], [41]. However, consistent with the onset of differentiation, a majority of expression profiles were coordinately up-regulated in the facial prominences, including genes that encode specific cell adhesion molecules, cadherins, and protocadherins, which exhibited 2 to 10-fold increases (Figure 16C).

Prominence specific expression profiles were observed for several additional genes involved in cell adhesion (Figure 16D). In particular, in the MdP expression of Cdh15 (GeneID: 12555), a gene associated with muscle differentiation [42], showed a sharp increase compared to the other prominences. Relatively low, stable, levels of Pcdh8 (GeneID: 18530) and Pcdh19 (GeneID: 279653) transcripts were also characteristic of the MdP. Distinctive expression profiles were not observed for the MxP, a finding in common with the data obtained from extracellular matrix gene analysis (see above). However, transcript levels for Itga8 (GeneID: 241226) and Itgb5 (GeneID: 16419) were slightly higher in the MxP compared to the other prominences. Several genes showed higher basal expression levels in the FNP, including those encoding the cell adhesion molecules Pcdh20 (GeneID: 219257), Pcdh21 (GeneID: 170677), Cadm1 (GeneID: 54725), and Epcam (GeneID: 17075). In this context, Pcdh21 is associated with olfactory function [43], and it is likely that some of the differences between the FNP and the two other prominences reflect the specialized function of the olfactory epithelium.

Transcription Factors

Although we observed significant changes in the expression of genes involved in cell cycle, cell adhesion, cell structure, and the extracellular matrix, by far the most dramatic molecular signatures of the prominences were provided by genes encoding sequence-specific DNA binding transcription factors (Figures 17A–D, 18A–C). Multiple transcription factors displayed coordinated gene expression in the three prominences over the time course of this analysis (Figure 17A–C and Figure S6). A small number of genes were down-regulated (Figures 8A and 17A) including those encoding p53 (GeneID: 22059), E2F5 (GeneID: 13559), and Sp6 (GeneID: 83395) whose expression dropped 2 to 4-fold correlating with the reduction in cell proliferation and general cellular metabolism noted above. Transcription factors implicated in maintenance of stem cell function, including Sall4 (GeneID: 99377) and Zbtb16 (Plzf, GeneID: 235320) [44], [45], [46], or in the expansion of precursor populations, such as Arid3b (GeneID: 56380) [47], were also down-regulated ∼5-fold. Several genes were expressed at relatively stable levels in the three prominences, including members of the Stat and Smad transcription factor families and Cited2 (GeneID: 17684), which encodes a transcriptional modulator (Figure 17B). The expression of a much larger group of genes was up-regulated in all prominences (Figure 17C). This class of genes includes Fosl2 (GeneID: 14284), Runx2 (GeneID: 12393), and Sp7 (osterix, GeneID: 170574), which are required for development of skeletal elements [48], [49], [50], and whose transcript levels rose 5 to 15-fold during this developmental window. Similarly, expression of Klf4 (GeneID: 16600), a gene responsible for skin barrier function in the differentiating epidermis [51], was increased ∼4 to 8-fold in all three prominences. Irf6 (GeneID: 54139), a gene expressed in facial ectoderm and critical for patterning mouse orofacial development as well as mutated in human CL/P [52], [53], [54], was also present in this category. Several genes that encode components of signal transduction pathways in many cell types, including Stat3 (GeneID: 20848), and the nuclear hormone receptors, corepressors, and coactivators Ncoa1 (GeneID: 17977), Ncor2 (GeneID: 20602), Nr3c1 (GeneID: 14815) and Thra (GeneID: 21833), were also coordinately up-regulated. In addition, all members of the NFI transcription factor gene family demonstrated coordinate 10 to 50-fold increases in expression over this time period in the facial prominences. A further group of transcription factor genes displayed distinctive patterns of expression in each of the three prominences that did not fall into the simple categories noted above (Figure 17D). This group was quite small and included genes such as Sox8 (GeneID: 20681) and Irx4 (GeneID: 50916), which had transcript levels rising and falling in different prominences. Another example was Foxg1 (GeneID: 15228), with stable expression in each prominence, but with a low, intermediate, and high level of expression in the MdP, MxP and FNP respectively.

The microarray data were also processed to ascertain those transcription factor genes that exhibited prominence-specific expression or were coordinately regulated in any two of the three prominences (Figure 18, and Figures S7 and S8). This analysis revealed a complex series of expression profiles for each prominence, with the MdP and FNP showing the greatest differences. With respect to the MdP, our analysis showed that several genes required for muscle differentiation [29] including Tcf21 (capsulin, GeneID: 21412), Myf5 (GeneID: 17877), Myod1 (GeneID: 17927), Myog (GeneID: 17928) and Msc (GeneID: 17681), were expressed at higher levels in this prominence (Figure 18A). Expression of Tcf21 was essentially level throughout the time course, with levels 2 to 3-fold higher than in the MxP and FNP. Expression of the other genes increased between 4 to 30-fold between E10.5–E12.5, presumably reflecting development of the tongue and the muscles of mastication [55]. Additional genes more highly expressed in the MdP included Dlx5, Dlx6, Gsc, Hand1 (GeneID: 15110), and Hand2, consistent with their developmental roles in patterning the mandible [24], [26].

The MxP shared multiple expression profiles with either the MdP or FNP, and only a few genes showed expression differences specific to this prominence – and most of these were subtle (Figure 18B). Thus, Foxc2 (GeneID: 14234) and Pknox2 (GeneID: 208076) showed a slight decrease in expression in the MxP compared with the other prominences, whereas levels of Foxl2 (GeneID: 26927), Irx3 (GeneID: 16373) and Irx5 (GeneID: 54352) were slightly higher. The FNP displayed higher expression for genes encoding transcription factors involved in specification and development of ectodermal placodes and in neurogenesis [56], [57], including Pax6 (GeneID: 18508), Sox2 (GeneID: 20674), Ascl1 (GeneID: 17172), Olig1 (GeneID: 50914) and Olig2 (GeneID: 50913) (Figure 18C). This observation again presumably reflects the development of the olfactory placodes and specialized olfactory epithelium in the FNP. Transcripts of genes involved in patterning the FNP, such as Pax7 (GeneID: 18509), were also expressed at a higher level in this prominence (Figures 17D). The integration of the multiple differences in transcription factor gene expression profiles is likely to provide a powerful driving force for discriminating the different prominences during development.

Genes Associated with Chromatin Dynamics

We also examined genes associated with chromatin structure (Figure 19A–D) since alterations in DNA packaging would be expected to accompany the rapidly changing profiles of transcription factors, as well as the switch between genes required for cell cycle progression versus differentiation. As indicated by the trajectory clustering analysis (Table S2), we observed dynamic changes in expression of genes associated with chromatin remodeling, DNA methylation, and gene regulation. Our findings imply that there is a switch from a progenitor to a differentiated cell program of gene accessibility during the period under analysis in all prominences. Thus, expression of genes encoding particular members of the Chd family of chromatin remodeling proteins was decreasing (Chd1, GeneID: 12648; Chd7 GeneID: 320790; Figure 19A), while others remained relatively constant (Chd8, GeneID: 67772; Figure 19B), and at least one was increasing (Chd9, GeneID: 109151; Figure 19C). These findings may reflect the difference between maintaining chromatin assembly functions and epigenetic marks during DNA replication versus the established role of Chd9 in regulating genes required for skeletal tissue development [58]. An alteration in chromatin dynamics was also indicated by the increasing expression of genes encoding many components of the SWI/SNF remodeling complex, including Smarca1 (GeneID: 93761), Smarca2 (GeneID: 67155), Smarcc2 (GeneID: 68094) and Smarcd3 (GeneID: 66993) (Figure 19C), whereas others such as Smarcc1 (GeneID: 20588) and Smarce1 (GeneID: 57376) showed relatively constant expression levels (Figure 19B). The majority of these chromatin remodeling genes showed coordinate regulation of expression, with the exception of Baz1a (GeneID: 217578), which was preferentially down-regulated in the MdP (Figure 19D). There was also a general increase in transcripts for genes involved in histone acetylation, including Kat2b (Pcaf, GeneID: 18519) and Myst4 (GeneID: 54169), suggesting an overall increase in chromatin accessibility. Simultaneously, expression of genes encoding specific components of the histone deacetylase and methyltransferase machinery declined over this developmental window, including Sap18 (GeneID: 20220), Sap30 (GeneID: 60406), and Suv39h2 (GeneID: 64707). Interestingly, though, transcripts for the HDAC components Hdac5 (GeneID: 15184) and Chd3 (GeneID: 216848) were rising. There was also a 2-fold increase in transcripts for the deacetylase Phf21a (GeneID: 192285), which may reflect the importance of this gene in the repression of neural genes in non-neural tissues [59].

With respect to polycomb group and trithorax genes, expression of Asxl3 (GeneID: 211961), Cbx4 (GeneID: 12418), Cbx6 (GeneID: 494448) and Mll5 (GeneID: 69188) increased during orofacial development (Figure 19C). These genes encode the archetypal epigenetic marker proteins responsible for repression of inappropriate genes in differentiated cells [60]. In contrast, transcripts corresponding to Pcgf6 (GeneID: 71041), Cbx1 (GeneID: 12412), and Cbx2 (GeneID: 12416) decreased between E10.5–E12.5 in the facial prominences. These data imply that such polycomb proteins could be responsible for maintaining progenitor cells by repression of genes required for commitment to differentiation. Dnmt3b (GeneID: 13436), encoding an enzyme involved in de novo CpG DNA methylation, which is associated with transcriptional repression and heterochromatin formation [61], was also down-regulated between E10.5–E12.5 (Figure 19A). In contrast, Dnmt3a (GeneID: 13435) showed flat or increasing expression depending on the probe used (data not shown). The timing of the switch between these related DNA methylases, also observed in embryos at the protein level [62], potentially indicates that they have different roles in early progenitor versus differentiated cells. The dynamics of Dnmt3b expression are also of interest given its mutation in human immunodeficiency-centromeric instability-facial anomalies (IFC) syndrome (OMIM: 242860) [61]. Taken together, the changes in the expression of genes involved in DNA modification and chromatin structure demonstrate that the genome is undergoing rapid and dynamic changes that will greatly alter transcription factor accessibility.

Analysis of Novel Facial Gene Expression

Our studies highlighted a large number of genes not previously associated with the developing face. This group of genes included well-studied genes as well as novel genes for which little or no functional information exists. We selected a number of these genes for further analysis based, in part, upon their high levels of expression in discrete prominences as determined by the microarray analysis (Figure S9). RNA in situ probes were designed for each gene and WMISH was employed to determine the location of expression within each prominence between E10.5–E12.5. The spatial distribution of transcripts observed using WMISH displayed excellent agreement with the microarray expression results, providing additional verification of our dataset (Figures 20, 21, and Figure S9). Figure 20A–D shows the WMISH data obtained for cDNA 2610017I09Rik (GeneID: 66297), a transcript of unknown function [63], that was readily detected by microarray analysis in each prominence at all time points studied. WMISH confirmed that this gene was expressed highly in all three prominences, but also indicated that the distribution was not uniform. Expression of 261007I09Rik showed a proximal distribution in the MdP and MxP, whereas in the FNP it was more concentrated in the medial rather than lateral nasal prominences. 2610017I09Rik maps to mouse chromosome 1, in close vicinity to Pou3f3 (GeneID: 18993) (Figure S10), with which it shares a similar expression profile in the facial prominences and CNS.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 20. Expression of representative genes in the facial prominences I: global and FNP-specific genes.

Whole mount in situ expression analysis of: 2610017I09Rik (A–D); 9430073N08Rik, Fam162b (E, F); Elav14 [HuD] (G); Prr15 (H, I); Fam107a (J–L); Wscd1 (M–P); and Image 6309403 (Q, R) showing lateral (A, D, E, G, H, J, M, P, Q) or ventral (B, C, F, I, K, L, N, O, R) views of E10.5 (A, B, E–K, M, N, Q, R) or E11.5 (C, D, L, O, P) mice. Note for panel G that despite the clear expression of Elavl4 in the trigeminal ganglia, we did not detect expression in the MxP and MdP RNA samples we had purified for microarray analysis. This finding provides additional proof that the samples we have purified were not significantly contaminated with surrounding tissue. Abbreviations: drg, dorsal root ganglia; e, eye; fn, frontonasal prominence; h, hyoid arch; l, lateral nasal prominence; m, medial nasal prominence; md, mandibular prominence; ov, otic vesicle; p, nasal pit; t, telencephalon; tg, trigeminal ganglion; x, maxillary prominence.

https://doi.org/10.1371/journal.pone.0008066.g020

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 21. Expression of representative genes in the facial prominences II: MdP and MxP-specific.

Whole mount in situ expression analysis of: Ifitm3 (A, B); Lrba (C–E); AI606473 (F, G); A930038C07Rik (H–K); Gpr50 (L–N); Rgs5 (O–Q); and AV026068 (R–T); showing lateral (A, C, E, F, H, K, L, O, T) or ventral (B, D, G, I, J, M, N, P–S) views of E10.5 (A–D, F–I, L, M, O, P, R) or E11.5 (E, J, K, N, Q, S, T) mice. Abbreviations: e, eye; fl, forelimb bud; fn, frontonasal prominence; h, hyoid arch; hl, hindlimb bud; l, lateral nasal prominence; m, medial nasal prominence; md, mandibular prominence; mg, midgut; ov, otic vesicle; p, nasal pit; rp, Rathke's pouch; t, telencephalon; x, maxillary prominence.

https://doi.org/10.1371/journal.pone.0008066.g021

The additional transcripts highlighted in Figure 20 show preferential expression in the FNP over the other prominences, in agreement with the microarray data. One of these, corresponding to the cDNA 9430073N08Rik (GeneID: 77296), was expressed in a portion of both the lateral and medial components of the nasal prominence proximal to the oral cavity (Figure 20E, F). This latter gene, Fam162b, is predicted to encode an integral membrane protein, but little else is know about its function [63]. For the rest of the genes that display FNP-specific expression, we utilized the distribution of Elavl4 (Hud) transcripts as a landmark ([27]; Figure 20G). Expression of Elavl4 was elevated in the FNP and localized to the developing olfactory epithelium of the nasal pit consistent with its association with neurogenesis. The other FNP-enriched transcripts analyzed (Prr15, GeneID: 78004; Fam107a, GeneID: 268709; Wscd1, GeneID: 216881; and Image:6309403, GenBank gi: 21854763) were similarly localized in the nasal pit (Figure 20H–R), presumably reflecting a role for these genes in olfactory neurogenesis. Expression of Prr15 (also termed G90), which may encode a proline-rich protein or act as a non-coding RNA, has previously been observed in post-mitotic cells, including within the olfactory epithelium ([64]; Figure 20H, I). Fam107a (Figure 20J–L) (also known as BC055107, DRR1 and TU3A) may encode a nuclear protein linked to growth control and tumorigenesis [65], [66]. Wscd1 (Figure 20, M–P) is a wsc domain-containing protein predicted to be integral to the cell membrane [63]. The Image:6309403 transcript (Figure 20Q, R) is predicted to encode microRNA 429 (mir429, GeneID: 723865) of unknown function [63].

WMISH was also used to examine genes that showed high expression in the MdP and/or MxP in the microarray analysis, and once again we observed good agreement between the data obtained using the two methods (Figure 21 and Figure S9). Ifitm3 (GeneID: 66141), Lrba (GeneID: 80877), and AI606473 (GeneID: 99686) were expressed in both the MdP and MxP, although with clearly different spatiotemporal distribution patterns. Ifitm3 (interferon induced transmembrane protein 3) [67] was expressed primarily in the ectoderm of the cleft between the MdP and MxP (Figure 21A, B). Lrba, which encodes an LPS-responsive beige-like anchor protein that may be involved in vesicle trafficking and Notch regulation [68], was expressed more diffusely in the MdP and MxP as well as in the more caudal branchial arches (Figure 21C–E). Expression of AI606473, a gene of unknown function [63], occurred primarily in the MxP and the distal portion of the MdP, although a discrete domain of expression was also seen in the medial nasal prominences (Figure 21F, G). This expression pattern showed strong similarity with that of Lhx8 ([25]; Figure S9). AI606473 is located just upstream of Lhx8 on mouse chromosome 3 and is transcribed in the opposite direction (Figure S10), suggesting that these two genes may share common cis-regulatory elements and/or a bidirectional promoter.

Other genes analyzed displayed more prominence-specific gene expression patterns. A930038C07Rik (GeneID: 68169) showed the highest levels of expression in the MxP, particularly between E11–E12, and serves as a valuable new marker for this prominence (Figure 21, H–K). This gene, known as C4orf31 (GeneID: 79625) in human, encodes a hypothetical membrane protein and is negatively regulated by Hox genes in the developing kidney [69]. This latter observation is intriguing given that the anterior-posterior Hox code terminates just caudal to the MxP [70]. Specific expression in the MdP was seen for two genes linked with G protein coupled signaling [71], [72], Gpr50 (GeneID: 14765) (Figure 21, L–N) and Rgs5 (GeneID: 19737) (Figure 21, O–Q). Expression of Gpr50 was concentrated in the rostral MdP, whereas Rgs5 was expressed more distally. At E10.5, both genes were expressed at high levels, but expression decreased considerably from E11.5 onwards. Expression of Gpr50, a melatonin-related orphan receptor, was also associated with Rathke's pouch in the roof of the oral cavity, suggesting a role in the development of the pituitary gland. AV026068 (GeneID: 102265), termed NBLA00301 in human (GeneID: 79804), encodes a short hypothetical protein of unknown function [63], and its transcripts were expressed at elevated levels in the MdP throughout the developmental time course (Figure 21, R–T). The expression domain of AV026068 in the distal MdP, as well as in the posterior region of the limb bud, resembled that of the aforementioned Hand2 gene (Figure 4). AV026068 is positioned adjacent to Hand2 on mouse chromosome 8, but is expressed from the opposite strand beginning several kb upstream of the latter's transcriptional start site (Figure S10). This arrangement is remarkably similar to those described above for 2610017I09Rik and Pou3f3 as well as Lhx8 and AI606473. Again, the similarity between the expression patterns of AV026068 and Hand2, as well as the bidirectional nature of these two gene's transcription units, suggests that they share common mandible-specific cis-regulatory sequences. The finding that this arrangement occurs for these three new pairs of genes indicates that it may be a common mechanism of coordinating gene expression in the developing mouse face.