MOZ directs the distal-less homeobox gene expression program during craniofacial development

Craniofacial development is a complex process involving the outgrowth of five distinct prominences: the frontonasal prominence, the paired maxillary processes and the paired mandibular processes. These processes need to grow in a concerted manner at a specific rate and then fuse in the midline to form normal facial structures (reviewed by Dixon et al., 2011; Lan et al., 2015). A reduction in growth, an increase or decrease in cell death (e.g. Grabow et al., 2018; Ke et al., 2018) or a failure of midline fusion can result in craniofacial anomalies, including small upper and lower jaws, misalignment of teeth, cleft lip and/or palate.

The palate divides the nasal and oral cavities, and so separates the airway from the mouth. Its development commences with the formation of the primary palate from the medial nasal prominences (reviewed by Bush and Jiang, 2012; Gritli-Linde, 2008). The secondary palate arises from the maxillary component of the first pharyngeal arch, develops as horizontal outgrowths on the oral side of the maxillary processes [in mice at embryonic day (E) 11.0] and continues with vertical growth down alongside the tongue (at E13.5). The palatal shelves elevate (at E14.5) until they lie horizontally above the tongue and then grow horizontally and fuse in the midline (Gritli-Linde, 2007; Walker and Fraser, 1956; Yu and Ornitz, 2011).

Failure of secondary palate development leads to cleft palate. Oral clefts, including cleft lip and cleft of the secondary palate, occur at a rate of 1 in 500-2500 live births, causing problems with feeding, speaking, hearing and social integration (Dixon et al., 2011). Phenotypic concordance of oral cleft in monozygotic twins has been reported as 31% (Grosen et al., 2011), indicative of genetic as well as non-genetic contributions to the cause of the disorder. Chromatin modifiers have been proposed to mediate environmental effects on gene expression. A small number of genes encoding chromatin-modifying or chromatin-associated proteins has been associated with cleft palate, suggesting that they may be involved in epigenetic or chromatin-mediated mechanisms that affect palate development. The proteins include MEN1 (Engleka et al., 2007), which is associated with the MLL/SET1 H3K4 methyltransferase complex, the H3K9me1/2 methyltransferase PHF8 (Fortschegger et al., 2010; Laumonnier et al., 2005; Loenarz et al., 2010), the chromodomain helicase DNA binding protein CHD7 (Sperry et al., 2014), the H3K4me and H3K9me demethylase KDM1A (Tunovic et al., 2014), the H3K27me demethylase KDM6A (Lindgren et al., 2013) and the lysine acetyltransferase KAT6B (Campeau et al., 2012; Clayton-Smith et al., 2011). Although their biochemical function is known, their molecular effects on palate development remain unclear except in the case of MEN1 (Engleka et al., 2007).

Heterozygous mutations of the chromatin modifier gene KAT6A (also known as MOZ and MYST3) have been found to cause a congenital disorder characterised by craniofacial anomalies (including cleft palate at low penetrance), heart defects and developmental delay (Arboleda et al., 2015; Millan et al., 2016; Tham et al., 2015). KAT6A encodes the histone lysine acetyltransferase monocytic leukaemia zinc-finger protein (MOZ) and was originally identified as the target of recurrent chromosomal translocations causing an aggressive form of acute myeloid leukaemia (Borrow et al., 1996). Null mutation of the Moz (Kat6a/Myst3) gene in mice causes haematological anomalies, characterised by a complete absence of definitive haematopoietic stem cells (Katsumoto et al., 2006; Sheikh et al., 2016; Thomas et al., 2006), and specific defects in B-cell (Good-Jacobson et al., 2014; Sheikh et al., 2015b) and T-cell development (Newman et al., 2016). In addition, loss of MOZ causes developmental defects, including body segment identity, craniofacial and heart defects. Mice lacking MOZ have micrognathia, a complete secondary palate cleft (100% penetrance), an extensive anterior homeotic transformation affecting 19 body segments, as well as aortic arch and cardiac ventricular septum defects (Vanyai et al., 2015; Voss et al., 2009, 2012). Similarly, MOZ is required to determine pharyngeal segmental identity in zebrafish (Crump et al., 2006; Miller et al., 2004).

MOZ target genes relevant for the correct specification of body segment identity and for normal heart development have been identified. MOZ is required for normal levels of histone 3 lysine 9 acetylation (H3K9ac) and for transcription of Hox genes, which are crucial for correct specification of segment identity. Induction of Hox gene expression with retinoic acid rescues the segment identity defects in Moz^–/– mice, but not the craniofacial defects (Voss et al., 2009). Similarly, for cardiac septum development, MOZ is essential for normal levels of H3K9ac and transcription at the Tbx1 gene locus; transgenic overexpression of Tbx1 rescues cardiac septum defects in Moz^–/– mice, but not the craniofacial defects (Voss et al., 2012). Thus, while MOZ target genes required for body segment identity specification and heart development are known, the effects of MOZ on the transcriptional program in the complex process of craniofacial and palate development are unknown. Unlike many other aspects of facial dysmorphogenesis, oral clefts (cleft lip/cleft palate) are clearly defined birth defects, and so we concentrated on this particular aspect of the Moz mutant craniofacial defects. Here, we report that MOZ is an activator of gene expression, including the Dlx5 gene, and of the DLX transcription factor-induced gene expression program during craniofacial development.

Moz mRNA is expressed widely in most or all cells at all stages of mouse embryonic development examined (E7.5, E8.5, E9.5, E10.5 and E13.5), including in the maxillary component of the first pharyngeal arch, which gives rise to the secondary palate (Sheikh et al., 2015a; Thomas et al., 2006; Voss et al., 2009). We have previously reported craniofacial defects in MOZ-deficient mouse foetuses, including micrognathia and a complete cleft of the secondary palate (Voss et al., 2012).

Frontal serial sections of the heads prior to (E13.5), during (E14.5) and immediately after (E15.5) palatal shelf elevation revealed that the palatal shelves of Moz^+/+ littermate control and Moz exon 16 truncated Moz^Δ/Δ mutants (no protein detectable, therefore null) occupied a similar position at E13.5 (Fig. 1A,B), but that the Moz^Δ/Δ shelves had slightly reduced cross-sections of the intermediate region (105,985±2266 µm² versus 115,543±1569 µm², P=0.026, n=3 animals per genotype). At E14.5, the palatal shelves of wild-type foetuses had completed elevation and fused along the midline (Fig. 1C). In contrast, the palatal shelves of the Moz^Δ/Δ mutants were smaller and had not elevated to a horizontal orientation above the tongue (Fig. 1D), indicating a defect in growth and elevation of the shelves at this stage. At E15.5, the palatal shelves of Moz^Δ/Δ foetuses had elevated, but remained apart, resulting in a wide cleft (compare Fig. 1F with E,G).

Cell death within the palatal shelves was not affected by the absence of MOZ (Fig. 1H). Assessment of BrdU incorporation at rostral, intermediate and caudal levels revealed a modest decrease of ∼15% in cells in S-phase of the cell cycle at the intermediate level, where palate fusion first occurs (P=0.017; Fig. 1G,I). Similarly, Moz-deleted E10.5 pharyngeal arch cells displayed a reduction in the percentage of cycling cells to ∼65% of controls, as assessed by flow cytometric analyses (P=0.03; not shown).

Moz^–/– pups have a range of craniofacial anomalies affecting the palate, the upper and the lower jaw. Cleft palate occurs with 100% penetrance (Voss et al., 2012). Development of the secondary palate begins with the outgrowth of palatal shelves as horizontal ridges on the oral side of the maxillary part of the first pharyngeal arch between E10.5 and E11.5 in mouse development. We wanted to examine the effects of loss of MOZ on gene expression before palate development was initiated to avoid possible confounding effects of already changed tissue and cell type compositions in the samples, reasoning that any effects of MOZ on gene expression should be detectable prior to its effects on phenotypic outcome. We therefore began our investigation at E10.5.

To enrich for facial, maxillary and mandibular structures [while at the same time performing RNA-seq on a tissue type that would also provide enough material for chromatin immunoprecipitation (ChIP) experiments], we dissected the maxillary and mandibular component of the first pharyngeal arch, as well as the second pharyngeal arch at E10.5, and used these for RNA-seq expression profiling (Fig. 2A). Among the more than 500 genes that were downregulated in the Moz^–/– 1st and 2nd pharyngeal arches (P values 0.002 to 6×10⁻²²; FDR<0.05; listed in Table S1) were 19 genes that had previously been associated with cleft palates [Barx1, Col2a1, Ctgf (Ccn2), Dlx1, Dlx2, Dlx3, Dlx5, Efnb2, Gbx2, Glce, Ift140, Osr2, Pkdcc, Pygo2, Shh, Sim2, Smoc1, Sox9 and Tshz1; Fig. 2B,C]. Among these, the distal-less homeobox (Dlx) gene family stood out as five members of this family were affected by the loss of MOZ (Dlx1, Dlx2, Dlx3, Dlx4, Dlx5; P values 0.003 to 6×10⁻¹³; FDR<0.05; Fig. 2B; Table S1). Dlx1, Dlx3, Dlx4 and Dlx5 were among the top 30 differentially expressed genes (Fig. 2D; P values and FDRs in Table S1). In addition to the Dlx genes, the gastrulation brain homeobox 2 (Gbx2) and odd-skipped related 2 (Ors2) genes were particularly affected by loss of MOZ (Fig. 2C). Dlx1, Dlx3, Dlx4, Dlx5 and Gbx2 mRNA levels in Moz^+/– heterozygous pharyngeal arches were intermediate between Moz^–/– and Moz^+/+ (Fig. S1), suggesting a tight regulation of these genes by MOZ. Six genes that had previously been associated with cleft palates were upregulated in the Moz^–/– pharyngeal arches (P values 0.001 to 10⁻⁵; FDR<0.05; Fig. 2E; Table S1). Among these, the runt-related transcription factor 2 gene (Runx2) was most prominently affected by loss of MOZ (Fig. 2E).

Based on their prominence among the top 30 differentially expressed genes, we investigated the effects of loss of MOZ on the Dlx gene family and downstream target genes of DLX transcription factors further.

The effects of loss of MOZ on Dlx gene family mRNA and antisense RNA levels were confirmed by RT-qPCR on independent E10.5 1st and 2nd pharyngeal arch samples (Fig. 2F). Together, the RNA-seq and RT-qPCR data showed that Dlx1, Dlx2, Dlx3, Dlx4 and Dlx5 gene expression was lower in Moz^–/– embryos compared with controls. However, it did not distinguish whether the mRNA levels were lower because the expression levels were lower or because the expression domains were smaller. To distinguish these two possibilities, we conducted whole-mount in situ hybridisation with a shorter (dynamic range; Fig. 2G) and a longer (endpoint; Fig. S2) staining to reveal expression levels and full expression domains, respectively. We observed that Dlx gene expression was impaired in both ways, the expression levels of Dlx3 and Dlx5 were lower (Fig. 2G), and the expression domains of Dlx1, Dlx2, Dlx3, Dlx5 and Dlx6 appeared reduced in size (Fig. S2).

To investigate the effects of MOZ for the duration of palate development, we dissected the maxillary component of the first pharyngeal arch, at E10.5 and E11.5, as well as the palatal shelves at E13.5 and E14.0, and determined Dlx gene expression levels by RT-qPCR. In addition, we conducted RNA-seq experiments on E13.5 palatal shelves. Overall, Dlx gene expression was lower in the Moz^–/– tissues than in controls (Fig. S3). These data indicate that normal levels of mRNA expression of the Dlx1, Dlx2, Dlx3 and Dlx5 genes depend on MOZ from the earliest stages of craniofacial and palate development, and that the importance of MOZ for Dlx1 and Dlx2 gene expression persists throughout palate development (Figs S1 and S3; Tables S1 and S2).

Homozygous mutation of Dlx1 (Qiu et al., 1997), Dlx2 (Qiu et al., 1995) or Dlx5 (Acampora et al., 1999; Depew et al., 1999) causes cleft palate with different levels of penetrance; the highest penetrance results from loss of Dlx5 on a mixed genetic background, which at 88% is similar to the penetrance of cleft palate in Moz-null pups on a mixed background (100%). Moreover, Dlx5 gene expression was strongly affected by the absence of MOZ (Fig. 2B,F). To determine whether the Dlx5 gene was likely to be a direct MOZ target, we examined MOZ occupancy and histone acetylation by chromatin immunoprecipitation (ChIP) followed by qPCR or sequencing. ChIP-qPCR of microdissected pharyngeal arches derived from embryos homozygous for a Flag-V5-BIO tag fused to the C-terminal end of MOZ (Moz^V5/V5; Vanyai et al., 2015) revealed enrichment for the MOZ-V5 signal in Moz^V5/V5 at two sites within the Dlx5 gene compared with wild-type controls (Fig. 3A,B). Histone 3 lysine 9 acetylation (H3K9ac) was reduced in the absence of MOZ (Fig. 3C). In contrast, H4K14ac was unaffected by Moz gene status (Fig. 3D). Furthermore, ChIP-seq experiments suggested further areas of enrichment of MOZ-V5 in the Dlx5/Dlx6/Dlx6os locus (Fig. 3E). ChIP-qPCR experiments confirmed three of these regions of enrichment (Fig. 3F). In comparison with the Dlx5 locus, the Dlx1, Dlx2, Dlx3 and Dlx4 loci did not show similarly strong MOZ occupancy by ChIP-seq (not shown). However, owing to the small amount of input material (1st and 2nd pharyngeal arches), we cannot exclude the possibility that MOZ also occupies other Dlx genes at lower levels.

A number of DLX transcription factor target genes, as well as regulatory pathways have been proposed as relevant to craniofacial development by previous studies (summarised in Fig. 4A). Expression of the zinc-finger transcription factor genes odd-skipped related 1 and 2 (Osr1 and Osr2), the BarH-like homeobox gene 1 (Barx1), the forkhead box transcription factor gene Foxl2, the basic helix-loop-helix (bHLH) transcription factor single-minded homolog 2 gene (Sim2) and the POU domain transcription factor gene Pou3f3 is reduced in Dlx2^–/– and/or Dlx1^–/–;Dlx2^–/– mutant embryos (Jeong et al., 2012, 2008). Similarly, mRNA levels of the goosecoid homeobox gene (Gsc) and Gbx2, as well as the heart and neural crest derivatives expressed 2 gene (Hand2) are decreased in Dlx5^–/– and/or Dlx5^–/–;Dlx6^–/– mutants (Depew et al., 1999, 2002; Jeong et al., 2008). Of these, Gbx2 and Hand2 have been confirmed as direct targets of DLX transcription factors (Barron et al., 2011; Jeong et al., 2008).

The expression of the DLX target genes Osr1, Osr2, Gbx2, Sim2 and Barx1 were significantly downregulated in the E10.5 Moz^–/– pharyngeal arches RNA-seq data compared with Moz^+/+ controls (FDR 0.03 to 9×10⁻⁶; Fig. S4; Table S1). Moreover, the expression of DLX transcription factor targets was also downregulated in Moz^–/– mutants compared with Moz^+/+ embryos as assessed by whole-mount in situ hybridisation (Hand2, Gbx2, Gsc, Osr2 and Sim2, Fig. 4; Osr1 and Pou3f3, Fig. S5). In contrast, pituitary homeobox 2 (PITX2) has been implicated in the regulation of Dlx2, i.e. acting upstream of Dlx gene expression (Green et al., 2001; Venugopalan et al., 2011). Consistent with this, Pitx2 expression in the pharyngeal arches was not affected by the loss of MOZ (Fig. S5). Similarly, other genes reported to operate upstream of Dlx gene expression or at the same level were not differentially expressed (Fig. S4).

These data indicate that, along with Dlx gene expression, the mRNA levels of genes downstream of DLX transcription factors were also downregulated in the absence of MOZ, whereas upstream factors were not affected.

The palate is formed by intramembranous ossification from mesenchymal condensations (Bush and Jiang, 2012). During this process, mesenchymal cells under the influence of the transcription factors paired box 1 and 9 (PAX1 and PAX9), as well as NK3 homeobox 1 and 2 (NKX3-1 and NKX3-2) give rise to skeletal precursor cells that, in the presence of RUNX2, form osteoprogenitor cells and then osteoblasts (reviewed by Hartmann, 2009; Lefebvre and Bhattaram, 2010). In contrast, SRY-box 9 (SOX9) induces the formation of chondroblasts in regions of endochondral ossification from cartilage precursors. RUNX2 (Komori et al., 1997; Otto et al., 1997) and PAX9 (Peters et al., 1998) are essential for skeletal development. PAX9 induces the differentiation of mesenchymal cells to bipotential skeletal precursor cells, and RUNX2 supports differentiation towards the osteoblast lineage.

Apart from the downregulation of genes required for palate development, the Moz^–/– E10.5 pharyngeal arches displayed an upregulation of genes involved in osteoblast differentiation. Runx2 and Nkx3-1 were upregulated (Fig. S6; Table S1). Bone collagen genes Col6a1 and Col6a3 were upregulated in Moz^–/– E10.5 pharyngeal arches. Furthermore, the osteoblast-specific factor/periostin gene (Postn) was upregulated. In contrast, Sox9 and the SOX9 target genes encoding cartilage collagens, Col2a1 and Col9a1, were downregulated. A similar pattern was observed in the Moz^–/– E13.5 palatal shelves, where the bone collagen gene Col1a2 was upregulated (Fig. S6; Table S2), which is a RUNX2 target gene (Kern et al., 2001). At E13.5, additional genes promoting bone formation and encoding bone matrix proteins were upregulated in the Moz^–/– palatal shelves. These included genes encoding elements of the bone morphogenetic protein (BMP) and transforming growth factor β (TGFβ) signalling pathways [BMP5, ALK1/ACVRL1 (activin A receptor-like type 1), TGFB2], as well as the matrix-associated secreted protein acidic and cysteine-rich protein (SPARC; Fig. S6; Table S2). The upregulation of genes that promote osteogenesis and of bone collagen genes suggests that the Moz^–/– E10.5 pharyngeal arches and E13.5 palatal shelves may undergo premature osteoblast lineage differentiation.

The bulk of the cells expressing Dlx genes in the pharyngeal arches are neural crest cells. An ability of MOZ to directly activate one or more of the Dlx gene family loci would necessitate a cell-autonomous requirement for MOZ within the neural crest for Dlx gene expression. To test this requirement, we used a Wnt1 promoter-driven cre-recombinase transgene [Wnt1-cre^T (Danielian et al., 1998)] to delete the loxP-flanked exons 3 to 7 of the Moz conditional mutant allele (Moz^lox) in the neural ectoderm, which gives rise to the neural crest cells that migrate to the pharyngeal arches.

Deletion of Moz in the Wnt1-cre lineage (Moz^lox/lox;Wnt1-cre^T/+) resulted in neonatal lethality (data not shown), leading to underrepresentation at weaning (Fig. 5A). Gross examination, skeletal preparations and serial sectioning (Fig. 5B, Fig. S7, Table S3) of E18.5 heads revealed cleft of the soft palate in 10 of 10 Moz^lox/lox;Wnt1-cre^T/+ pups, which was never seen in control Moz^+/+;Wnt1-cre^T/+ animals. An abnormal additional bone caudal to the alisphenoid bone (similar to the os paradoxicum observed in Dlx5 null mice) was observed in 4 of 4 Moz^lox/lox;Wnt1-cre^T/+ animals examined for its presence.

All except one Moz^lox/lox;Wnt1-cre^T/+ pup displayed a normal external appearance. Only this one severely affected Moz^lox/lox;Wnt1-cre^T/+ pup displayed shortening of the upper and lower jaw (Fig. S7). RT-qPCR analysis of E10.5 pairs of microdissected first and second pharyngeal arches confirmed that loss of MOZ in the neural crest resulted in the downregulation of Dlx genes compared with control samples (Fig. 5C) to levels similar to those observed in Moz germline-deleted (Moz^–/–) pharyngeal arches (Fig. 2F).

Dlx5 homozygous mutant animals have numerous craniofacial defects (Acampora et al., 1999; Depew et al., 1999, 2005; Han et al., 2009), which appear to vary depending on the genetic background and include cleft of only the soft palate (Han et al., 2009) or the bony palate (88%) (Acampora et al., 1999; Depew et al., 1999), and asymmetric hypoplasia (88%) or near-complete agenesis of the nasal capsule and cavity (Depew et al., 1999). In addition to these defects, Dlx5^–/– mutants have an ectopic bone structure originating from the ectotympanic process and extending toward the pterygoid process, termed the os paradoxicum (Acampora et al., 1999; Depew et al., 1999).

As our results thus far suggested that the Dlx gene family requires MOZ for normal mRNA expression levels in E10.5 embryos and that the Dlx5 locus was directly activated by MOZ, genetic interaction between the Moz and Dlx5 genes was examined in the offspring of Moz^+/– by Dlx5^+/– matings. Moz^+/–; Dlx5^+/– animals were underrepresented at weaning (3 weeks of age; Fig. 6A), indicating synergy between Moz and Dlx5 heterozygosity. Analysis of pups revealed that 11 out of 17 Moz^+/–;Dlx5^+/– animals had a cleft of the soft palate, which was not observed in wild type or either single heterozygous animals (asterisk in Fig. 6F), but similar to the cleft soft palate produced through Wnt1-Cre mediated deletion of Moz (see overview in Table S3). Although clefts of the soft and/or bony palate were a common feature of Moz^–/–;Dlx5^–/– and Moz^+/–;Dlx5^+/– pups, Moz^+/–;Dlx5^+/– pups displayed certain similarities with Dlx5^–/– animals that were not present in Moz^–/– animals. Asymmetric hypoplasia of the nasal capsules was observed in Moz^+/–;Dlx5^+/– (arrows in Fig. 6H,I) and Dlx5^–/– pups, but not in Moz^–/– animals. Externally, a depression in the dorsal surface of the rostrum reminiscent of Dlx5^–/– pups was observed in three out of 14 Moz^+/–;Dlx5^+/– pups. This depression was not observed in wild-type, Moz^–/– or single heterozygous animals (arrows Fig. 6M,N).

Analysis of skeletal preparations of wild-type, Moz^+/–, Dlx5^+/–, Moz^+/–;Dlx5^+/–, Moz^–/– and Dlx5^–/– E18.5 pups revealed an ossified strut in all six Moz^+/–;Dlx5^+/– pups examined, albeit less extensive than the os paradoxicum observed in Dlx5^–/– mutants (indicated in red in Fig. 6S′,T′, Table S3). This ectopic bone was not observed in any wild-type or Moz single heterozygous animal, whereas one out of six Dlx5^+/– single heterozygous pups had a very small fragment of ectopic bone in approximately the same region (indicated in red in Fig. 6R′). Moz-deficient pups occasionally had a supernumerary bone structure in a similar position (indicated in blue in Fig. 6U′). Lastly, Moz^+/–;Dlx5^+/– animals, and animals that had lost one or both copies of the Dlx5 gene or the Moz gene displayed varying degrees of shortening of the mandible, as well as a shortening of individual bony protrusions from the mandible (Fig. S8).

In order to determine whether the synergistic effects of the combined loss of Moz and Dlx5 alleles extended to transcriptional regulation, we performed RNA-seq of microdissected first and second pharyngeal arches from individual stage-matched female E10.5 embryos of the following six genotypes: wild type, Moz^+/–, Dlx5^+/–, Moz^+/–;Dlx5^+/–, Moz^–/– and Dlx5^–/– (Fig. 7A). A greater-than-additive increase in the number of up- and downregulated genes in Moz^+/–;Dlx5^+/– tissues compared with the sum of each of the heterozygotes suggested synergy between the Moz and Dlx5 genes in transcriptional regulation (Fig. 7B). Dlx3 and Dlx4 mRNA levels were downregulated in animals that had lost one of two copies of Moz and/or Dlx5 (Fig. 7C; FDRs in Table S4; all differentially expressed genes in Table S5; Fig. S9). In Dlx5^–/– mutant mice, where exons 1 and 2 are deleted, the Dlx5 locus demonstrated substantial reads in the third exon, but negligible reads in exons 1 and 2 (Fig. S9). The other genotypes showed decreased numbers of reads in exons 1 and 2 to varying magnitudes when compared with wild type, with Moz^–/– and Moz^+/–;Dlx5^+/– showing the lowest number and Dlx5^+/– expression intermediate between wild type and these genotypes. Most Dlx gene mRNA levels were also reduced in the Dlx5^–/– and the Moz^+/–;Dlx5^+/– pharyngeal arches (Fig. S9; FDRs in Table S4).

A multiple dimension scaling (MDS) plot showed that the biological replicates of each genotype clustered together (Fig. 7D). In addition, wild-type and single heterozygous samples grouped together. Interestingly, Moz^+/–;Dlx5^+/– samples displayed expression profiles more similar to Dlx5^–/– mutants when compared with either single heterozygote, indicating that loss of one allele of Moz shifted the expression profile of Dlx5^+/– samples to become more Dlx5^–/– mutant like. Overall, the MDS plot suggests a progression in transcriptional distance from wild type to single heterozygous mutants to double heterozygous mutants and finally to Dlx5^–/– homozygous mutants, with the Moz^–/– homozygous mutant expression profiles being separate and unlike the other genotypes (Fig. 7D).

A heatmap of the top 30 differentially expressed genes in Moz^+/–;Dlx5^+/– samples compared with wild type showed that the Moz^+/–;Dlx5^+/– resembled the Dlx5^–/– patterns more closely than the Moz^–/– and wild-type patterns (Fig. 7E). Correspondingly, a heatmap of the top 30 differentially expressed genes in Dlx5^–/– samples compared to wild type also showed the similarity between the Moz^+/–;Dlx5^+/– pattern and the Dlx5^–/– samples (Fig. S9). In contrast, a heatmap of the top 30 differentially expressed genes in Moz^–/– samples compared with wild type displayed a pattern for Moz^–/– that was distinct from the other genotypes, as expected because MOZ directly regulates more genes than Dlx5 alone (Figs 2, 7B; Tables S1 and S2). Even so, many genes in Moz^+/–;Dlx5^+/– and Dlx5^–/– samples showed changes in the same direction as the Moz^–/– samples (Fig. S9).

The confirmed direct DLX transcription factor target gene Gbx2 was downregulated in the Moz and/or Dlx5 mutant genotypes (Fig. 7F; FDRs in Table S4). In contrast, the mRNA levels of upstream regulators of Dlx gene expression were not significantly affected, including the endothelin receptor type A (Ednra), fibroblast growth factor 8 (Fgf8), myocyte enhancer factor 2C (Mef2c) and Pitx2 genes (FDRs in Table S4).

The wingless-like protein 6 (Wnt6), a gene encoding a growth factor promoting the self-renewal of palatal mesenchymal cells (Jiang et al., 2017), was downregulated in Moz and/or Dlx5 mutant genotypes compared with the wild-type controls (Fig. 7G; FDRs in Table S4). In contrast, bone development-inducing genes were upregulated in the Moz^+/–;Dlx5^+/– and the Dlx5^–/– samples compared with wild-type controls (Fig. 7H; Fig. S9). In particular, Runx2 and Pax9 mRNA were elevated in Moz and/or Dlx5 mutant genotypes compared with the wild-type controls (Fig. 7H; FDRs in Table S4). Moreover, genes encoding proteinaceous extracellular matrix, including collagens, were significantly upregulated in the Moz^+/–;Dlx5^+/– samples compared with wild-type controls (Fig. 7I). These included the bone-specific type I and type VI collagen genes [Col1a2 (a RUNX2 target gene; Kern et al., 2001) and Col6a3], and the osteoblast-specific factor periostin (Postn; Fig. 7I; FDRs in Table S4). In particular, with respect to expression changes in Runx2, Pax9, Col1a2, Col6a3 and Postn, loss of Moz and Dlx5 appeared to synergise (Fig. 7H,I; Table S4), in that single knockouts had a lesser effect than loss of one allele of both Moz and Dlx5. Overall, the expression profiles suggested premature onset of osteoblast differentiation in Dlx5, Moz and Dlx5;Moz compound mutant embryos.

Heterozygous mutation of the human MOZ (KAT6A) gene causes facial dysmorphogenesis and can result in palate defects (Arboleda et al., 2015; Millan et al., 2016; Tham et al., 2015) and homozygous deletion of Moz in mice causes cleft palate with 100% penetrance, indicating that MOZ is important for craniofacial development in humans and mice. In this study, we have identified the molecular mechanism of MOZ function during craniofacial development. Of a large number of genes affected by the lack of MOZ, we examined the distal-less homeobox genes, because seven of the top 30 differentially expressed genes were Dlx genes or their antisense RNAs, making them the most prominently affected gene family with reduced expression in the Moz-deleted embryos. Moreover, downstream target genes of DLX transcription factors were also affected. Genomic occupancy of the Dlx5/Dlx6 locus by MOZ, tissue-specific deletion in the neural crest cells and genetic interaction between Moz and Dlx5 indicated a direct effect of MOZ on Dlx5 gene expression. MOZ is one of only a few chromatin modifiers implicated in palate and craniofacial development, and therefore joins the ranks of candidates proposed to account for the integration of genetic and environmental effects on this process. In this context, it should be noted that Moz^+/– heterozygous mice are more susceptible to cleft palate caused by excessive supply of retinoic acid (Voss et al., 2012), suggesting that two wild-type alleles of Moz render palate development more robust.

Interestingly, the methyl-CpG-binding protein 2 has been found to recruit histone deacetylase 1 (HDAC1) to the Dlx5/Dlx6 locus, promoting the repressive mark H3K9 methylation, at least in the brain (Horike et al., 2005). If this was also the case in craniofacial structures, Dlx5/Dlx6 locus gene expression in craniofacial development may be regulated by the opposing functions of MOZ and HDAC1.

The DLX transcription factors are crucial for craniofacial development (Acampora et al., 1999; Depew et al., 1999, 2005; Jeong et al., 2012, 2008; Qiu et al., 1997, 1995). One-third of individuals with copy number variation at the DLX5/DLX6 locus present with mild to severe palatal abnormalities, such as high palate or cleft palate and cleft lip (Elliott and Evans, 2006), and a DLX4 mutation has been described in a child with bilateral cleft lip and palate (Wu et al., 2015). Single homozygous mutant mice for Dlx1, Dlx2 and Dlx5 display cleft palate with varying penetrance (Acampora et al., 1999; Depew et al., 1999, 2005; Han et al., 2009; Qiu et al., 1997, 1995).

Loss of MOZ causes a reduction in Dlx1, Dlx2, Dlx3, Dlx4, Dlx5 and Dlx6 mRNA levels by 25% to 55% (rather than complete loss). Compound heterozygous loss of Dlx1, Dlx2, Dlx3, Dlx5 and Dlx6 causes cleft palate with 100% penetrance (Depew et al., 2005). Although Depew and colleagues did not present the expression levels of the Dlx genes in the compound heterozygous mice, based on studies of other genes in which mRNA levels in heterozygous mice have been determined, it can be hypothesised that mRNA levels may be halved in the Dlx1, Dlx2, Dlx3, Dlx5 and Dlx6 compound heterozygotes. These levels are similar to those in the Moz-null embryos and represent limiting levels of Dlx gene mRNA. Although the effects of MOZ on the Dlx gene family alone could account for the complete penetrance cleft palate phenotype, direct effects of MOZ on other genes may also contribute to the cleft palate phenotype. The downregulation of 15 cleft palate genes in addition to the Dlx genes in the Moz-null embryos would be expected to contribute to the Moz-null phenotype.

Supporting a role for MOZ in promoting transcription of Dlx genes, we found that loss of MOZ also affected the expression of target genes of DLX transcription factors, Hand2 (Depew et al., 2002), Gsc (Depew et al., 1999), Gbx2 (Jeong et al., 2008), Osr1, Osr2 and Sim2 (Jeong et al., 2012). The fact that the Wnt1-cre deletion of Moz caused a less severe phenotype than germline deletion suggests a role for MOZ outside the neural crest, in addition to the pronounced effects on Dlx gene expression within the neural crest. A role for MOZ in the epithelium remains to be examined. In contrast, mesoderm-specific deletion of Moz does not result in cleft palate or other craniofacial anomalies (Vanyai et al., 2015).

Mesenchymal cells of neural crest or mesodermal origin (reviewed by Lefebvre and Bhattaram, 2010) give rise to skeletal mesenchymal precursor cells (also termed osteochondral progenitors) under the influence of PAX9, NKX3-1 and other transcription factors (reviewed by Hartmann, 2009). Skeletal precursors can undergo self-renewing divisions and give rise to either chondrocytes or osteoblasts directed by opposing functions of SOX9 or RUNX2, respectively (Zhou et al., 2006). While RUNX2 is required for the differentiation of both chondrocytes and osteoblasts (Ducy et al., 1997; Komori et al., 1997; Otto et al., 1997; Stricker et al., 2002), it is considered a ‘master regulator’ of intramembraneous ossification (Bhatt et al., 2013; Takarada et al., 2016).

Genetic deletion of Runx2 results in the complete failure of ossification (Komori et al., 1997; Otto et al., 1997); the combined absence of NKX3-1 and the related NKX3-2 causes skeletal anomalies, in particular in the cranial region (Herbrand et al., 2002), and the lack of PAX9 results in multiple skeletal defects, including cleft palate (Peters et al., 1998). Runx2, Nkx3-1, Pax9 and bone-specific downstream effector genes were upregulated in Moz^–/– single knockout, Dlx5^–/– single knockout and Moz^+/–;Dlx5^+/– double heterozygous embryos. Upregulated collagen genes, Col1a1, Col6a1 and Col6a3, encode structural proteins expressed during osteogenesis (Shoulders and Raines, 2009). Both neural crest- and mesoderm-derived cells produce a collagen-rich matrix in the early condensation state of skeletal development (reviewed by Lefebvre and Bhattaram, 2010). Several collagen genes, including Col1a2, Col3a1 (Diez-Roux et al., 2011) and Col6a2 (Reymond et al., 2002) are expressed during palate development. Mutations in human or mouse collagen genes, including COL1A2 (Garofalo et al., 1991) and COL2A1 (Snead and Yates, 1999), cause cleft palate. Collagen genes are directly regulated by SOX9 (Lefebvre et al., 1997) and DLX5 (Hojo et al., 2016), and by RUNX2 (Kern et al., 2001).

The downregulation of the gene encoding WNT6, which promotes palatal mesenchymal cell proliferation, and the modest reduction in cell proliferation combined with the upregulation of the bone differentiation-inducing genes Pax9 and Runx2, as well as the bone extracellular matrix genes Col1a2, Col6a3 and Postn in the Moz mutants and in the Moz and Dlx5 compound mutants suggest a gradual premature cell cycle deceleration and premature onset of osteogenic lineage differentiation. Premature cell cycle deceleration and premature onset of differentiation would be expected to result in undersized palatal shelves that cannot meet and fuse in the midline, as observed in the Moz^–/– embryos.

Cleft palate can occur as an intrinsic effect of palate development or as a consequence of other craniofacial anomalies. Moz-deficient mice display mild craniofacial defects, which may contribute to the cleft palate. The regulation of the Dlx5 locus by MOZ may direct attention to the derivatives of the mandibular arch and raise the issue of whether the cleft palate observed in Moz^–/– pups may be secondary to defects in lower jaw development, because severe shortening of the mandible, as seen in Dxl5^–/– foetuses, can affect tongue placement and, secondarily, palate development. However, Moz^–/– and Moz^+/–;Dlx5^+/– mandibles were only slightly shorter (∼10%) and were similar to Dlx5^+/– mandibles, which do not have palate clefting. Therefore, the only slightly shortened Moz^–/– mandibles cannot explain the interesting relationship between MOZ activation of the Dlx5 locus and clefting of the palate. We propose that the combined effect of loss of MOZ on the expression of the entire Dlx gene family and other genes results in cleft palate and other craniofacial anomalies observed in the Moz^–/– mice. In addition, a number of findings suggest a palate intrinsic effect of MOZ. Moz single homozygous mutant animals displayed a delay in palatal shelf elevation and a reduction in horizontal growth, while surrounding structures appeared comparably normal. This was accompanied by a small reduction in BrdU-positive cells in the neural crest-derived mesenchymal cells of palatal shelves of about 15% at a rostrocaudal level, where the palatal shelves first meet in the midline and commence fusion. Although 15% may appear to be a modest reduction in cells in S phase, compounded over only four cell cycles, a 15% reduction in cells in each S phase results in a ∼50% reduction in tissue size [(1-0.15)⁴=0.52]. These findings suggest that intrinsic defects in mesenchymal cell proliferation contribute to the cleft palate phenotype. In conclusion, we have identified changes in the molecular program underlying craniofacial development that manifest in the absence of the histone acetyltransferase MOZ, providing evidence for a molecular mechanism that potentially causes the craniofacial anomalies seen in individuals with the KAT6A (MOZ) heterozygous mutation syndrome.

All animal experiments (Mus musculus) were approved by the Walter and Eliza Hall Institute Animal Ethics Committee and conducted in accordance with the Australian code of practice for the care and use of animals for scientific purposes. For timed matings, midday following the detection of a vaginal plug in the morning was designated embryonic day (E) 0.5. Four previously described Moz alleles were used for this study: (1) exon 16 truncated Moz^Δ mice (no protein detectable, therefore null) on a mixed FVB×BALB/c background (Thomas et al., 2006), and (2) exons 3 to 7 loxP-flanked Moz^lox mice (null allele; Voss et al., 2009), (3) exons 3 to 7 deleted Moz^– mice (Voss et al., 2009) and (4) Flag-V5-Bio 3′-tagged Moz^V5 mice (Sheikh et al., 2015c; Vanyai et al., 2015) maintained on a C57BL/6 background. Depending on the genetic background, animals homozygous for the Moz^– allele die between E15.5 and birth with 100% penetrance due to developmental defects that are individually lethal, namely cleft palate, a ventricular septum defect and an interrupted aortic arch (Voss et al., 2012). In contrast, Moz^V5/V5 mice are viable, healthy and fertile, indicating that the Moz^V5 allele produces functional protein. C57BL/6 Moz^–/– were used for mid-gestational molecular work, while Moz^–/– on a mixed genetic background, as well as FVB×BALB/c Moz^Δ/Δ animals were used to characterise morphological defects at late gestation. The Moz^– allele on a mixed genetic background was also used in whole-mount in situ hybridisation and RT-qPCR experiments. Moz^– mice on an inbred C57BL/6 background were used for the genetic interaction experiments with Dlx5^– mice (Depew et al., 1999), also on the C57BL/6 background. Moz^lox mice on a C57BL/6 background were used for neural crest-specific deletion using Wnt1-cre^T/+ transgenic mice (Danielian et al., 1998), also on a C57BL/6 background.

Skeletal preparations were performed as previously described (Thomas et al., 2000). Whole embryo heads were fixed overnight in 4% paraformaldehyde (E13.5-15.5) or for 3 days in Bouin's fixative (E18.5) before serial sectioning and Hematoxylin and Eosin staining according to standard histological protocols.

Whole-mount in situ hybridisation was conducted using standard protocols (Thomas et al., 2007) using sense and antisense cDNA probes detailed in Table S6. qPCR was performed on ChIP genomic DNA or cDNA following reverse transcription (RT) using Superscript III (Life Technologies) on the LightCycler 480 (Roche) using SYBR-green technology (SensiMix, Bioline). qPCR primers are listed in Table S7.

Primary murine pharyngeal arch cells were isolated in the following manner. E10.5 first and second pharyngeal arches were dissected from individual embryos, transferred into 25 µl of 0.125% trypsin/EDTA solution and incubated for 10 min on ice, then for 2 min at 37°C. Culture medium [175 µl neural stem cell medium (Merson et al., 2006) with 5% ESC-qualified FBS] was added. The pharyngeal arches were dissociated by gentle pipetting. Cells were centrifuged at 200 g and resuspended in culture medium.

For cell cycle analysis by fluorescence-activated flow cytometry, cells were fixed in 500 μl of cytofix/cytoperm solution (BD Pharmingen) for 30 min on ice, washed in Perm/Wash solution (BD Pharmingen) and then incubated for 5 min in Perm/Wash solution. Cells were centrifuged, resuspended in residual buffer and split into two: one aliquot was incubated with 20 μl anti-Ki67 (FITC-conjugated anti-Ki-67, BD Pharmingen, 556026) and the other was incubated with isotype control antibody at 4°C overnight. Cells were washed twice in KDS-BSS +2% FCS buffer and resuspended in 200 μl with a 1 μg/ml final concentration of DAPI in deionised H₂O and incubated at room temperature for 30 min. Cells were washed twice in KDS +2% FCS buffer then sieved into a single cell suspension and analysed using a flow cytometer (LSRIIC Becton, Dickinson and Company).

E13.5 pregnant mice were injected with a single dose of 100 mg/kg BrdU intraperitoneally and embryos recovered 1 h later. BrdU-treated E13.5 heads were fixed in 4% paraformaldehyde, embedded in paraffin wax and serially sectioned. Frontal sections of rostral, intermediate and caudal palatal shelves were processed for anti-BrdU immunostaining and detection of genomic DNA fragmentation by TUNEL staining on adjacent sections as described previously (Thomas et al., 2000). BrdU- and TUNEL-positive cells were counted on both sides of the entire cross-section of the palatal shelves at three rostro-caudal levels.

ChIP followed by qPCR was essentially performed as previously described (Voss et al., 2012) with the modifications described (Vanyai et al., 2015). Specifically, microdissected pharyngeal arch tissue was prepared and sheared under the conditions described previously (Vanyai et al., 2015). ChIP was performed using antibodies against H3K9ac (Cell Signaling, #9649S) and H3K14ac (Cell Signaling, #7627S) in conjunction with magnetic beads or using agarose bead-bound anti-V5 antibody (Sigma, A7345). Anti-H3K9ac and anti-H3K14ac data are displayed as ChIP signal enrichment over a precipitated internal control region (B2M); anti-V5 is displayed relative to input fraction. Gene loci were drawn based on the UCSC Genome Browser (genome.ucsc.edu).

Single-end ChIP sequencing reads were aligned to the mouse mm10 genome using Rsubread (Liao et al., 2013), and BAM files were generated and processed using SeqMonk Mapped Sequence Data Analyser Version 0.29.0 (Brabraham Institute, Cambridge, UK). Peaks were identified using MACS peak caller.

For RNA-sequencing and RT-qPCR, total RNA was isolated from the maxillary component of the 1st pharyngeal arch at E10.5 and E11.5, from the palatal shelves at E13.5 and E14.0, and from E10.5 1st and 2nd pharyngeal arches (RNAeasy, Qiagen) from embryos on a C57BL/6 background.

Uniquely indexed libraries were generated per sample with the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina), according to manufacturer's instructions. Indexed libraries were sequenced on the NextSeq platform (Illumina), generating 75 bp single end reads, yielding a minimum of ∼15 million total reads per sample.

Sequence reads were aligned to the mouse mm10 genome using Rsubread (Liao et al., 2013). Mapped reads were assigned to Entrez genes using featureCounts (Liao et al., 2014) and the in-built RefSeq annotation in Rsubread. Gene annotation was downloaded from the NCBI (ftp://ftp.ncbi.nlm.nih.gov). Differential expression analyses were undertaken using the edgeR (Robinson et al., 2010) and limma (Ritchie et al., 2015) software packages. Genes were filtered from the analysis if they failed to achieve at least 0.5 counts per million reads (CPM) in at least four samples. Genes without current annotation were also removed. Libraries were scale normalised using the trimmed mean of log expression ratios (TMM) method (Robinson and Oshlack, 2010). Counts were transformed to log₂ counts per million with associated precision weights using voom (Law et al., 2014). Differential expression was assessed using linear models and robust empirical Bayes moderated t-statistics (Phipson et al., 2016). P values were adjusted to control the false discovery rate (FDR) below 5% using the Benjamini and Hochberg method. To increase precision, the linear models incorporated corrections for a batch effect and two surrogate variables. The batch effect adjusted for two mouse colony locations. The surrogate variables were computed from a singular value decomposition of the linear model residuals. Expression signature analysis of the differential expression results was conducted using rotation gene set tests (ROAST; Wu et al., 2010). Gene ontology analyses were conducted using limma's goana function.

To make multidimensional scaling plots (MDS) and heatmaps, the expression of each gene was summarised as a log₂-CPM value using the cpm function of edgeR, with a prior count of 5, and batch effects were removed using the removeBatchEffects function in limma. The plots themselves were made using the plotMDS and heatmap.2 functions. Distances on the MDS plot represent leading log₂ fold-change, which is the root-mean-square average of the largest 500 log₂ fold-changes between each pair of RNA samples. Heatmap z-values were standardized to have mean 0 and standard deviation 1 for each gene, and hierarchical clustering was applied to the genes. RNA-seq read distribution over the exons of the Moz locus are displayed in Fig. S1 and show loss of reads over the deleted exons 3 to 7 in the Moz^–/– E10.5 pharyngeal arches.

The number of observations used was calculated to enable detection of a 25% difference between genotypes at alpha 0.05 with 80% power. Embryos and pups were examined blinded to genotype, i.e. before genotyping, and machine data generation was used to avoid bias where possible (RNA-seq, ChIP-seq, RT-qPCR, FACS). No animals were excluded. The number of observations and the statistical tests used are specified in the figure legends.

We thank the following researchers for mice and reagents: K. Rajewsky, J. L. R. Rubenstein and A. McMahon for making cre-deleter, Dlx5^+/– and Wnt1-cre mice, respectively, available to the scientific community; M. Busslinger and S. L. Nutt for the FLAG-V5-BIO triple tag; and L. Robb for the Pitx2 cDNA probe. We are grateful for excellent technical support from S. Eccles, N. Downer, C. Gatt and F. Dabrowski.