We obtained wnt1 ^cre mice from the Jackson laboratory. The Rosa26 ^{CAG-LSL-Dlx2−3xFlag} mouse was constructed by the Shanghai Model Organisms Center, Inc. (Shanghai, China). To generate wnt1 ^cre ; Rosa26 ^Dlx2/- mice, which could specifically overexpress Dlx2 in neural crest cells, we mated wnt1 ^cre mice with Rosa26 ^{CAG-LSL-Dlx2−3xFlag} mice. Wildtype C57BL/6J mice were purchased from Shanghai Jihui Laboratory Animal Care Co. Ltd. (Shanghai, China). All mice were maintained under SPF conditions at the Animal Center of the Ninth People’s Hospital affiliated with Shanghai Jiao Tong University School of Medicine. The day of the appearance of a vaginal plug was defined as E0.5 in all timed pregnancies. Embryos at the E10.5 and E12.5 stages (12:00 h of the day when the vaginal plug was detected was counted as E0.5) and P0 pups were collected for subsequent experiments. All animal experiments were approved by the Animal Care and Usage Committee of the Ninth People’s Hospital affiliated to Shanghai Jiao Tong University School of medicine.

Micro-CT was performed using a SkyScan 1176 (Bruker, Germany). Micro-CT images were acquired from P0 mice, with an x-ray source voltage of 45 kV and current of 550 µA. The data were collected at a resolution of 18 µm. Volume rendering in 3D was achieved using Mimics Medical 21.0 (Materialize). We evaluated micro-CT scans from three replicates per genotype. All landmarks were determined based on Mouse Development (Eds. J Rossant and P.L.Tam, 2002) and www.getahead.la.psu.edu.

All bones used in this study were manually segmented. Micro-CT scanning data were uploaded to Mimics as DICOM files. The background noise from these segmentations and bones outside the scope of this study were manually removed using Mimics’ editor tools. The remaining craniofacial bones were isolated and labeled using pre-scale thresholds that allowed only bones to be labeled. Reconstruction data were then rendered using Mimics’ 3D calculation tools and analysis tools were used for the measurements of isolated bones. The mean measurements of the maxillary bones were compared between the P0 wildtype and Dlx2 overexpression groups.

GraphPad Prism v.8 for Windows (GraphPad Software, La Jolla, CA, United States) was used for the statistical analysis. For all graphs, error bars represent standard deviations. Independent two-tailed Student’s t-tests were applied for comparisons between two groups. Differences were considered to be statistically significant at p-values < 0.05.

In E10.5 and E12.5, pregnant C57BL/6 females were euthanized using isoflurane and cervical dislocation. The embryos were removed from the uterus into cold PBS and transferred into a 6 cm Petri dish, using a disposable glass straw. For sequences library construction, complete maxillary processes were carefully dissected out from embryos using micro tweezers under a stereomicroscope.

Total RNA was extracted using Trizol from freshly dissected E10.5 maxillary process tissues. Three independent RNA samples were prepared for each genotype (WT and wnt1 ^cre ; Rosa26 ^Dlx2/- ). We used 2 μg total RNA as input material for the library preparations for each sample. Sequencing libraries were generated using the NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (#E7530L, NEB, United States) following the manufacturer’s recommendations. The libraries were sequenced on an Illumina HiSeq X ten platform and 150 bp paired-end reads were generated. Sequenced reads were mapped to the mm 10 genome using STAR aligner version 2.7.3a. Comparisons between the RNA-seq datasets were performed using the DESeq2 package in R. Enrichment analyses and visualization of functional profiles of differentially expressed genes (DEGs) were performed using the clusterProfiler package in R.

Fresh maxillary process tissues were conserved in the GEXSCOPE^® Tissue Preservation Solution (Singleron) until library preparation. The scRNA-Seq libraries were constructed in accordance with the Singleron GEXSCOPE™ protocol from the GEXSCOPE™ Single-Cell RNA Library Kit (Singleron Biotechnologies). Pools were sequenced on the Illumina HiSeq X to generate 150 bp paired-end reads. Unsupervised clustering of cell populations was performed using the tSNE and UMAP analysis from the Seurat R package.

After obtaining fresh cells from the maxillary processes of E12.5 wildtype mice, the CUT&Tag libraries were prepared using the Hyperactive In-Situ ChIP Library Prep Kit for Illumina (Vazyme Biotech, TD901) as previously described (Zuo et al., 2021). Approximately 50,000 cells were used per sample. The Anti-DLX2 antibody (Abcam, ab272902) was used as the primary antibody and goat anti-rabbit IgG (Vazyme, Ab206-10-AA) was used as the secondary antibody. All CUT&Tag libraries were sequenced on the Illumina Nova Seq 6000 platform at PE150 mode. Low-quality reads and adapters were trimmed by Trim Galore (v0.6.5). The clean reads were mapped to the mm 10 genome using bowtie2 (v2.4.2).

The wnt1 ^cre ; Rosa26 ^Dlx2/- mice, with Dlx2 overexpressed in neural crest-derived cells exhibit craniofacial deformities such as cleft palate (Sun et al., 2022). To quantify the effect of Dlx2 overexpression on craniofacial bone formation, we performed Micro-CT scanning on the head of P0 wild-type and wnt1 ^cre ; Rosa26 ^Dlx2/- mice. Micro-CT analysis provides comprehensive information on anatomical landmarks and the size of each craniofacial bone (Ho et al., 2015). The premaxilla and nasal bone, maxilla, palatine bone, frontal bone, parietal bone, interparietal bone, occipital bone and mandible were isolated for analysis (Figure 1A). Obvious ectopic bone formation and absorption were found in the premaxilla and nasal bone, frontal bone and parietal bone of the wnt1 ^cre ; Rosa26 ^Dlx2/- mice. The cranial anteroposterior diameter of the wnt1 ^cre ; Rosa26 ^Dlx2/- mice was significantly smaller than that of wild-type mice. We isolated the maxilla from wild-type and wnt1 ^cre ; Rosa26 ^Dlx2/- mice and defined the anatomical landmarks (Figure 1B). We next quantitatively compared the sizes of the maxilla using the landmarks (Figures 1C–G). The full width and half-width of the maxilla of the wnt1 ^cre ; Rosa26 ^Dlx2/- mice were significantly decreased (Figures 1D, E) but there was no significant difference in the length and height (Figures 1C–F). The distance between the left and right halves of the maxilla in wnt1 ^cre ; Rosa26 ^Dlx2/- mice was significantly increased (Figure 1G), which was consistent with the cleft palate phenotype. In summary, the quantitative comparison of Micro-CT scans revealed that Dlx2 overexpression had a teratogenic effect on the mouse maxilla.

In previous research, it was found that the overexpression of Dlx2 had an impact on gene expression in E12.5 maxillary processes. However, the temporal and spatial expression analysis of Dlx2 showed that the overexpression was already evident in the earliest stage (E10.5) of maxillary process formation (Sun et al., 2022). In order to further understand how the overexpression of Dlx2 affects the development of maxillary processes, bulk RNA-Seq was performed on the maxillary processes of E10.5 wild-type and wnt1 ^cre ; Rosa26 ^Dlx2/- mice. The individual replicates exhibited a high degree of correlation among the same genotype but a lower correlation was observed between replicates of different genotypes (Figure 2A), suggesting the overexpression of Dlx2 already induced transcriptome changes at this stage.

Comparisons between the wnt1 ^cre ; Rosa26 ^Dlx2/- and wild-type mice revealed that 6,230 genes exhibited significant expression changes. Of these genes, 2,192 genes were significantly upregulated and 1,762 genes were significantly downregulated (Figure 2B). The Gene Ontology (GO) enrichment analysis of DEGs that were significantly upregulated revealed that they are involved in critical biological processes and molecular pathways, such as organic cyclic compound catabolic process, peptidyl-lysine modification, ncRNA processing and RNA catabolic process. The upregulated genes were also involved in a variety of RNA metabolic processes, which included ncRNA metabolic process, regulation of mRNA metabolic process, tRNA metabolic process and mitochondrial RNA metabolic process (Figure 2C). Notably, the downregulated DEGs are enriched for functional terms related to neuronal development, such as synapse organization, axonogenesis, dendrite development and axon guidance (Figure 2D), which reflects the neural crest-origin of the maxillary processes.

We found the DEGs of E10.5 are significantly different from the previously reported bulk RNA-Seq DEGs of the mouse maxillary process in E12.5 wildtype and wnt1 ^cre ; Rosa26 ^Dlx2/- mice (Sun et al., 2022). Between the 6230 E10.5 DEGs and the 2428 E12.5 DEGs, only 780 genes are common (Figure 2E). Among the 5,450 genes that were unique to E10.5, the most enriched GO terms are related to cell proliferation and transcription, such as mRNA processing, DNA repair, regulation of DNA metadata process and mitotic nuclear division (Figure 2F). The 780 DEGs that were shared by E10.5 and E12.5 are enriched for GO terms involved in neuronal development, such as axonogenesis, regulation of neurogenesis, neuron project guidance and axon guidance (Figure 2G). Notably, the 1648 DEGs that were unique to E12.5 mice are enriched for ossification related genes (Figure 2H). Thus, the overexpression of Dlx2 can lead to different transcriptional responses and physiological outcomes at different stages of craniofacial development.

Our bulk RNA-Seq analyses reveal pronounced transcriptional regulatory effects of Dlx2 during early maxillary development. However, the inability to distinguish among different cell subpopulations within maxillary processes precludes further dissection of transcriptome changes associated with the differentiation of mesenchyme. Overcoming these limitations requires transcriptome profiling at single-cell resolution.

In the maxillary process at E10.5, the differentiation of most tissue types has not occurred and the mesenchymal cell population is relatively homogeneous. In order to more clearly reveal whether overexpression of Dlx2 affects the differentiation trajectory of cells, we isolated the maxillary process tissues of E12.5 wild-type and wnt1 ^cre ; Rosa26 ^Dlx2/- mice for single-cell RNA sequencing. The two scRNA-Seq datasets were combined and further analyzed. After dimensional reduction, the main cell populations from the two different samples largely overlapped with each other, indicating the batch effect was minimal (Figure 3A). The combined data were further clustered into 14 cell populations, with the largest cell populations corresponding to mesenchymal cells (Figure 3B). The nine mesenchymal clusters (clusters 0, 1, 2, 3, 4, 6, 7, 8, and 9) were selected for further analysis (Figure 3C).

In order to identify the cell populations most affected by Dlx2 overexpression, we quantified the relative proportions of wnt1 ^cre ; Rosa26 ^Dlx2/- mice samples to wild-type cells in each of the 9 mesenchymal cell clusters. In each cell cluster, the number of cells from wnt1 ^cre ; Rosa26 ^Dlx2/- mice samples was divided by the number of cells from the wild-type samples. We found that the wnt1 ^cre ; Rosa26 ^Dlx2/- cells were relatively enriched in clusters 0, 3, and 6, while depleted in the other clusters, particularly for clusters 1 and 8 (Figure 3D). To further understand the identities of clusters 0, 3, and 6, we examined their marker genes. For each of these clusters, several marker genes related to different tissue types can be found (Figure 3E), suggesting these cells may represent various precursor cells that have not fully committed to a specific cell type.

Notably, we also found that there are a large number of marker genes related to the cell cycle in each cell population (Supplementary Figure S1). We assessed the cell cycle stages for each mesenchymal cell and found that clusters 0, 3, and 6 mainly consisted of cells in the G1 phase, while clusters 1 and 8 consisted of cells in the G2M phase (Figure 3F). When comparing the cell cycle composition of mesenchymal cells from the two genotypes, we found that the proportion of cells in the G1 phase markedly increased, while the proportion of cells in G2M and S phase decreased in the wnt1 ^cre ; Rosa26 ^Dlx2/- mice cells (Figure 3G). These results suggest the overexpression of Dlx2 led to a slowdown of cell cycle progression and inhibition of cell proliferation.

To further understand how the overexpression of Dlx2 affects the differentiation trajectory of maxillary mesenchymal cells, we performed pseudotime developmental trajectory analysis on the combined scRNA-Seq dataset. While the cells from wild-type and Dlx2-overexpressing mice exhibit similar trajectories without obvious divergence, the cells from the wnt1 ^cre ; Rosa26 ^Dlx2/- mice were located at more downstream positions on the pseudotime trajectory compared to the wild-type cells (Figure 3H). This difference was further confirmed by quantifying the pseudotime scores for the cells from the two genotypes (Supplementary Figure S2). These analyses thus suggest the overexpression of Dlx2 caused the mesenchymal cells within the maxillary process to enter a more differentiated state.

Taken together, our scRNA-Seq result suggests that overexpression of Dlx2 had two related effects: inhibition of cell proliferation and promotion of differentiation. In Dlx2-overexpressing mice, the maxillary process cells may have precociously entered a more downstream differentiation state before they had sufficient proliferation, thereby impairing the development of the maxillary bone and ultimately causing the phenotypes of narrowing width, widening distance and cleft palate.

To uncover the regulatory mechanism of Dlx2 in early maxillary process development, CUT&Tag analysis was performed. CUT&Tag is a novel and highly sensitive method used to identify transcription factor occupancy sites (Kaya-Okur et al., 2019; Kaya-Okur et al., 2020). We used this method to identify candidate direct targets of DLX2. We performed DLX2 CUT&Tag on two replicates of wild-type mice maxillary processes and identified 14,738 and 6899 peaks in each replicate. Intersection was used to obtain 3518 peaks that were common to both replicates. Through annotation of these peaks, we found that less than 7% were located in the promoter area (within 2 kb from the TSS) (Figure 4A). The largest proportion of DLX2 occupancy sites was located between genes, which indicated that DLX2 may bind to potential enhancer regions to regulate the expression of protein-coding genes (Figure 4A).

The ten most enriched known motifs identified by HOMER software are listed in Figure 4B. Among these enriched motifs, Runx2 was of particular interest because substantial in vivo and in vitro studies have shown that this gene is strongly associated with osteogenesis (Tosa et al., 2019; Deiana et al., 2020). This suggested that Dlx2 may function in collaboration with Runx2 to reshape the transcriptome when Dlx2 is overexpressed. Mnt is likely to be a transcriptional repressor and an antagonist of Myc-dependent transcriptional activation and cell growth (Hurlin et al., 1997), which may explain in part the inhibition of cell proliferation found by scRNA-Seq.

The CUT&Tag peak annotation identified 2,511 genes that are associated with DLX2 peaks. By cross-referencing these genes with bulk RNA-Seq upregulated DEGs, 83 upregulated Dlx2 target genes were obtained (Figure 4C). These genes were enriched for genes involved in the regulation of binding, regionalization, negative regulation of cell development and regulation of Notch signaling pathway (Figure 4E). Interestingly, the Notch signaling pathway has been shown to play an important role in palatal development (Casey et al., 2006). The 254 downregulated Dlx2 target genes were enriched for genes involved in axonogenesis and synapse organization (Figures 4D, F), consistent with the earlier analysis results. Among these downregulated genes, key osteogenic genes such as Zeb2 (Wang et al., 2022) and Rora (Tao C et al., 2022) were significantly expressed in mesenchymal cell clusters 0 and 3 of scRNA-seq respectively, and these two clusters of cells constituted the majority of the mesenchymal cell group of wnt1 ^cre ; Rosa26 ^Dlx2/- . Overexpression of Dlx2 affects the osteogenesis of most mesenchymal cells. These putative Dlx2 target genes may be the most direct effectors in the downstream regulatory network of Dlx2.