AP-2α and AP-2β cooperatively orchestrate homeobox gene expression during branchial arch patterning

Development of the facial complex in mammals begins as an outgrowth of structural protrusions, the facial prominences, surrounding the nascent oral cavity. These outgrowths include a separate pair of mandibular (MdP) and maxillary (MxP) prominences that constitute the first branchial arch (BA1) with a single frontonasal prominence (FNP) superior to these structures (see Fig. S1) (Trainor, 2014). The mesenchyme of these prominences is largely composed of neural crest cells (NCCs), which migrate from the dorsal mid/hindbrain and are involved in outgrowth and morphogenesis associated with facial development (Trainor, 2014). Signals from the adjacent ectoderm, endoderm and neural tube are key sources of positional identity for the invading NCCs, e.g. regulating the establishment of dorsal-ventral identity within BA1 (Trainor, 2014). For BA1, this ultimately leads to formation of both maxillary (upper-jaw) and mandibular (lower-jaw) units, with distal cap regions and a hinge domain at the articulation point of the MdP and MxP, typifying the bauplan of the vertebrate jaw (Fish et al., 2011) (see Fig. S1B-E). Notably, modulation of these signals, or the ability of NCCs to respond to them, can ultimately influence final adult structures emerging from this bauplan. Such plasticity is crucial for the array of facial features observed among vertebrate species but also provides an opportunity for disruption in a variety of craniofacial disorders.

Numerous studies have shown that homeobox transcription factors (TFs) are key players in the regulation of BA1 patterning. Specifically, nested expression of the DLX family of homeobox TFs (DLX1-DLX6), is crucial for establishing axes essential for craniofacial patterning within the jaws (Fig. S1F,G) (Cobourne and Sharpe, 2003). Additional homeobox TFs, including Msx-family members (Alappat et al., 2003), Alx-family members (McGonnell et al., 2011), Six1 (Tavares et al., 2017), Emx2 (Compagnucci et al., 2013), Gsc (Rivera-Perez et al., 1995; Yamada et al., 1995) and Barx1 (Barlow et al., 1999; Sperber and Dawid, 2008), have been implicated in associated aspects of BA1 patterning. Thus, a craniofacial gene regulatory network (GRN) is emerging for BA1 composed of numerous homeobox TFs that drive NCC patterning. However, the upstream regulators of homeobox TF expression are less well understood, and are likely key to priming cranial NCCs for the cues they will encounter while invading BA1. Identifying these regulators would provide additional details into autonomous nodes that can be used to influence craniofacial form during normal and abnormal development.

One gene family postulated to be crucial throughout NCC development are the AP-2 TFs, which consist of five paralogs in mammals: Tfap2a to Tfap2e, encoding the proteins AP-2α to AP-2ε, respectively (Eckert et al., 2005; Green et al., 2015; Simoes-Costa and Bronner, 2015). These proteins all recognize a similar consensus sequence (Eckert et al., 2005), which raises the possibility of functional redundancy in regulating target gene expression when multiple family members are expressed simultaneously (Bassett et al., 2012; Hoffman et al., 2007; Knight et al., 2005; Li and Cornell, 2007; Schmidt et al., 2011; Seberg et al., 2017; Van Otterloo et al., 2010; Wang et al., 2008). The potential for coordinated action is particularly pertinent to the early NCCs, where both Tfap2a and Tfap2b are highly expressed (Fig. S2A-E) (Eckert et al., 2005). Currently, Tfap2a is the best characterized with respect to NCC development in multiple model organisms (Brewer et al., 2004, 2002; de Croze et al., 2011; Knight et al., 2003; Luo et al., 2003; O'Brien et al., 2004; Schorle et al., 1996; Shen et al., 1997; Zhang et al., 1996). However, NCC roles for additional paralogs, notably AP-2β (Tfap2b), have also been identified in mouse and non-mammalian species (Knight et al., 2005; Martino et al., 2016; Schmidt et al., 2011; Simoes-Costa and Bronner, 2016). Further evidence suggesting crucial roles for AP-2α and AP-2β in NCCs has been obtained from the analysis of human branchio-oculo-facial (Milunsky et al., 2008) and Char syndromes (Satoda et al., 2000), respectively, that display defects in NCC derivatives in the face, heart and skin. Finally, NCC enhancers mapped in human and chimp contain an over-representation of AP-2-binding sites and such sites correlate with active histone modifications and elevated gene expression (Rada-Iglesias et al., 2012). Together, these findings raise the possibility that the AP-2 genes have been integral to the evolution of the NC and that modulation of AP-2 regulatory circuits has potentially influenced species-specific craniofacial traits during evolution (Prescott et al., 2015). Despite strong evidence that AP-2 is crucial within the hierarchy of NCC development, single mouse NCC-specific gene knockouts for Tfap2a and Tfap2b – the most pertinent members expressed in early NCCs – produce only partially penetrant (∼44%) cleft secondary palate or eye defects, respectively (Brewer et al., 2004; Martino et al., 2016) – with some animals remaining viable into adulthood. To address this discrepancy, we generated mice lacking both Tfap2a and Tfap2b within the NCCs. Our results reveal that these two AP-2 TFs act in concert within NCCs in a cell-autonomous manner to regulate major aspects of craniofacial development. These data resolve the paradox of the mild phenotypes seen with single conditional mutants, demonstrate the importance of AP-2 in the cranial BA1 NCC GRN, and provide an important framework for analysis of vertebrate evolution and human craniofacial birth defects.

Although three of the five mouse Tfap2 genes, Tfap2a, Tfap2b and Tfap2c, are expressed at appreciable levels in the cranial NCC, only mRNAs from Tfap2a (encoding AP-2α) and Tfap2b (encoding AP-2β) are found at high levels at E8.5 and E9.5 when the NC is emerging from the margins of the neural tube and migrating to the facial prominences (Fig. S2A-E). Surprisingly, previous NCC-specific deletion of either AP-2α or AP-2β alone resulted in relatively limited craniofacial defects (Brewer et al., 2004; Martino et al., 2016), and therefore we considered that they might be functionally redundant. To address this issue, we generated a Tfap2a/Tfap2b NCC-specific conditional double mutant using a new Tfap2b-conditional allele (Fig. S3A-C), in conjunction with the Tfap2a-floxed allele (Brewer et al., 2004). Next, we bred mice heterozygous for Tfap2a- (Zhang et al., 1996) and Tfap2b-null (Fig. S3D) alleles that also contained the Wnt1:CRE transgene (Danielian et al., 1998) with mice homozygous for both conditional alleles, to target Tfap2a/Tfap2b in pre-migratory NCCs. The expected offspring from these matings included the desired double conditional mutants, Wnt1:CRE; Tfap2a^null/flox; Tfap2b^null/flox, that we term DCMs for the remainder of this study. In addition, other allelic combinations were generated with different Tfap2a and Tfap2b gene doses (Table S1A). These included mice that have lost both Tfap2a alleles in the NCC but retain one functional allele of Tfap2b (single conditional mutant AP-2α; abbreviated SCM-A), or those missing both Tfap2b alleles but possessing one functional copy of Tfap2a (single conditional mutant AP-2β: abbreviated SCM-B). The availability of the DCM, SCM-A, SCM-B and double-conditional heterozygous (DCH) mice enabled us to compare how various Tfap2a/Tfap2b allelic combinations impacted development. Western blot analysis revealed that, compared with controls, SCM-A and SCM-B embryos showed a modest reduction in AP-2 protein levels, whereas DCMs displayed in the presence of a CRE recombinase transgene, allowing genetic dissection of their redundant function in the NC (Fig. S2F,G).

We next addressed the developmental consequences of NCC-specific deletion of Tfap2a and Tfap2b on gross craniofacial development by examining embryos at E18.5 (Fig. 1A-H). Three phenotypic classes were identified, each corresponding to particular genotypes: (1) embryos with relatively normal craniofacial structures that were either control, DCH or SCM-B (Fig. 1A,B,E,F, and not shown); (2) embryos with a slightly malformed nasal bridge, slight micrognathia and a cleft secondary palate that were SCM-A (Fig. 1C,D); and (3) embryos with a complete midface cleft, shortened and cleft mandible, low set ears with microtia, microglossia and a cleft secondary palate that were DCM (Fig. 1G,H). In SCM-A embryos, the cleft secondary palate was 100% penetrant, in contrast to the limited penetrance (44%) of this phenotype in Tfap2a-Wnt1:CRE mutants (Brewer et al., 2004). In addition, a small proportion (∼22%) of SCM-B embryos developed cleft secondary palate (Table S1A), a feature not yet observed in Tfap2b-Wnt1:CRE mutants (Martino et al., 2016). The increased incidence of clefting in SCM-A and SCM-B strains compared with the corresponding single knockouts presumably reflects the detrimental consequence from loss of one additional Tfap2 allele. Collectively, several conclusions can be derived from the analysis of gross morphology. First, AP-2α and AP-2β cooperatively regulate craniofacial development cell-autonomously in NCCs. Second, AP-2α may play a more prominent role than AP-2β in this process. Third, there is a phenotypic spectrum, dependent on the number of functional Tfap2a/Tfap2b alleles present in the NCC – implying dose sensitivity to AP-2 levels during craniofacial development.

Given the cell-autonomous role for AP-2α and AP-2β within cranial NCCs, and the emphasis placed on AP-2α as a ‘master regulator’ of NCCs (Green et al., 2015; Hong et al., 2014; Rada-Iglesias et al., 2012, 2013; Simoes-Costa and Bronner, 2015; Trainor, 2014), we next examined NCC developmental progression in SCM-A, SCM-B and DCM embryos when compared with controls. To assess initial specification, emigration and migration of NCCs into the craniofacial complex, we conducted the genetic cross described above but also included the Rosa26 reporter allele in the breeding scheme (Soriano, 1999). In conjunction with the Wnt1:CRE transgene (Danielian et al., 1998), this reporter allele allows precise indelible labeling of NCCs, detected by β-galactosidase (β-gal) staining. At E9.0, streams of migratory β-gal⁺ NCCs were clearly identifiable in both the first and second BA of control embryos with no obvious differences noted in either SCM or DCM embryos (Fig. S4A,B and data not shown). However, by E10.0-E10.5, although similar labeling patterns were detected between controls (Fig. 2A, Fig. S4C), SCM-A (Fig. 2B), SCM-B (Fig. 2C) and DCM embryos (Fig. 2D, Fig. S4D), differences in the size of the BAs were becoming apparent in SCM-A and DCM embryos (Fig. 2B′,D′, Fig. S4C,D and not shown). We then examined whether differences in cell proliferation and/or cell death could account for the observed morphological changes in the DCMs using α-phospho-histone H3 (pH3) and α-cleaved caspase 3 (CC3) immunofluorescence, respectively. Grossly, cell proliferation was not dramatically altered in mutant versus control embryos at E10.5 in BA1 (Fig. S4E,F) or within the anterior dorsal neural tube at E8.5 (data not shown). However, analyses of cell death in adjacent sections identified localized groups of cells that were CC3⁺ in DCM embryos. Such cells were largely observed around the intersecting hinge domain of the developing upper and lower jaw in DCM embryos (Fig. 2E,F), consistent with the observed BA1 hypoplasia. These latter findings indicate that AP-2α/β are crucial for preventing the death of a small subset of BA1 NCCs, but do not appear to have a global role in NCC survival.

We next examined specification and differentiation of cranial NCC derivatives, including the peripheral nervous system (PNS), odontoblasts, bone and cartilage. With respect to the PNS, we employed anti-neurofilament immunostaining (Dodd et al., 1988) at E10.5 to visualize the cranial ganglia and associated cranial nerves. These studies demonstrated that the gross development of the cranial ganglia, including the trigeminal ganglia that are partially derived from the cranial NC (Begbie and Graham, 2001; Hamburger, 1961), occurred similarly in DCM embryos when compared with controls (Fig. 2G,H). Although we did note subtle alterations in the distal nerve projections at E10.5 and slightly hypoplastic trigeminal ganglion by E17.5 (Fig. S4G,H), we postulate that these differences were secondary to morphological changes in the face of DCM embryos. Furthermore, comparative Hematoxylin and Eosin analysis of the mandible at E17.5 revealed that DCM embryos and their wild-type counterparts were similar in the histological appearance of various NC-derived skeletal elements, including the lower incisor (Fig. 2I,J), the dentary bone and Meckel's cartilage (Fig. 2K,L).

In sum, the cranial NCC properties of specification, migration, proliferation and differentiation were relatively normal in the absence of AP-2α and AP-2β, except for a limited increase in apoptosis of post-migratory NCCs within the hinge domain. Except for this small domain of apoptosis, our studies nonetheless reveal that cranial NCCs from DCM embryos can differentiate into major derivatives, including cranial ganglia, odontoblasts, bone and cartilage. The finding that many facets of NCC biology can still occur in DCM embryos is surprising considering the substantial role AP-2 paralogs play in NCC specification in non-mammalian species (de Croze et al., 2011; Hoffman et al., 2007; Hong et al., 2014; Li and Cornell, 2007; Luo et al., 2003). Nevertheless, craniofacial development is severely disrupted in the DCM embryos (Fig. 1), suggesting that the crucial function of these AP-2 TFs is to cooperatively control within the NCCs the GRNs that are involved in morphogenesis, growth and patterning of craniofacial structures.

Previous studies have derived a series of interconnected GRNs that are necessary for the appropriate formation and patterning of the vertebrate mandible (Chai and Maxson, 2006; Clouthier et al., 2010; Parada and Chai, 2015). Therefore, given the cleft and hypoplastic mandible, we focused our remaining studies on how the combined function of AP-2α and AP-2β impacted the GRNs required for development of this BA1 derivative. Specifically, we conducted RNA profiling (RNA-seq) of the MdP in control versus DCM embryos at E10.5, a period when significant patterning is occurring and morphological changes are becoming more apparent in DCM embryos (see Fig. 2D). Triplicate biological replicates of individual MdPs were isolated from three control and three DCM embryos (Fig. 3A). After sequencing, two separate bioinformatic pipelines were implemented to identify differentially expressed genes (Fig. S5). The use of two pipelines was employed to reduce false-positives and false-negatives inherent to different bioinformatic approaches (Fonseca et al., 2014). These analyses resulted in either 279 (class I) or 673 (class II) genes showing 1.25-fold or greater differential expression between control and DCM MdP samples (Fig. 3B, Table S2). Overlapping these two datasets identified 182 ‘high-confidence’ differentially expressed genes (class III; 89 down, 93 up) in DCM embryos compared with controls (Fig. 3C, Fig. S6 and Table S2), and this smaller dataset showed the most robust correlation coefficient between fold-change values compared with the individual pipelines (Fig. S5).

We next used these gene lists to perform several bioinformatic analyses to assess the types of genes enriched and their relationship to known AP-2 targets and GRNs. First, we used both functional annotation clustering (FAC) (Huang da et al., 2009a,b) as well as geneset enrichment analyses (GSEAs) (Mootha et al., 2003; Subramanian et al., 2005) to initiate a systems-level analysis and extract over-represented molecular pathways from our datasets of altered gene expression. Using the 182 class III geneset, FAC identified the number one and two over-represented clusters to contain the terms ‘homeobox transcription factors’ (Fig. 3D, cluster 1) and ‘transcription’ (Fig. 3D, cluster 2), respectively. These two terms were also highly enriched when the analysis was run using the larger class I and II geneset (Fig. 3D), and were also significantly enriched in all three geneset using an alternative clustering software package [Enrichr (Chen et al., 2013; Kuleshov et al., 2016), data not shown]. Interestingly, all six Dlx genes, which are central to BA1 patterning, were present in the class III geneset and all six were downregulated in the DCMs (Fig. 3C, Fig. S6, Table S2). Notably, Dlx genes occur as bi-gene clusters in the mouse genome and these loci also produce opposite strand transcripts from Dlx1, Dlx4 and Dlx6. These Dlx1as, Dlx4os and Dlx6os1 transcripts were also represented in the class III dataset and were downregulated similarly to their adjacent Dlx genes. Therefore, of all the genes that were expressed at a significantly lower level in the DCMs, 10% mapped to the three Dlx gene clusters. Also in relation to the DLX-GRN, we observed modulation of genes associated with the activity of the homeobox TF Six1, which is notable as SIX1 has recently been shown to antagonize Dlx gene expression and regulate development of the intermediate ‘hinge’ region (Tavares et al., 2017). Thus, in class III geneset from the DCMs, we recorded increased expression of Six1, its partner Eya1, and potential co-factors Sox2, and Tlx1, but a reduction in the levels of Aes, a co-factor that can act to repress SIX1 function (Ahmed et al., 2012; Bajoghli et al., 2005; Riddiford and Schlosser, 2016) (Fig. S6 and Table S2).

Next, using the entire gene expression dataset, we implemented GSEAs with a variety of defined ‘geneset’ (Table S3) – based on a priori inferences about our dataset – along with a Q-value ranked list of most-altered genes in DCM embryos versus controls. Supporting FAC analysis, homeobox TFs were again found to be significantly enriched within our most-altered transcripts (Fig. 3E, Homeobox TFs). In addition, using a previously published microarray of altered gene expression in the MdP of Dlx5^−/− Dlx6^−/− double mutants (Jeong et al., 2008), we found significant enrichment of a similar set of genes to that discovered in DCM embryos versus controls (Fig. 3E, Dlx5^−/− Dlx6^−/−). Importantly, no such enrichment was observed with randomly generated geneset, including those based on known expression in BA1 (Fig. 3E, Random). Thus, both FAC and GSEAs identified misregulation of key homeobox TFs, especially Dlx family members, within the MdP of DCM embryos – placing AP-2α and AP-2β at a key node within the GRN patterning BA1 NCCs.

We further examined the intersection of AP-2 protein binding, enhancer activity and NCC specification by using our list of genes misregulated in the DCM mandible to mine recently published datasets relating to NCC chromatin and transcriptional signatures. First, we examined information on chromatin accessibility in migrating mouse NCCs (Minoux et al., 2017) and determined that cis-regulatory elements for the majority of the 182 genes we had identified were marked with either active (H3K4me2) or poised (H3K4me2/H3K27me3) chromatin marks in E10.5 mandibular NCC populations (Fig. 3F). Second, we analyzed a dataset derived from human induced NCCs (Rada-Iglesias et al., 2012) that examined the overlap between enhancer chromatin marks and the presence of AP-2-binding sites. Here, again, we found that a majority (∼78%, 143 of 182 genes) of the gene list we had generated had an associated AP-2-binding peak, with a large subset (∼77%, 110 of 143 peaks) of these found within active or poised enhancer marks (Fig. 3G). Two additional transcriptional regulators, NR2F1 and NR2F2, that co-occupied a subset of these human NCC enhancers with AP-2α (Rada-Iglesias et al., 2012), were also significantly downregulated in our 182 geneset. As previous studies have demonstrated that knockdown of Nr2f1 in NCCs decreases both TFAP2A and NR2F2 levels, (Rada-Iglesias et al., 2012), these findings suggest that the AP-2 and NR2F protein families may act in a regulatory circuit to mark active NCC enhancers. Taken together, examination of the gene expression differences between control and DCM mouse mandibles, in combination with mining of relevant enhancer datasets, reinforces the contention that the AP-2 TFs form a crucial aspect of the GRN that is required for cranial NCC development.

We extended the gene expression analysis by conducting whole-mount in situ hybridization on control and DCM embryos focusing on both downregulated (Dlx1, Dlx3, Dlx6, Emx2, Hand1) and upregulated (Six1, Fgf8) genes, many of which are part of a positional GRN that patterns the jaws. First, examining the ‘DLX code’, we assessed expression of one paralog from each cis first-order linked pair (Dlx1/2, Dlx3/4 and Dlx5/6). These bi-genic Dlx clusters are normally expressed in a nested fashion in the E10.5 mesenchyme (Fig. 4A,E,I and Fig. S1F,G). Dlx1 expression is found throughout the paired jaws, as well as where they intersect at the hinge region (Fig. 4A,C). Dlx6 expression is restricted to the hinge region and the MdP (Fig. 4E,G). Finally, Dlx3 expression is even more restricted to aboral regions of the distal region of the MdP (Fig. 4I,K). In contrast, at E10.5 in DCM embryos, Dlx1 expression was significantly reduced throughout both the MdP and MxP, including the hinge domain (Fig. 4B,D). Dlx6 expression in DCM embryos was also drastically reduced throughout most of the MdP and the adjacent hinge-region, although limited expression was detected at more distal locations (Fig. 4F,H). Finally, expression of Dlx3, which is a target of DLX5/6 activity, was also reduced, with expression nearly absent throughout the MdP (Fig. 4J,L). Regions of Dlx expression in BA2 mesenchyme, or for Dlx3 in the BA1 ectodermal domain, were less affected. Consistent with alteration to the DLX code in the NCC component of BA1, as well as the agreement between the GSEA profiles of DCM and Dlx5^−/−Dlx6^−/− datasets, additional downstream targets of DLX activity (Depew et al., 2002; Jeong et al., 2008) were reduced in the MdP of DCM embryos when compared with controls [e.g. Alx3, Alx4, Cited1, Col8a2, Dlx4, Gbx2, Hand1 and Hand2, see Table S2 and Jeong et al. (2008)]. For example, Hand1, which is expressed at the distal midline of the converging MdPs in controls at E9.5 and E10.5 (Fig. 4M,O), lacked this expression domain completely in DCMs (Fig. 4N,P). Moreover, expression of the homeobox TF Emx2, which was specifically detected in the hinge region in controls (Fig. 4Q), was abolished in the DCMs (Fig. 4R) whereas expression within the telencephalon and in BA2 remained relatively intact. In contrast to the loss of Emx2 signal and the downregulation of Dlx transcripts, the expression of Six1 expanded in the DCMs throughout BA1, including the hinge domain (Fig. 4, compare S with T). In addition, Fgf8, a signaling molecule normally expressed in more oral regions of the MdP (Fig. 4U) was expanded aborally in DCM embryos (Fig. 4V), again indicating significant changes in the GRN underlying jaw patterning.

Next, to address whether these gene expression changes also occurred earlier during NCC migration into BA1, we conducted quantitative real-time PCR (qRT-PCR) on E9.5 BA1 RNA of control and DCM embryos (Fig. 4W). These experiments revealed that, although some genes (e.g. Rspo2) that showed altered expression in E10.5 DCM embryos were unchanged compared with controls at E9.5, other crucial regulators, including Six1 and Dlx-family members were still significantly upregulated or downregulated, respectively at this earlier timepoint (Fig. 4W). Note that these early changes in the BA1 GRN presage the increase in cell death noted in the DCM hinge region at E10.5 that correlates with the loss of the Emx2 expression domain. In summary, RNA-seq, in situ hybridization and qRT-PCR, all highlight an important function of AP-2α and AP-2β activity in regulating expression of multiple genes that provide positional information to the developing jaws. Importantly, this patterning role for AP-2α/β is likely operational at the initial stages of establishing the proximodistal axis.

We next performed a detailed analysis of the Tfap2a/Tfap2b-Wnt1:CRE skeletons. The mammalian cranial skeleton has viscero/splanchnocranial, neurocranial and chondrocranial components. Whereas the first, which comprises the bones of the facial skeleton, is NCC derived, the neurocranial and chondrocranial components are derived from both NCCs (largely anterior) and mesoderm (largely posterior) (Couly et al., 1993; Gross and Hanken, 2008; Jiang et al., 2002; Le Douarin and Kalcheim, 1999; McBratney-Owen et al., 2008). Gross examination of E18.5 skeletons revealed that, as with the external morphology, mutant phenotypes ranged from mild in SCM-B skeletons, intermediate in SCM-A skeletons, to severe in DCM skeletons (Fig. 5 and Fig. S7). SCM-B skeletons (Fig. S7) were largely indistinguishable from controls (Fig. 5A,D), whereas SCM-A embryos had noticeable shortening of the premaxilla and mandible, as well as defects in the palatine bones consistent with cleft palate (Fig. 5B,E). DCM embryos had major defects throughout the NCC-derived skeleton, including the mandible, maxilla and premaxilla (Fig. 5C,F). In contrast, more posterior, mesoderm-derived elements in the skull, such as the basioccipital and exoccipital bones, developed normally in all mutant embryos (Fig. 5A-F). Overall, these observations are consistent with a redundant and cell-autonomous role for AP-2α and AP-2β in craniofacial skeletal development from the NC.

Strikingly, the analysis of the SCM-A and DCM craniofacial skeletons (Fig. 5 and 6) also revealed remarkable similarity to defects previously reported in mice with certain Dlx mutant combinations (Depew et al., 2005). We describe the most notable phenotypic overlaps below, concentrating on the SCM-A and DCM animals, and summarize our comprehensive findings in Fig. S8. We first examined the upper- and lower-jaw in greater detail, including the maxilla and mandible, as well as bones of the zygomatic arch (ZGA), which are partly derived from BA1 regions associated with the hinge domain (Fig. 5G-Q, Table S1B). In general, Dlx1 and Dlx2 have important functions with respect to the upper jaw and ZGA, Dlx5 and Dlx6 are more crucial for lower-jaw development, whereas loss of unlinked Dlx paralogs (e.g. Dlx2 and Dlx5) affect structures at the intersection of the upper and lower jaw (Depew et al., 2005; Jeong et al., 2008). The ZGA normally consists of the zygomatic process of the maxilla, the jugal and the zygomatic process of the squamosal (Fig. 5G). The latter process was often absent in SCM-A skeletons (Fig. 5H,K, asterisk in H) whereas DCM embryos had loss or transformation of the squamosal and jugal as well as a highly truncated thickened zygomatic process of the maxilla (Fig. 5I,L). With respect to the lower jaw, the most severe defects were observed in the proximal domain (Fig. 5M-Q). SCM-A skeletons displayed hypoplastic development of the coronoid, angular and condylar processes (Fig. 5N), elements that are exquisitely sensitive to reduced Dlx dose (Depew et al., 2005; Jeong et al., 2008). DCM mandibles had complete ablation of these proximal processes and/or fusion with components of the ‘hinge’ and upper jaw (syngnathia) (Fig. 5O,O′, Fig. S9A-C). Fusions included loss of the primary and secondary jaw joint articulations that varied among and within (left to right) DCM skeletons. Thus, some DCM embryos showed a nearly complete loss of patterning identity for the structures involved (Fig. 5O′), whereas others had apparent fusions between squamosal, maxillary and mandibular elements (Fig. 5O, Fig. S9B,C). DCM mandibles also showed defects along the mediolateral axis, with ectopic ossifications sometimes observed projecting orally from the mis-patterned mandible (Fig. 5P,Q). Jaw regions more distal to the hinge, such as the frontal process, the foramen of the maxilla (Fig. 5K,L) and lower incisors of the mandible, were less affected (Fig. 5M-O). Thus, a modest shortening of the distal mandible was apparent in SCM-A (Fig. 5N) and this defect was slightly more severe in DCM skeletons (Fig. 5O,O′).

We next examined additional skeletal elements derived from the BA1 hinge-associated regions that form the NC-derived chondrocranium and neurocranium (i.e. cranial base and cranial sidewalls) (Fig. 6). These regions are susceptible to ablations and ectopias in Dlx loss-of-function mutants (Depew et al., 2005; Qiu et al., 1997, 1995). We first studied E15.5 cartilage preparations to determine how chondrocranial elements, such as Meckel's cartilage, the malleus and the ala temporalis were affected by Tfap2 mutations (Fig. 6A-F, Fig. S10). When compared with controls, SCM-A embryos displayed a shorter Meckel's cartilage that failed to make a continuous connection at the proximal end with the malleus (Fig. S10B,E). Defects were more severe in DCM embryos with Meckel's cartilage that was about one half of its normal length, not fused at the midline, often deviating in its projection and failing to connect to a definitive malleus (Fig. 6C,F, Fig. S10C,F). The ala temporalis, the development of which is sensitive to Dlx dose, is a derived palatoquadrate element that, together with the lamina obturans, eventually forms the alisphenoid bone (Cerny et al., 2004; Depew et al., 2002). In SCM-A embryos, this cartilage was malformed and lacked the foramen (Fig. 6E). In DCM embryos, the ala temporalis failed to form, replaced by rudimentary cartilage elements that appeared to articulate with – or be non-continuous pieces of – Meckel's cartilage (Fig. 6F). Consistent with these earlier cartilage defects, E18.5 SCM-A and DCM skeletal preparations displayed a markedly affected alisphenoid bone (Fig. 6G-I). Thus, the alisphenoid bone was smaller in the SCM-A embryos than in controls and lacked the normal foramina (Fig. 6G,H). Interestingly, however, in a significant fraction of these SCM-A mice (∼40%, Table S1B), the adjacent squamosal bone had developed foramina, a feature absent from controls. In DCMs, the alisphenoid was severely mis-shapen and hypoplastic, and the squamosal bones were missing or severely transformed (Fig. 6I, also see Fig. 5I,O).

Dlx mutants exhibit ectopic skeletal struts and ectopic lateral cartilage and bony elements associated with the cranial base (Depew et al., 2005; Qiu et al., 1997, 1995). Examination of the E15.5 DCM chondrocranium also revealed that these had ectopic cartilage struts attaching to the cranial base, accompanied by loss of the alicochlear commissure (Fig. 6F). Although such ectopic struts were not seen in the E15.5 SCM-A preparations, by E18.5 they were present in both SCM-A and DCM skeletons (Fig. 6K,L, Fig. S9D-G). Like the reported ectopias in some Dlx loss-of-function mutants (Depew et al., 1999, 2005; Qiu et al., 1997, 1995), these struts varied in appearance between SCM-A and DCM skeletons but were never observed in controls (Fig. 6J). Specifically, SCM-A embryos generally displayed a single ossified strut extending laterally from the basisphenoid and often fusing with the gonial (Fig. 6K). In contrast, DCM skeletons had multiple and more irregular ectopic struts extending from the cranial base, including some with apparent intermittent cartilaginous ‘joints’ (Fig. 6L, Fig. S9D-G). In addition, unidentifiable cartilage and bony elements of variable morphology, reminiscent of the ectopic ‘palatoquadrate-like’ structure observed in Dlx mutants (Depew et al., 2005), were noted at the lateral aspects of DCM skulls (yellow arrowheads in Fig. 5I, also see Fig. S9H-J).

Additional defects in derivatives of BA1 and more caudal BAs, such as the bones and cartilages of the ear, hyoid, and thyroid, again mirrored those seen in Dlx mutants (Fig. 6M-T). With regard to the ear, these defects were most severe in DCMs, with loss of the gonial and tympanic bones (Fig. 6I,L,N), disrupted projection of the styloid process (Fig. 6N) and loss or transformation of the cartilage primordia of the middle ear bones (Fig. 6N). A range of severity was also apparent with respect to the hyoid and its association with the cranial base (Fig. 6O-T). Unlike control mice (Fig. 6O,R), both SCM-A and DCM embryos exhibited cleft or ‘pinched’ hyoids that were often fused to the thyroid cartilage, as well as displaying abnormal pterygoids that were positioned more rostrally (Fig. 6P,Q,S,T). Moreover, in DCMs, the abnormal hyoid-thyroid cartilage was also fused to both the pterygoids and the basisphenoid (Fig. 6Q,T).

Taken together, these findings indicate that both the Tfap2 and Dlx gene families show a dose effect on the craniofacial skeleton that is dependent on the number and type of mutant alleles present. The loss of AP-2α and AP-2β in the NC does not recapitulate the phenotype of the more extreme Dlx loss-of-function allelic combinations (e.g. Dlx5^−/− Dlx6^−/−), but instead mimics intermediate Dlx mutant phenotypes that are consistent with the observed ∼50% reduction in Dlx gene expression seen in DCMs (Fig. 3, Table S2). Furthermore, the similarity in skeletal defects between the Dlx and Tfap2 mutants extended to some external features, including low-set ears, microtia and microglossia (Fig. 1). Interestingly, for some Dlx loss-of-function mutant combinations, mandibular clefting was often associated with the homeotic transformation of lower-jaw structures into more ‘upper-jaw-like’ structures, including the emergence of ectopic vibrissae on the lower jaw (Beverdam et al., 2002; Depew et al., 2002, 2005). In this respect, the presence of structures with the superficial appearance of rows of ectopic vibrissae on the lower jaw of E14.5 DCM embryos is suggestive of a possible homeotic transformation (Fig. S11A-C). Similarly, the presence of foramina in the SCM-A squamosal bone may indicate a partial transformation towards the normal morphology of the adjacent alisphenoid bone (Fig. 6H). However, it was not possible to assess this feature in DCMs as these bones are missing or severely transformed (Fig. 6I). Overall, though, in agreement with the molecular analysis, morphological features in DCMs are consistent with altered Dlx signaling resulting in a major disruption to the BA1-MxP/MdP ‘hinge-associated’ structures. Moreover, these observations indicate that regulation of Dlx gene expression represents a major function of AP-2α/β in the BA1 NCCs.

Here, we have shown that AP-2α and AP-2β function cooperatively, and cell autonomously, within the cranial NCCs to orchestrate craniofacial development (see Fig. 7). Previous studies have shown that loss of each protein alone in the NCC results in limited defects (Brewer et al., 2004; Martino et al., 2016). Here, we demonstrate that pathology worsens as an additional allele is deleted from the remaining paralog, with loss of Tfap2a alleles showing greater effects than removal of Tfap2b alleles. Ultimately, deletion of both paralogs simultaneously using Wnt1:CRE results in major craniofacial defects that impact the BAs and the midface. These findings are consistent with the two AP-2 paralogs exerting cooperative and redundant roles in face development, particularly impacting late cranial NCC development. Such findings help resolve the apparent paradox that each individual AP-2 paralog has a limited role in mammalian NCC development (Auman et al., 2002; Brewer et al., 2004; Martino et al., 2016), whereas in other gnathostomes AP-2 is postulated to act as a conserved master-regulator of the NCC lineage (Green et al., 2015; Hong et al., 2014; Rada-Iglesias et al., 2012, 2013; Simoes-Costa and Bronner, 2015; Trainor, 2014). These findings also have substantial bearing on the interpretation of results derived from the study of gene expression in human and chimp NCC-derivatives. Genome-wide analysis of active enhancers in such NCC lines demonstrated an over-representation of AP-2-binding site sequences, and further indicated using ChIP that AP-2 proteins were detected at many of these sites. Additionally, nucleotide changes within AP-2-binding sites between human and chimp genomes often reduced AP-2 binding, with concomitant loss of H3K27ac levels, and reduced gene expression levels. Overall, these primate studies suggested that the AP-2 proteins were central to NCC enhancer function and that changes in AP-2 binding may drive craniofacial morphology differences between species (Prescott et al., 2015; Rada-Iglesias et al., 2012). Although these studies focus mainly on AP-2α, it is now clear that AP-2 function in NCC development should be considered in the context of all expressed paralogs. For example, although the over-representation of AP-2 binding at human NCC enhancers was attributed to AP-2α, AP-2β shares an identical DNA-binding motif (Mohibullah et al., 1999; Williams and Tjian, 1991a,b; Zhao et al., 2001) and many antibodies used for ChIP-Seq detect not only AP-2α, but also other AP-2 proteins. Thus, these sites – and associated chromatin features and attributed influences on primate craniofacial morphology – likely represent coordinated AP-2α and AP-2β targets. Indeed, RNAseq analysis of these human NCC cultures identified both TFAP2A and TFAP2B to be highly enriched relative to the embryonic stem cells and neurospheres from which they were derived (Rada-Iglesias et al., 2012), whereas TFAP2C to TFAP2E were either expressed at low levels or not at all, in a pattern highly reminiscent of the E8.5 mouse NCCs (Minoux et al., 2017). Thus, at least in vitro for human and in vivo for mouse, NCCs share common AP-2 paralog expression. Furthermore, the presence of AP-2α/AP-2β appears to be required for important aspects of NCC GRN expression in both mouse and primates. Such findings are important in the context of dominantly inherited human branchio-oculo-facial and Char syndromes as mutations in one allele of TFAP2A or TFAP2B, respectively, could interfere with global AP-2 function in the NCCs via a dominant-negative mechanism of action (Li et al., 2013; Milunsky et al., 2008; Satoda et al., 2000).

The conclusion that combinations of AP-2 proteins, notably AP-2α and AP-2β, are important in NCC development is also supported by the concerted action of these proteins in the generation of trunk NCC derivatives, including the enteric nervous system and melanocytes (Seberg et al., 2017). Given these observations, it will be important to determine whether these two AP-2 TFs are also crucial in other NCC populations that contribute to tissues such as the heart and adrenal gland. A deeper understanding of the molecular function of AP-2α and AP-2β in these additional NCC derivatives will establish whether regulation of key developmental genes, including homeobox TFs, is a conserved function of AP-2 in the NC or whether this mechanism is unique to the NCCs of BA1. Similarly, identifying whether AP-2α and AP-2β regulate survival of a subpopulation of post-migratory cells, in this case approximating to the hinge domain, in additional derivatives outside the craniofacial region will be an important facet of further understanding the cellular functions of AP-2.

Despite major NCC-derivative defects in the head and trunk of DCM embryos, early NCC specification in this mouse model appears unaffected. This is interesting as loss of AP-2 paralogs (singly or redundantly) in non-mammalian vertebrates, including Xenopus and zebrafish, disrupts NCC specification (de Croze et al., 2011; Hoffman et al., 2007; Hong et al., 2014; Li and Cornell, 2007; Luo et al., 2003). Several possibilities exist to explain why NCC specification is largely intact in DCM embryos. First, the Wnt1:CRE driver may mediate gene deletion post-NCC specification. However, we do not favor this model as preliminary data show that Sox2-CRE-generated DCM embryos, which lose both AP-2 genes in the early epiblast prior to NCC specification, still undergo NCC development (E.V.O., unpublished). Second, additional paralogs, particularly Tfap2c, may also function redundantly with Tfap2a/Tfap2b in the early stages of NCC specification such that loss of all three paralogs would be required to disrupt NCC formation. Third, AP-2 paralogs may not only act autonomously within the NCCs but also non-cell-autonomously, particularly through their expression in the embryonic ectoderm, with both modes impacting overall NCC formation and function (Knight et al., 2005; Meulemans and Bronner-Fraser, 2002). Finally, AP-2 proteins may have distinct roles in mammalian versus non-mammalian NCC development – a feature that has previously been posited for additional gene families in the mouse (Barriga et al., 2015). Future experiments will be necessary to distinguish between these models.

Mice with NCC-specific loss of Tfap2a/Tfap2b expression add to a small group of models that present with hinge defects, including syngnathia. Prominent among these are loss-of-function mutations in Dlx gene family that normally function to subdivide BA1 into upper- and lower-jaw identities (Acampora et al., 1999; Beverdam et al., 2002; Depew et al., 1999, 2002, 2005; Jeong et al., 2012, 2008; Qiu et al., 1997, 1995). Manipulation of the endothelin (EDN) signaling pathway also impacts the development of the upper and lower jaws, with the Dlx genes acting as crucial targets (Clouthier et al., 2010). Recent studies have further shown that reduced expression of Fgf8 or Foxc1 can cause this pathology (Inman et al., 2013), as can overexpression of BMP4 in NCCs (He et al., 2014). However, both Fgf8 and Foxc1 are upregulated in the DCM samples, while the expression of BMP pathway genes associated with facial development are not significantly altered in the DCM model (Table S2 and Fig. S12). Instead, molecular profiling of the MdP indicated that loss of AP-2α and AP-2β activity resulted in alterations of multiple transcriptional regulators, including the entire Dlx gene family. Notably, this reduction in Dlx gene expression occurred without changes in the mRNA levels of EDN pathway components (e.g. Edn1, Ece, Plc, Ednra, Mef2c; Table S2). Such findings imply that AP-2 functions to regulate Dlx expression through a node that is independent or downstream of EDN signaling in the NC.

How then might AP-2α/β cooperatively regulate craniofacial patterning autonomously within the NCCs to drive proper hinge formation, ultimately facilitating upper- and lower-jaw registration? One possibility is that AP-2 regulates the survival of a small population of cells within BA1, such that NCC-specific loss of Tfap2a/Tfap2b induces apoptosis at the intersection of the MdP and MxP. Loss of these BA1 hinge domain cells, which would normally express Emx2, would then subsequently result in syngnathia development. At the same time, the considerable overlap between AP-2 and DLX loss-of-function phenotypes suggests that there is also a more global impact on the expression of homeobox TFs that extends beyond the limited cell death observed in DCM embryos. This conclusion is also supported by the reduction of Dlx expression observed in DCM embryos, as well as the overlap between DCM and Dlx5^−/−Dlx6^−/− MdP mRNA expression profiles. Theoretically, such additional regulatory mechanisms can be split into direct (Fig. 7D) or indirect (Fig. 7E) modes of action. First, for the direct scenario, AP-2α and/or AP-2β would bind directly to critical cis-elements associated with relevant genes to either activate or repress transcription within this BA1 GRN. In support of this direct model, homologs of at least 143 of the 182 altered genes in DCMs have an associated AP-2-bound region in the human in vitro-derived NCCs, the majority of which (77%) are also associated with active or poised histone marks (Rada-Iglesias et al., 2012). The involvement of AP-2 with these cis-acting sequences may be crucial for NCCs to respond appropriately to external signals as they enter the facial prominences (Minoux et al., 2017). Various indirect models can also be proposed in which AP-2α/β may modulate the transcript levels of an intermediary TF or a chromatin effector, or alternatively act post-transcriptionally to alter the properties of a protein complex required for gene expression (Fig. 7E, ‘Z’). Ultimately, both direct and indirect mechanisms are likely involved in coordinating the GRNs responsible for BA1 development. Teasing apart these potential mechanistic roles of AP-2α/β will be a major priority in further understanding their function in normal craniofacial development, as well as how their aberrant expression results in human birth defects.

As previously stated (Van Otterloo et al., 2016), this study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Institutional Animal Care and Use Committee of the University of Colorado – Denver. Noon on the day that a copulatory plug was present was denoted as embryonic day 0.5. For the majority of experiments, littermate embryos were used when comparing between genotypes. Yolk sacs or tail clips were used for genotyping. DNA for PCR was extracted using DirectPCR Lysis Reagent (Viagen Biotech) plus 10 µg/ml proteinase K (Roche), incubated overnight at 65°C, followed by heat inactivation at 85°C for 45 min. Samples were then used directly for PCR-based genotyping with primers (Table S4) at a final concentration of 200 nM using the Qiagen DNA polymerase kit, including the optional Q Buffer solution (Qiagen).

Sox2:CRE (Hayashi et al., 2002), Wnt1:CRE (Danielian et al., 1998) and Gt(ROSA)26Sor^tm1Sor/J [Rosa26 reporter/R26r/; (Soriano, 1999)] mice were obtained from the Jackson Laboratory. The Tfap2a alleles used in this study have been previously described, including the Tfap2a-null (Zhang et al., 1996), Tfap2a^tm2Will/J [Tfap2a floxed conditional; (Brewer et al., 2004)] and Tfap2a^tm1Hsv/J [Tfap2a-lacZ knock-in; (Brewer et al., 2002)] alleles. Once established, lines were maintained on a Black Swiss (strain code 492, Charles River) background.

To generate the Tfap2b alleles used in this study, a P1 genomic clone from 129 DNA corresponding to the Tfap2b locus was isolated by library screening and subsequently an ∼1.65 kb BamHI fragment containing exon 5 and flanking sequences was isolated. This clone was extended by 70 nucleotides at the 5′ end using a PCR-based strategy to introduce an NheI site ∼340 bp upstream of exon 5, and by ∼6 kb at the 3′ end by the addition of a BamHI-HindIII fragment from the P1 clone that extends across exon 6 and into exon 7. An Acc65I site 70 nt upstream of exon 6 was repaired with T4 DNA polymerase and used to introduce a LoxP FRT neo FRT cassette derived from a modified version of plasmid PK11 (Lewandoski et al., 1997) containing a LoxP at the 5′ end rather than at the 3′ end. A second LoxP site was introduced at a PshAI site ∼50 nt downstream of exon 6, such that an ∼260 bp fragment containing exon 6, which encodes a vital part of the DNA-binding and dimerization domain of AP-2β, is flanked by LoxP sites. Two copies of the thymidine kinase selection cassette were then introduced downstream of the HindIII site, generating a targeting vector that was verified by sequence analysis (Fig. S3). The construct was then linearized with NheI prior to electroporation into 129S1/SV W9.5 ES cells. Clones showing appropriate recombination at both the 5′ and 3′ end were identified using long-range PCR and Southern blotting using primer sequences or probes outside the region of homology. Three positive clones were identified and these were karyotyped. Subsequently, a euploid ES clone was used to generate chimeras, and these were bred with Black Swiss mice to obtain germline transmission of the Tfap2b-neo allele. Next, the neo cassette was removed by breeding to (ACTFLPe)9205Dym/J (Jackson Labs) to generate the Tfap2b-flox allele (Fig. S3). The Tfap2b-null allele was then generated by breeding with Tmem163^{Tg(ACTB-cre)2Mrt}/J mice (Jackson Labs). Successful modification of offspring was confirmed by PCR analysis and mice containing these alleles have been maintained for multiple generations on a Black Swiss background.

The following genetic scheme was generally implemented to generate embryos for experiments. The male adult, which was Tfap2a^null/wt;Tfap2b^null/wt;CRE+ (Wnt1:CRE for example) was crossed with a female that was, Tfap2a^flox/flox;Tfap2b^flox/flox (and R26r homozygous for appropriate experiments). This cross would generate SCM-A, SCM-B and DCM, each at an ∼1:8 ratio (12.5%), with the remainder being various combinations of alleles serving as controls. Note, SCM embryos would be either Tfap2a^null/flox;Tfap2b^wt/flox;CRE+ (SCM-A) or Tfap2a^wt/flox;Tfap2b^null/flox;CRE+ (SCM-B), and DCM embryos would be Tfap2a^null/flox;Tfap2b^null/flox;CRE+. By virtue of this breeding scheme, all SCMs would be conditionally heterozygous for the alternate paralog, as described in the Results. The exception to the males listed above included experiments using the Tfap2a-lacZ knock-in allele (Brewer et al., 2002), in which case males used for breeding were Tfap2a^KI/wt;Tfap2b^null/wt;CRE+.

For western blot analysis, E12.5 whole heads (Sox2:CRE experiments) or E10.5 craniofacial prominences (Wnt1:CRE experiments) were lysed, on ice, in RIPA buffer containing protease inhibitors (ThermoScientific Protease Inhibitor Cocktail). Briefly, as described previously (Van Otterloo et al., 2016), tissue was first minced with a disposable pestle followed by further disruption using a tissue homogenizer (Pro200 homogenizer, PRO Scientific). Following homogenization, samples were allowed to lyse on ice for ∼30 min. After lysis, 6× Laemmli Buffer (plus β-mercaptoethanol) was added to a 1× final concentration and samples were stored at −80°C until use. On the day of use, samples were boiled at ∼100°C for 10 min, spun at 16,000 g for 10 min, and then ∼20 µl of lysate loaded per lane on a 12% stacking SDS polyacrylamide gel and separated at 100 V. Once samples were resolved, they were transferred onto a nitrocellulose membrane overnight at 4°C and 40 mA. Following transfer, membranes were blocked for 1 h with 3% powdered milk in TBST (TBSTM), and then incubated with the primary antibody, diluted in TBSTM, overnight at 4°C (rabbit anti-TFAP2 polyclonal, H-87/sc-8976, 1:1000, Santa Cruz Biotechnology). Following primary antibody incubation, membranes were washed four times for 20 min with TBSTM and then incubated with an infrared (IR)-labeled (goat anti-rabbit) secondary antibody for 1 h (1:20,000 in TBST, LI-COR Biosciences). The membrane was then washed four times for 20 min with TBSTM, for 5 min with PBS and then imaged on a LI-COR Odyssey imager. Serving as a loading control, post-imaging, the membrane was again blocked with TBSTM, four times for 20 min, and an additional primary antibody was added to the membrane (anti-GAPDH, mouse monoclonal G8795 at 1:10,000 from Sigma-Aldrich or rabbit anti-ACTIN polyclonal, H-196 from Santa Cruz Biotechnology at 1:1000) and a similar washing, secondary antibody and imaging procedure carried out, using the appropriate species-specific secondary antibody.

Embryos were collected at appropriate time points and processed using Alizarin Red and Alcian Blue stains as previously described (Van Otterloo et al., 2016).

Embryos were collected at appropriate time points and processed using Alcian Blue stain as previously described (Van Otterloo et al., 2016).

β-Galactosidase (β-gal) staining of whole embryos was carried out as previously described (Seberg et al., 2017). Briefly, dissected samples were fixed for 30 min in 0.2% glutaraldehyde (in PBS) at room temperature. Following fixation, embryos were washed three times for 30 min with a ‘lacZ’ rinse buffer (see Seberg et al., 2017) at room temperature. Following washes, embryos were developed in a ‘lacZ staining’ solution (‘lacZ’ rinse solution plus 1 mg/ml X-Gal, Invitrogen/Life Technologies) in the dark at 37°C until desirable β-gal signal was produced (∼1-2 h for R26r×Wnt1:CRE experiments). Post-staining, embryos were briefly rinsed multiple times in PBS and then post-fixed in 4% paraformaldehyde overnight at 4°C, while protected from light.

Whole-mount immunohistochemistry was carried out essentially as described previously (Pabst et al., 2003), replacing 10% horse serum with 5% milk as a blocking reagent. The anti-neurofilament antibody (Dodd et al., 1988) was used at a 1:100 dilution (IgG clone 2H3, obtained from the Developmental Studies Hybridoma Bank, created by the NICHD of the NIH and maintained at The University of Iowa, Department of Biology, Iowa City, IA, USA). Hematoxylin and Eosin (H&E) staining was carried out as previously described (Cardiff et al., 2014; Van Otterloo et al., 2016). In situ hybridization was carried out as previously described (Simmons et al., 2014; Van Otterloo et al., 2016).

For analysis of cell proliferation and cell death, embryos were collected at various developmental time points and fixed overnight in 4% PFA at 4°C. Following fixation, tissue was dehydrated in a graded series of ethanol and xylene, and subsequently embedded in paraffin wax. After embedding, samples were sectioned with a Leica RM 2235 microtome at 10-12 µm onto charged glass slides. After drying, sections were processed for immunofluorescence. Briefly, paraffin wax was washed from tissue sections with three 10 min washes in xylene, followed by rehydration in a graded series of ethanol and water. Following rehydration, samples were washed in PBS for 5 min at room temperature. Subsequently, an antigen retrieval procedure was carried out, which involved placing samples in 10 µM citric acid (pH 6.0) and boiling samples in a Cuisine Art electric pressure cooker, on high for 5 min. Following antigen retrieval, tissue was washed once in TBST, followed by blocking for 1 h in TBSTM. Sections were then incubated overnight in primary antibody, diluted 1:250 in TBSTM at 4°C in a humidified chamber. Primary antibodies included anti-p-Histone H3 (sc-8656-R, Santa Cruz Biotechnology, rabbit polyclonal) or anti-cleaved caspase 3 (5A1E, Cell Signaling Technology, rabbit monoclonal), for cell proliferation or death, respectively. Following primary antibody incubation, samples were washed twice for 10 min in TBST at room temperature, followed by a 30 min wash in TBSTM. Samples were then incubated for 1 h with a secondary antibody (goat anti-rabbit IgG, Alexa Flour 488 conjugate, ThermoFisher Scientific/Invitrogen, R37116) and DRAQ5 (Abcam, ab108410) nuclear stain, diluted 1:250 and 1:5000, respectively, in TBSTM. Processed samples were imaged on a Leica TCS SP5 II confocal microscope and individual images taken for visualization.

For the MdP RNA-seq, E10.5 DCM and litter-matched control embryos were dissected in ice-cold PBS. Subsequently, the mandibular prominences were carefully dissected from each embryo and stored in RNAlater (Ambion/Life Technologies) until later use. Following genotyping, RNA was extracted, essentially as previously described (Van Otterloo et al., 2016), from paired (i.e. left and right) mandibular prominences (three control and three DCM embryos; three biological replicates) using the microRNA Purification Kit (Norgen Biotek) and following manufacturer's protocol. Following elution, mRNA was further purified using the Qiagen RNAeasy Kit according to the manufacturer's protocol. The quality of extracted mRNA was assessed using DNA Analysis ScreenTape (Agilent Technologies) to ensure that it was of sufficient quality for library production. Following validation of extracted mRNA, cDNA libraries were generated using the Illumina TruSeq Stranded mRNA Sample Prep Kit. Following library generation and subsequent quality control assessment, cDNA was sequenced using the Illumina HiSeq2500 platform and single-end reads (1×150).

1. Pipeline 1: For an initial ‘pipeline’ to identify mRNA-based gene expression changes between control and DCM mandibular prominences, derived fastq sequences were analyzed by applying a custom computational workflow (Baird et al., 2014; Bradford et al., 2015; Henderson et al., 2015; Maycotte et al., 2015) consisting of the open-source gSNAP (Wu and Nacu, 2010), Cufflinks (Trapnell et al., 2010), and R for sequence alignment and ascertainment of differential gene expression. In short, reads generated were mapped to the mouse genome (Mm10) by gSNAP, expression (FPKM) derived by Cufflinks and differential expression analyzed with ANOVA in R.

2. Pipeline 2: For a secondary ‘pipeline’ to further refine mRNA-based gene expression changes between control and DCM mandibular prominences, fastq sequence reads were first trimmed using the Java software package Trim Galore! (Babraham Bioinformatics, Babraham Institute, Cambridge, UK) and subsequently mapped to the Mm10 genome using the HISAT2 software package (Pertea et al., 2016) (both with default settings). Following mapping, RNA expression levels were generated using StringTie (Pertea et al., 2016) and differential expression computed between genotypes using CuffDiff2 (Trapnell et al., 2013), with a significance cut-off value of Q<0.05 (FDR-corrected P-value).

Genes identified as significantly altered in their expression values that were present in both bioinformatic pipelines (class III, ‘182 geneset’) were used for clustering analysis based on their expression values calculated with pipeline 1 or pipeline 2. Hierarchical clustering and heatmap generation based on these gene expression changes were produced using the pheatmap R-package (written by Raivo Kolde, cran.r-project.org/package=pheatmap).

For functional annotation clustering of significantly upregulated and downregulated genes, the recently updated DAVID (v6.8) platform (Dennis et al., 2003), with default parameters, was used. Following functional annotation clustering, ‘annotation clusters’ above a twofold enrichment threshold were considered. Gene lists from both bioinformatic pipelines 1 and 2 were used, along with the overlapping 182 geneset. In addition to DAVID, the online pathway analysis suite, Enrichr (Chen et al., 2013; Kuleshov et al., 2016) (amp.pharm.mssm.edu/Enrichr/), was used with the available gene lists generated.

For geneset enrichment analysis (GSEA), the GSEA application (Subramanian et al., 2005) was downloaded and implemented from the Broad Institute website (software.broadinstitute.org/gsea/index.jsp). First, using gene expression and Q-values calculated from bioinformatic pipeline 1, a ranked gene dataset file was generated (.rnk file). Subsequently, a curated dataset that included multiple gene lists was generated (Table S3). Lists included, genes upregulated, downregulated, or both, in the MdP of Dlx5^−/−Dlx6^−/− mutants (Jeong et al., 2008) as well as a curated list of homeobox transcription factors (Holland et al., 2007). The ‘Dlx-associated’ lists were calculated by re-analyzing data deposited in the GEO repository (accession number GSE4774) using the GEO2R software package (Davis and Meltzer, 2007) and retaining genes that showed over a 1.5-fold change between control and Dlx5^-/−Dlx6^−/− mutant samples, and P<0.002.

For E8.5 pre-migratory cranial NCCs, pre-processed RNAseq datasets (Minoux et al., 2017) were downloaded from the GEO repository (accession number GSE89434, file GSE89434_Minoux_RNAseq_gene_counts.csv) and average normalized mRNA expression values were calculated from the three biological replicate values listed. For E10.5-E12.5 cranial NCCs, facial prominences (MdP, MxP and FNP) were micro-dissected from E10.5, E11.5 and E12.5 embryos, which were generated from in-crossing C57BL/6J wild-type adults, and superficial ectodermal layers were separated from underlying mesenchyme, as described previously (Li and Williams, 2013). mRNA was extracted from each individual tissue source, and cDNA was generated and sequenced following procedures listed above. mRNA expression levels were calculated, as described above, for ‘Pipeline 1’. A more-detailed and comprehensive analysis of the generation and analysis of the E10.5-E12.5 datasets will be presented elsewhere (H.L. and T.W., unpublished) and these datasets are now available through the Facebase Consortium website (www.facebase.org) with the accession number FB00000867.

To generate expression percentiles for the various RNAseq datasets, a normalized ranked gene list was calculated for each. Briefly, after calculation of RNA expression values, the entire dataset was sorted based on expression levels and each gene ascribed a ranking within the dataset (i.e. rank of 1, being the gene with the highest mRNA expression, rank of 2, the next highest, etc.). mRNAs with an RPKM of 0 would all correspond to an equal, and lowest, ranking within the dataset. Next, to normalize between RNAseq experiments, the rank value of each gene was then divided by the total number of genes within the dataset, generating a ‘normalized-rank value’. Data were subsequently plotted inversely (so that the most highly expressed genes were near the top of the y axis), and location of Tfap2 paralogs within the normalized ranked list were highlighted.

E10.5 MdP histone datasets were obtained from the GEO repository (accession number GSE89435). Raw fastq sequencing reads for H3K4me2, H3K27me3 and input were downloaded and reprocessed, essentially as described previously (Minoux et al., 2017). Briefly, raw reads were mapped to the mouse GRCm38/mm10 reference genome using the QuasR Bioconductor package (Gaidatzis et al., 2015), implementing the qAlign function. Subsequently, the qCount function was used to assess the number of reads per promoter region (transcription start site, ±1 kb). ChIP enrichments were calculated using the formula described (Minoux et al., 2017), and resulting data graphed using Microsoft Excel 2016. Subsequent to generating a scatterplot of H3K4me2-H3K27me3 enrichments, the 182 class III geneset was mapped on top of these data points in Excel.

To identify genes associated with TFAP2 ChIP-Seq data points, the previously published, pre-processed, ChIP-Seq datasets (Rada-Iglesias et al., 2012) were downloaded from the GEO repository (accession number GSE28876, files NCC_TFAP2A, NCC_H3K4me3 and NCC_H3K27ac) in bed file format. Subsequently, datasets were uploaded to the University of California Santa Cruz Genome Browser (UCSC genome browser) and the Table Browser function (genome.ucsc.edu/goldenpath/help/hgTablesHelp.html) used to transfer data to the online GREAT algorithm, version 3.0.0 (bejerano.stanford.edu/great/public/html/index.php) (McLean et al., 2010). Using default parameters, GREAT was implemented to generate a list of genes associated with AP-2 ChIP-Seq peaks (i.e. AP-2 binding) or AP-2 ChIP-Seq peaks that occurred within genomic regions containing active histone marks (e.g. H3K4me3 or H3K27ac). Overlapping of TFAP2A ChIP-Seq and histone mark datasets was carried out with the ‘intersection’ function in Table Browser prior to GREAT analysis. After gene lists were generated, the VLOOKUP function in Excel was used to identify common genes between the 182 significantly altered genes in DCM embryos and the GREAT-generated gene lists.

For real-time PCR experiments, embryos were harvested at E9.5 and BA1 dissected off for RNA isolation. Tissue was stored in RNAlater at −20°C until genotyping was completed on samples. Three biologically independent control and DCM embryos were collected. Following positive identification of genotypes, tissue was equilibrated at 4°C for 1 day, RNAlater removed and RNA extracted from tissue samples using the RNeasy Mini Kit (Qiagen). Following RNA isolation and quantification, cDNA was generated using 200 ng/µl of RNA and the SuperScript III First-Strand Synthesis Kit (Invitrogen/ThermoFisher Scientific), along with the optional on-column DNase digestion. Once cDNA was generated, quantitative real-time PCR analysis was conducted using a Bio-Rad CFX Connect instrument, Sybr Select Master Mix (Applied Biosystems, ThermoFisher Scientific) and 20 µl reactions (all reactions performed in triplicate). All primers (sequences are in Table S5) were designed to target exons flanking (when available) large intronic sequences. Relative mRNA expression levels were quantified using the ^ΔΔCt method (Dussault and Pouliot, 2006) and an internal relative control (β-actin). Significance (standard Student's t-test) was calculated using a previously published Statistical Analysis Software (SAS) package (Yuan et al., 2006).

We thank Dr Timothy Nottoli (Yale University) for targeted ES-cell injections in generation of Tfap2b mice, Drs Sonia Leach and Joan Hooper (University of Colorado – Denver) for additional bioinformatic support, Irene Choi and Dr Rebecca Green (University of Calgary) for assistance with the mouse lines, the DSHB at the University of Iowa for antibodies used in these studies, the Facebase consortium, and the University of Colorado – Denver Microarray and Genomics Core for RNA sequencing. In addition, we thank the many members of the Department of Craniofacial Biology at the University of Colorado – Denver who provided critical review of initial drafts of this manuscript, particularly Drs Katherine Fantauzzo, Kristin Artinger and Andre Tavares.