Runx2-Twist1 interaction coordinates cranial neural crest guidance of soft palate myogenesis

Cranial neural crest (CNC) cells give rise to bone, cartilage, tendons, and ligaments of the vertebrate craniofacial musculoskeletal complex, as well as regulate mesoderm-derived craniofacial muscle development through cell-cell interactions. Using the mouse soft palate as a model, we performed an unbiased single-cell RNA-seq analysis to investigate the heterogeneity and lineage commitment of CNC derivatives during craniofacial muscle development. We show that Runx2, a known osteogenic regulator, is expressed in the CNC-derived perimysial and progenitor populations. Loss of Runx2 in CNC-derivatives results in reduced expression of perimysial markers (Aldh1a2 and Hic1) as well as soft palate muscle defects in Osr2-Cre;Runx2fl/fl mice. We further reveal that Runx2 maintains perimysial marker expression through suppressing Twist1, and that myogenesis is restored in Osr2-Cre;Runx2fl/fl;Twist1fl/+ mice. Collectively, our findings highlight the roles of Runx2, Twist1, and their interaction in regulating the fate of CNC-derived cells as they guide craniofacial muscle development through cell-cell interactions.


Introduction
The craniofacial musculoskeletal complex is an important evolutionary innovation in vertebrates that facilitates feeding, breathing, facial expression, and verbal communication. One unique component of this complex is the cranial neural crest (CNC) cells. CNC cells and their derivatives give rise to all facial bones, ligaments, and muscle connective tissues including tendons and fascia that directly surround muscle cells (Chai and Maxson, 2006;Heude et al., 2010;Le Douarin et al., 2004). Recently, CNC cells have been shown to regulate formation of mesoderm-derived craniofacial muscles through cell-cell interactions. Mouse genetic studies have further shown that CNC cells and their derivatives surround myogenic cells, facilitate myogenic cell migration, and establish cellular scaffolding at future myogenic sites to regulate muscle morphogenesis (Grimaldi et al., 2015;Han et al., 2014;Rinon et al., 2007). For instance, disruption of Dlx5/6, which is specifically expressed by CNC-derived cells in the mouse, leads to the loss of all first pharyngeal arch-derived masticatory muscles and second pharyngeal arch-derived muscles (Heude et al., 2010). Proliferation and survival of CNC-derived cells and fourth to sixth pharyngeal arch-derived myogenic cells in the soft palate are also affected, resulting in a truncated soft palate in Dlx5 -/mice (Sugii et al., 2017). Similarly, TGFb signaling in CNC-derived cells is critical for proliferation and differentiation of tongue and masseter muscle cells (Han et al., 2014;Hosokawa et al., 2010;Iwata et al., 2013). It is important to note that the transcription factors and signaling pathways critical for the role of CNC-derived cells in myogenesis are not restricted in their expression to merely the CNC-derived cells surrounding the muscle, known as perimysial cells; they are also expressed in other CNC-derived musculoskeletal tissues (e.g. bones, bone eminences, and tendons) and regulate their development (Depew et al., 2002;Hosokawa et al., 2010;Zhao et al., 2008). This suggests that the same transcription factors and signaling pathways could activate cell-type-specific responses in multiple components of the musculoskeletal complex that may help coordinate the development of this intricate system. Therefore, it is important to investigate the cell-type-specific signaling mechanisms that regulate the heterogeneous CNC-derived cells and reveal their impact on craniofacial musculoskeletal development.
The soft palate is a muscular structure that comprises the posterior third of the palate. Its movement opens and closes the nasopharynx and oral cavity to direct air and food into different passages, as well as during speech. Several components of the soft palate are CNC-derived, including perimysial cells, palatal stromal cells that constitute the majority of palatal shelf mesenchyme, and tendons. In contrast, the soft palatal muscles are derived from pharyngeal mesoderm (Grimaldi et al., 2015). Five muscles are involved in the human soft palate. They include the tensor veli palatini (TVP) and levator veli palatini (LVP), which descend from the skull base and elevate the soft palate, and the palatoglossus (PLG) and palatopharyngeus (PLP), which ascend from the tongue and the pharyngeal wall, respectively, and depress the soft palate (Li et al., 2019). The fifth muscle, the musculus uvulae, which is specific to humans, is located at the end of the soft palate. Patients with cleft palate often have multiple types of tissue abnormalities including bone defects and insufficient, misoriented muscle fibers (Dixon et al., 2011;Li et al., 2019). Functional restoration of cleft soft palate is challenging because the muscles have limited ability to regenerate after surgical repair of the cleft (Von den Hoff et al., 2019). Therefore, comprehensive understanding of the growth and transcription factors that regulate the coordinated development of the distinct tissues in the soft palate is of both scientific and clinical significance.
Runx2, a known regulator of skeletogenesis and odontogenesis, is a Runt DNA-binding domain family transcription factor and contains multiple activation and repression domains. Patients with haploinsufficiency of RUNX2 exhibit cleidocranial dysplasia, which is associated with specific skeletal and dental phenotypes. During osteoblast differentiation, Runx2 acts as a master organizer, recruiting phosphorylated Smad1/5, c-Fos, and c-Jun to activate expression of osteoblast-specific collagen and fibronectin upon receiving BMP signals and parathyroid hormones; it also binds histone deacetylases to repress cell cycle inhibitors and stimulate proliferation (Schroeder et al., 2005). Despite its well-known roles in regulating hard tissue development, the importance of Runx2 in soft tissue development has not been studied. Interestingly, several clinical case reports reveal that some RUNX2-deficient patients have thin masseter muscles, cleft lip, or high-arched palate (Furuuchi et al., 2005;Sapp et al., 2004;Sull et al., 2008;Yamachika et al., 2001). These studies hint that Runx2 may regulate the development of the palatal muscles and other components in sync with the bone to form the intricate craniofacial musculoskeletal complex by performing multiple tissue-specific roles.
In this study, we performed an unbiased transcriptional profile analysis of the developing soft palate using single-cell RNA-seq (scRNA-seq). We identified cellular-level heterogeneity in the CNCderived soft palate mesenchyme, associated with distinctive cell fates: perimysial and midline mesenchymal lineages, as well as previously unknown cell types associated with putative progenitors. In addition, we found Runx2 was expressed in non-osteochondrogenic cells in the perimysial populations and in CNC-derived progenitor cells. Consistent with its expression pattern, loss of Runx2 in CNC-derived cells resulted in a soft palate cleft along with tendon, bone, and muscle differentiation defects. We further revealed that loss of Runx2 led to ectopic expression of Twist1 and reduction in the expression of perimysial marker genes (Aldh1a2 and Hic1) in CNC-derived perimysial cells. We also identified that suppression of Twist1 expression by Runx2 is important for the development of palatal muscles and for maintaining the expression of the perimysial marker and myogenic-promoting gene Aldh1a2, thus coordinating soft palate morphogenesis by orchestrating the fate determination of CNC-derived mesenchymal lineages. Taken together, our findings reveal that Runx2 regulates distinct downstream targets in different subgroups of CNC-derived cells to fine-tune the development of craniofacial structures.
To characterize the roles of these less known subpopulations, we analyzed the top 10 differentially expressed genes in each cluster and performed functional annotation for those highly specific markers using Ingenuity Pathway Analysis. We thus identified four major types of CNC-derived cells in the soft palate besides osteogenic (Cluster 2) and chondrogenic cells (Cluster 10) ( Figure 1B). Cluster 0 was enriched with early CNC marker genes such as Tfap2b, Six2, and Prrx1 (Simõ es-Costa and Bronner, 2015;Soldatov et al., 2019), so we suspected that this population might be an undifferentiated early progenitor population associated with early CNC cells, and accordingly we hypothesized that they were CNC-derived progenitors ( Figure 1B). Genes enriched in Cluster 1 (Tbx22, Wnt16, Meis2) were associated with the palatal shelf midline during development (Louw et al., 2015;Pauws et al., 2013;Warner et al., 2009; Figure 1B; Figure 1-figure supplement 1B); hence, we refer to this cluster as midline mesenchymal cells. Clusters 3, 4, and 7 expressed high levels of genes related to head and muscle morphogenesis (Cxcl12, Igf1, Aldh1a2); we refer to them as perimysial cells (Matt et al., 2008;Schiaffino and Mammucari, 2011;Vasyutina et al., 2005; Figure 1B; Figure 1-figure supplement 1C). Interestingly, Cluster 8 was strongly enriched in genes associated with mitosis (Top2a, Ccnb1, Ube2c) (Nielsen et al., 2020;Pines, 2011;Strauss et al., 2018) even after cell cycle regression was performed ( Figure 1B). We therefore refer to this cluster as mitotic cells.
To investigate the in vivo identities of each cluster, we performed RNAscope in situ hybridization of the soft palate at E14.5. Different soft palate myogenic sites develop sequentially from anterior to posterior direction. Specifically, in coronal sections, the unfused palatal shelves in the LVP region (posterior) protrude toward the midline and the myogenic cells grow in a lateral to medial direction along the palatal shelves at E14.5 ( Figure 1C), while the palatal shelves in the TVP region (anterior) are already fused and the myogenic cells wrap around the pterygoid plate ( Figure 1D). As the TVP and LVP myogenic sites are more identifiable than those of the PLG and PLP at E14.5, we used the former as reference points for the anatomical locations of each cluster in our analysis. Using the top enriched genes of each cluster, we identified their distinct anatomical locations in vivo. The Ibsp+ osteogenic ( Figure 1E,M) and Col2a1+ chondrogenic clusters ( Figure 1F,N) were mostly associated with part of the styloid process of the temporal bone in the LVP region and the pterygoid plate of the sphenoid bone in the TVP region. In the LVP region, the Tfap2b+ progenitor cluster was mainly located in the lateral portions of the palatal shelves ( Figure 1G). The majority of the Aldh1a2+ perimysial cluster was distributed in the lateral portion while only a small portion of this cluster appeared in the central myogenic sites ( Figure 1H). In contrast, the two other perimysial clusters (Hic1+ and Tbx15+) were most abundantly located in the central myogenic sites of the LVP ( Figure 1I-J). Midline mesenchymal cells (Tbx22+) were mainly located in the medial portions of the palatal shelves ( Figure 1K). The Top2a+ mitotic cells were distributed throughout the palatal shelves and adjacent to both early progenitors and committed CNC-derived cells ( Figure 1L). A similar distribution of Because the oropharyngeal muscles are only present in the soft palate, not the hard palate, we investigated whether the perimysial markers are also specific to the soft palate. Interestingly, Tbx15 expression was absent from the hard palate, but Aldh1a2 and Hic1 were expressed in the hard palate mesenchyme specifically surrounding the tooth germ ( Figure 1-figure supplement 3D-F). This suggests that Aldh1a2 and Hic1 might have different functions in the hard and soft palate. Interestingly, we also observed that Aldh1a2 was expressed in the medial mesenchyme of the tongue, while Hic1 and Tbx15 were expressed broadly in the mesenchyme of the tongue at E14.5 (Figure 1 As myogenic precursors started to appear in the center of tongue primordium (Han et al., 2012), the Aldh1a2+ population might be specifically associated with early myogenic populations, while Hic1+ and Tbx15+ populations may be associated with more general myogenic populations.

Runx2 is expressed in the perimysial populations and CNC-derived progenitor cells during soft palate development
To elucidate the dynamic process by which CNC-derived cells differentiate during soft palate development, we performed individual single-cell transcriptome analyses for E13.5, E14.5, and E15.5, then compared them. The pterygoid plate of the sphenoid bone and part of the styloid process of the temporal bone are not considered to be part of the palate, so we excluded the osteochondrogenic clusters belonging to these structures from further analysis. Interestingly, we observed decreased cell heterogeneity in CNC-derived soft palate mesenchymal populations during development. The number of CNC-derived clusters declined from seven at E13.5 to six at E14.5 and eventually five at E15.5 using the same unsupervised clustering settings (Figure 2-figure supplement 1A). In contrast, myogenic cells formed a single cluster at E13.5 but expanded to two clusters at E15.5 ( Figure 2-figure supplement 1A).
To further investigate how each cluster changed over time, we extracted and compared the CNC-derived and myogenic cells from E13.5 to E15.5 based on the earlier integration analysis ( Figure 2A). Consistent with previous observations, in the myogenic clusters we observed an increased number of both early myogenic precursors (Cluster 9, Msc+, Myf5+) and differentiated myocytes (Cluster 17, Myl4+) as development progressed (Figure 2A). The number of cells in Clusters 1 (Tbx22+), 3 (Hic1+) and 4 (Tbx15+) also increased from E13.5 to E15.5, but the number of cells in CNC-derived Clusters 0 (Tfap2b+), 7 (Aldh1a2+), and 8 (Top2a+) gradually decreased. This suggests that Clusters 0, 7, and 8 may be progenitors that are transiently present at early stages of soft palate development and give rise to Clusters 1, 3, and 4 ( Figure 2A-B). To test this, we also  computationally predicted the differentiation trajectory of CNC-derived cells using pseudotime analysis. Our results predicted Cluster 0 to be common CNC-derived progenitors that bifurcate into two more committed groups: perimysial progenitors (Cluster 7) for the later perimysial population (perimysial fibroblasts) (Clusters 3 and 4), and another group of progenitors (a subset of Cluster 0 and Cluster 8) for midline mesenchymal cells (Cluster 1) (Figure 2-figure supplement 1B-C). The integration analysis suggested that the fate decision between perimysial and midline mesenchymal cells happens at E13.5-E14.5. Cluster 8, predicted to be a more committed group of progenitors, represents the Top2a+ mitotic population. Because mitosis establishes a time window during which transcription factors can easily access and activate genes important for cell lineage determination (Gurdon, 2016), cells with high mitotic activity are likely undergoing cell fate transition. Therefore, those Top2a+ mitotic cells might be transitioning from early progenitor status to becoming more committed to a particular fate. Interestingly, the number of Aldh1a2+ cells in Cluster 7 gradually decreased, but the number of Aldh1a2+ cells increased in Clusters 3 and 4 from E13.5 to E15.5 ( Figure 2B), probably because Aldh1a2 labels both the majority of the early perimysial population (Cluster 7) and also some of the late perimysial populations (Clusters 3 and 4).
Notably, we observed that expression of Runx2 in the CNC-derived mesenchyme gradually decreased from E13.5 to E15.5, suggesting it might play a role in regulating CNC-derived cell differentiation during early soft palate development ( Figure 2B). Furthermore, Runx2 was expressed not only by the CNC-derived common progenitors Cluster 0 (Tfap2b+), but also by other more committed progenitor cells, Cluster 8 (Top2a+) and Cluster 7 (Aldh1a2+), which were mainly distributed around the bifurcation regions of different lineages; only a few midline mesenchymal cells in Cluster 1 (Tbx22+) expressed Runx2 ( Figure 2C-J).
To investigate the functional significance of Runx2 for soft palate development in vivo, we examined Runx2 expression in the TVP and LVP regions of control mice. Double staining of Runx2 and the myogenic marker MyoD/Myod1 from E13.5-E15.5 revealed changes in Runx2 expression in the myogenic region as development progressed. Runx2 expression was gradually restricted from most of the palate primordium at E13.5 to only the mesenchymal cells in the putative progenitor, perimysial, and osteogenic sites in the LVP region at E14.5; eventually, it was found only in the osteogenic regions at E15.5 ( Figure 2K-M). As there was no detectable Runx2 expression in the TVP perimysial site from E13.5 to E15.5 ( Figure 2N-P), we focused on the LVP region as we investigated the colocalization of Runx2 with markers of early CNC-derived progenitors and different lineages in vivo. Consistent with our single-cell analysis, Runx2 was predominantly expressed in the putative progenitor population (Tfap2b+), actively amplifying population (Top2a+) and perimysial cells (Aldh1a2+), with only a few in the midline mesenchymal cells (Tbx22+) ( Figure 2Q-X). Previous studies have shown that CNC-derived cells guide the migration and potentially regulate the maturation of mesodermderived myogenic precursors in the soft palate through tissue-tissue interactions (Grimaldi et al., 2015;Li et al., 2019;Sugii et al., 2017). We hypothesized that Runx2 may regulate differentiation of CNC-derived cells in a cell-autonomous manner at early stages, which may indirectly affect myogenesis in the soft palate mesenchyme.

Loss of Runx2 in CNC-derived cells results in soft palate development defects
To test the functional significance of Runx2 in regulating soft palate muscle development, we specifically targeted Runx2 in CNC-derived palate mesenchymal cells. We first tested whether Osr2-Cre, Figure 2 continued palate development. (C-J) Co-expression of Runx2 with cluster-specific markers Tfap2b, Top2a, Aldh1a2, Tbx22 in E13.5-E15.5 soft palate integration analysis. Boxed areas in (C-F) are enlarged in (G-J). Black arrows point to cells co-expressing Runx2 with individual cluster-specific markers. (K-P) Runx2 with myogenic markers MyoD or Myod1 on coronal sections of the tensor veli palatini (TVP) and levator veli palatini (LVP) regions of control mice at E13.5, E14.5, and E15.5. Boxes indicate regions shown at higher magnification in the insets. (Q-X) Co-localization of Runx2 with cluster-specific marker genes Tfap2b, Aldh1a2, Top2a, and Tbx22 on coronal sections of the LVP region of E14.5 control mice. Boxed areas in Q-T are enlarged in U-X. Yellow dashed lines in (K-T) outline the myogenic cells. White dashed lines outline the palatal shelf. Scale bars in K and Q indicate 100 mm for K-P and Q-T. Scale bar in U indicates 30 mm for U-X. The online version of this article includes the following figure supplement(s) for figure 2: which specifically labels the CNC-derived cell subset in the developing palatal mesenchyme from the beginning of palatal shelf outgrowth (Lan et al., 2007), could also label the CNC-derived population in the soft palate. We confirmed in Osr2-Cre;tdTomato mice that tdTomato+ cells indeed contribute to soft palate mesenchyme including the perimysial cells surrounding all soft palatal muscles as early as E14.5 (Figure 3-figure supplement 1A-G). Furthermore, co-expression of tdTomato and Runx2 in the soft palate suggested that we could use Osr2-Cre to specifically delete Runx2 in a subset of CNC-derived cells in the soft palate region (Figure 3-figure supplement 1A-G).
To test whether Runx2 is a key regulator of soft palate development, we next generated Osr2-Cre;Runx2 fl/fl mice, which showed cleft soft palate (5/10), misoriented muscle fibers and reduced muscle size (10/10) along with defects in hard tissues including the palatine bone (3/6) and pterygoid process ( . Notably, in analyzing the CT scans, we found that three out of six Runx2 mutant mice with missing palatine bones and more severe pterygoid plate defects also had soft palate clefts, while the other three Runx2 mutants without clefts exhibited palatine bones that were smaller, though not statistically significantly so, and less severe pterygoid plate defects, particularly shorter pterygoid plate height (Figure 3-figure supplement 2A-J), suggesting the severity of skeletal defects is associated with the variability of soft palate clefts in Osr2-Cre;Runx2 fl/fl mice. Consistent with the CT scans, histological analysis showed that the height of pterygoid plate was reduced and muscle attachment was abnormal in the TVP region of Osr2-Cre;Runx2 fl/fl mice ( Figure 3G-H,L-M). Because the aponeurosis serves the important function of attaching the hard tissue to the muscle, we also analyzed the fibrous tendon tissue marked by Scx by RNAscope in situ hybridization. The tendon tissue did not extend to the midline in the palate primordium in the TVP region in Osr2-Cre;Runx2 fl/fl mice as it did in the controls at E16.5 (Figure 3-figure supplement 3I-L). It could be seen more clearly at P0 that the aponeurosis in Osr2-Cre;Runx2 fl/fl mice was thinner and it did not stretch from the lateraloral side to the medial-nasal side as it did in control mice (Figure 3-figure supplement 4A-D), suggesting its attachment to the posterior bone probably was likely abnormal. As Runx2 is not expressed in the perimysial site of the TVP region, Therefore, this muscle attachment defect of the TVP might be due to disruption of the hard tissue and aponeurosis. In the LVP region, the muscles were reduced in size in Osr2-Cre;Runx2 fl/fl mice compared to controls ( Figure 3I-K,N-P). Interestingly, a significant number of LVP muscle fibers had anterior-posterior alignment in Osr2-Cre; Runx2 fl/fl mice ( Figure 3O-P), compared to the uniform lateral-medial alignment of LVP muscle fibers in controls ( Figure 3J-K). This suggests that the muscle fibers were mis-oriented, similar to the phenotype seen in patients with cleft soft palate. Additionally, we observed that the muscle fibers had centralized nuclei in the soft palate of Osr2-Cre;Runx2 fl/fl mice, which suggests that they had muscle differentiation defects ( Figure 3P). Similar muscle defects were also observed in other palatal muscles such as the PLP (Figure 3-figure supplement 3C-D,G-H). We concluded that loss of Runx2 leads directly to defects in CNC-derived cells and indirectly to muscle defects in the soft palate.
To examine soft palate muscle differentiation in Osr2-Cre;Runx2 fl/fl mice, we analyzed expression of myogenic markers at multiple developmental stages to identify the time point at which muscle defects began to appear using the LVP as an example. In the LVP, there was no apparent change of early myogenic marker MyoD expression between control and Osr2-Cre;Runx2 fl/fl mice at E13.5 (Figure 4-figure supplement 1A-B). MyoD staining revealed that defects started to appear at E14.5 ( Figure 4A-B), when the palatal shelves began to grow and protrude towards the midline. Expression of the late myogenic marker MHC was decreased in the soft palate of Osr2-Cre;Runx2 fl/fl mice compared to controls at E15.5 ( Figure 4C-D), suggesting delayed muscle differentiation. This reduced expression of MHC persisted in the soft palate of Osr2-Cre;Runx2 fl/fl mice at E16.5 (Figure 4-figure supplement 1C-D). MHC staining suggested that the myoblasts had fused to form myofibers, which were uniformly aligned in layers running in the lateral-to-medial direction in the LVP of control samples at E15.5 ( Figure 4C). However, more immature myoblasts and fewer differentiated myofibers were present in Osr2-Cre;Runx2 fl/fl mice ( Figure 4D) and the myofibers extended in different directions, potentially hindering further muscle development and compromising physiological function.
To investigate the cellular mechanism underlying soft palate defects in Osr2-Cre;Runx2 fl/fl mice, we analyzed cell proliferation, apoptosis, and differentiation. Consistent with the MyoD expression pattern, we did not detect any change in the number of BrdU+ proliferating cells in the LVP region of Osr2-Cre;Runx2 fl/fl mice compared to controls at E13.5 (Figure 4-figure supplement 1A-B). In the LVP region at E14.5 and E15.5, the proliferation rate of MyoD-CNC-derived cells did not have significant difference in the perimysial sites of Osr2-Cre;Runx2 fl/fl mice compared to controls ( Figure 4E-H,I,K). The proliferation rate of MyoD+ myogenic cells was not significantly different between controls and Runx2 mutants at E14.5 ( Figure 4E-F,J), but a significant reduction in the proliferation rate was observed in Runx2 mutants at E15.5 ( Figure 4G-H,L). We also performed cas-pase3 immunofluorescence staining to investigate cell apoptosis. The number of apoptotic cells was       indistinguishable between controls and mutants at E14.5 and E15.5 (Figure 4-figure supplement  1E-H). It is worth noting that although the proliferation rate of CNC-derived cells in the Runx2 mutant mice were not significantly different from that of the controls, we observed that they had fewer and less proliferative MyoD+ cells than Runx2 fl/fl control mice. These differences might be due to altered signaling in CNC-derived cells causing the reduction of MyoD expression as well as proliferation defects of myogenic cells in Osr2-Cre;Runx2 fl/fl mice.

Runx2 plays an important role in the lineage commitment of CNCderived cells in the soft palate
To investigate whether Runx2 regulates CNC-derived cell fate determination during soft palate development, we compared cell composition and gene expression profiles of E14.5 Osr2-Cre; Runx2 fl/fl and control soft palates using scRNA-seq, bulk RNA-seq, and in vivo expression analyses. Using integration analysis based on shared variance, we identified similar cell clusters in the soft palates of control and Osr2-Cre;Runx2 fl/fl mice at this stage. However, the composition of the CNCderived cells was altered in the Runx2 mutants compared to controls ( Figure 5A). Using markers of different subtypes of CNC-derived cells, we observed that the percentage of perimysial cells (Cluster 4) in the population decreased in Runx2 mutants, while the percentage of midline mesenchymal cells (Clusters 0 and 3) increased ( Figure 5B). Moreover, in situ RNAscope staining revealed decreased expression of perimysial markers in the soft palates of Runx2 mutants compared to controls at E14.5, suggesting the CNC-derived perimysial populations were affected, potentially leading to further myogenic defects ( Figure 5-figure supplement 1A-H). Consistent with the scRNA-seq results, bulk RNA-seq also identified that certain genes associated with specific types of CNC-derived cells were differentially expressed in the Osr2-Cre;Runx2 fl/fl mice ( Figure 5-figure supplement 1I). A number of genes not exclusively associated with specific types of CNC-derived cells, including Twist1 and Meox2, were also identified as being differentially expressed in the bulk RNA-seq analysis ( Figure 5-figure supplement 1I).
We focused our attention on Twist1, which inhibits binding of Runx2 to its downstream targets and antagonizes Runx2's function in osteoblasts (Bialek et al., 2004). We began by analyzing the expression pattern of Twist1 during soft palate development. Based on the integration analysis of E13.5-E15.5 single-cell transcriptomes from controls, we observed Twist1 was primarily expressed in midline mesenchymal cells (Tbx22+), while its expression in CNC-derived common progenitors and perimysial cells (Tfap2b+; Aldh1a2+) was relatively low ( Figure 5-figure supplement 2A). In addition, expression of Twist1 in the CNC-derived cells changed over time. In the LVP region, at E13.5 Twist1 was expressed at a low level in the palate primordium and perimysial sites (Figure 5-figure supplement 2B-C). At E14.5, expression of Twist1 in the palate primordium had increased, whereas its expression was maintained at a low level in the perimysial site ( Figure 5-figure supplement 2D-E). At E15.5, Twist1 showed a similar expression pattern to that of E14.5 ( Figure 5-figure supplement 2F-G). This spatiotemporally specific Twist1 expression in the palate primordium and myogenic regions of the soft palate was accompanied by an opposite trend in Runx2 expression in the same regions at the same stages. This is perhaps shown most clearly by the colocalization of Runx2 and Twist1 in the LVP region at E14.5 ( Figure 5D-G), which suggested that expression levels of Runx2 and Twist1 are tightly coordinated during soft palate development. Interestingly, we discovered that expression of Twist1 was upregulated in most of the palatal shelf region including the  Haploinsufficiency of Twist1 rescues soft palate defects in Osr2-Cre; Runx2 fl/fl mice Based on the complimentary expression patterns of Runx2 and Twist1 in the soft palate, we hypothesized that they may oppose each other in regulating their common downstream targets which are important for fate determination of CNC-derived cells. Therefore, we performed ATAC-seq and found that both Runx2 and Twist1 binding sites are present in the regulatory region located around 15-40 kb downstream of the genetic locus of perimysial marker Aldh1a2 ( Figure 6A), suggesting that both Runx2 and Twist1 might directly regulate the expression of Aldh1a2. However, as the binding sites of those two transcription factors are more than 20 kb apart, they are likely to regulate Aldh1a2 independently. Based on this finding, we sought to investigate whether haploinsufficiency of Twist1 may rescue the soft palate defects in Osr2-Cre;Runx2 fl/fl mice by generating Osr2-Cre; Runx2 fl/fl ;Twist1 fl/+ mice. Histological analysis confirmed that the palatal stromal mesenchyme, pterygoid plate, and muscle defects were all indeed rescued in these mice. None of the five newborn Osr2-Cre;Runx2 fl/fl ;Twist1 fl/+ pups we collected had palatal clefts, compared to the 50% penetrance of cleft soft palate in Osr2-Cre;Runx2 fl/fl mice ( Figure 6B-G). Pterygoid plate height was restored in Osr2-Cre;Runx2 fl/fl ;Twist1 fl/+ mice ( Figure 6B-D,H), and muscle fiber orientation and muscle size were also recovered ( Figure 6E-G). To confirm whether the rescue of these muscle defects was due to restoration of perimysial genes, we performed in situ hybridization expression analysis of the soft palate of Osr2-Cre;Runx2 fl/fl ;Twist1 fl/+ mice at E14.5 ( Figure 6J-R). Aldh1a2 was downregulated in the region ( Figure 6K) where Runx2 was deleted ( Figure 6N) and Twist1 expression was expanded in the soft palate of Osr2-Cre;Runx2 fl/fl mice ( Figure 6Q). Compared to Osr2-Cre;Runx2 fl/fl mice, the expression of Aldh1a2 was restored in the Osr2-Cre;Runx2 fl/fl ;Twist1 fl/+ mice at E14.5 ( Figure 6L). This suggests that Runx2 and Twist1 exhibit opposite regulatory effects on the expression of Aldh1a2 in a subset of CNC-derived cells, which may be important for regulating muscle differentiation in the soft palate. Since haploinsufficiency of Twist1 in Osr2-Cre;Runx2 fl/fl ;Twist1 fl/+ mice rescues the expression of Aldh1a, it is most likely that Runx2 activates the expression of Aldh1a2 through repressing Twist1 instead of directly activating Aldh1a2 during normal soft palate muscle development ( Figure 6-figure supplement 1).

Discussion
CNC-derived cells are essential for craniofacial musculoskeletal development, as they give rise to multiple hard and soft tissues in the system, and guide muscle development (Heude et al., 2011;Sugii et al., 2017;Tzahor, 2015). These multiple roles are likely achieved by different subtypes of CNC-derived cells. The heterogeneity of CNC-derived cells has long been studied in the palate based on its anatomical structures along the anterior-posterior, mediolateral, and oral-nasal axes (Bush and Jiang, 2012;Han et al., 2009;Li et al., 2017;Potter and Potter, 2015). To date, understanding of the molecular heterogeneity in different regions of the soft palate mesenchyme has mainly been based on location-specific genes that control local development and signal induction during palate outgrowth (Han et al., 2009;Potter and Potter, 2015). However, this does not    Figure 6 continued on next page completely explain how molecular heterogeneity contributes to the multiple roles played by CNCderived cells during palate formation, or specifically how they guide soft palate muscle development.
In this study, we have revealed the cellular-level heterogeneity in the soft palate and established that different subtypes of CNC-derived cells are associated with distinct differentiation potentials and functions. Functional analysis of each CNC-derived cluster shows previously unknown subtypes of CNC-derived cells, and computational analysis suggests previously unknown lineage differentiation trajectories. Tfap2b+ cells are the least differentiated subtype and give rise to Aldh1a2+ perimysial progenitor cells and Top2a+ transitioning cells, which further differentiate into perimysial fibroblasts and midline mesenchymal fibroblasts. Tfap2b+ common progenitors, transitioning cells and some Aldh1a2+ perimysial progenitor cells are transiently present only at early stages of soft palate development, consistent with their roles as progenitors. This transient presence of CNCderived progenitor cells is similar to that of neural crest cells, which are known to be pluripotent during embryonic development and disappear at later stages (Bronner and LeDouarin, 2012). In addition, we show that markers labeling soft palate perimysial populations, such as Hic1, Aldh1a2, and Tbx15, are also expressed by connective tissues in other craniofacial muscles, including the tongue and masseter muscles. Consistent with our findings, Aldh1a2 is known as an important enzyme for retinoic acid signaling, which is crucial for neural crest cells as they guide the positioning of extraocular muscles (Matt et al., 2008). The potential myogenic-supportive function of Hic1 has been confirmed by the recent finding that the Hic1+ population represents a source of quiescent mesenchymal progenitors that play important roles during the regeneration of skeletal muscles in the limbs (Scott et al., 2019). Hence, Aldh1a2 and Hic1 might be novel markers for CNC-derived perimysial tissues, which may perform important pro-myogenic functions during muscle development.
During embryonic development, the palatal shelves grow in a lateral-to-medial direction both before and after their elevation (Bush and Jiang, 2012). Consistent with this, our in vivo analysis shows that the least differentiated Tfap2b+ subpopulation is located in the lateral region of the soft palate during the early stages of its development, and the perimysial and midline mesenchymal populations reside in central myogenic sites and the medial region of the soft palate. Our results have revealed complex cellular heterogeneity and a differentiation hierarchy of cell populations that contribute to the craniofacial musculoskeletal system, which will require further analysis.
Recently, studies have found that several transcription factors regulate the development of different components of a musculoskeletal complex in a coordinated fashion to form a functional unit (Colasanto et al., 2016;Hasson et al., 2010;Mathew et al., 2011;Vickerman et al., 2011). Our study shows that Runx2 is expressed in CNC-derived cells involved in early cell fate determination and in perimysial cells in the soft palate mesenchyme. Loss of Runx2 in CNC-derived cells of the soft palate mesenchyme leads to multiple tissue defects in the soft palate, including fibrous tendon tissue, soft palate cleft and muscle defects. There is a fate change of CNC-derived cells from perimysial cells to midline mesenchymal cells in the soft palate of Runx2 mutant mice. As perimysial cells are closely associated with muscle development, loss of Runx2 affecting their differentiation may in turn affect their secretion of signaling cues that promote muscle proliferation and differentiation. Indeed, we show that multiple genes associated with pro-myogenic secreted factors specifically expressed by perimysial populations, such as Aldh1a2, Igf1, Cxcl12, and Cthrc1, are downregulated in Runx2 mutant mice (Matt et al., 2008;Schiaffino and Mammucari, 2011;Spector et al., 2013; Source data 1. Source data for Figure 6H.  Vasyutina et al., 2005). Our study provides clues as to how those transcription factors might play different roles in regulating multiple musculoskeletal system components, but how the development of multiple components is integrated still needs further investigation.
Transcription factors often regulate different downstream targets in distinct tissues. Previous studies have shown that Runx2 regulates the differentiation of CNC-derived cells during early tooth and intramembranous bone formation through distinct sets of downstream targets including Gli1, Lef1,Tcf1,Wnt10a,Wnt10b,and Tgfb1 in osteogenic cells and Dusp6,Enpp1,Igfbp3,and Fgf3 in dental mesenchyme (James et al., 2006). In this study, we have shown that downstream targets of Runx2 have differing responses to the loss of Runx2 in the soft tissue. Perimysial markers were specifically downregulated upon loss of Runx2, while genes expressed specifically in the midline mesenchymal cells and a set of more broadly expressed genes are upregulated. Although Runx2 has both transcriptional activation and repression domains, its different regulatory effects on distinct downstream targets in the soft palate mesenchyme could be direct or indirect. Our results thus reveal previously unknown roles of Runx2 in muscle development and help to elucidate the tissue-specific regulatory mechanisms by which Runx2 guides development.
Twist1 suppresses the function of its binding partner Runx2 through blocking the DNA binding domain of Runx2 to inhibit osteoblast differentiation and promote chondrocyte maturation (Bialek et al., 2004;Hinoi et al., 2006). In this study, we reveal complimentary expression patterns of Runx2 and Twist1 in the perimysial and midline mesenchymal populations during soft palate development, which seems to confirm their antagonistic interaction (Bialek et al., 2004). However, in contrast to the previously reported model of Twist1 and Runx2 antagonistic interaction, we show that loss of Runx2 in CNC-derived cells leads to abnormal upregulation of Twist1 in the perimysial population. Further analysis has shown that suppression of Twist1 in the perimysial population by Runx2 is necessary to maintain the expression level of perimysial marker gene Aldh1a2, which may be important for regulating muscle development. Our findings thus reveal a novel mechanism of Runx2-Twist1 genetic interaction that integrates the development of different types of CNC-derived cells with muscles to guide them to form a functional unit in the soft palate.
In summary, our study reveals a complex cellular heterogeneity within the developing soft palate and demonstrates that distinct subpopulations of CNC-derived cells are associated with distinct functions, which coordinate to form intricately connected components of the oropharyngeal complex. Moreover, the regulation of myogenesis by perimysial CNC-derived cells through Runx2-Twist1 interaction in the soft palate might also be shared by other craniofacial musculoskeletal structures. Our study highlights the complex regulatory roles of CNC-derived cells in the development of craniofacial musculoskeletal systems and provides knowledge that may lead to new strategies for craniofacial muscle regeneration.

Immunofluorescence staining
Sections were processed with antigen-retrieval buffer (Vector Labs, Burlingame, CA, H-3300) for 15 min at 100˚C, followed by 1% triton (Sigma Aldrich, St. Louis, MO, T8787) treatment for 10 min at room temperature. Afterwards, sections were incubated with blocking reagent (PerkinElmer, Waltham, MA, FP1012) for 1 hr at room temperature and the primary antibody overnight at 4˚C. Alexaconjugated secondary antibodies were used to show the fluorescence signal at 1:200 dilution. For myoblast determination protein 1 (MyoD), poly HRP-labeled goat anti-mouse IgG (ThermoFisher Scientific, Waltham, MA, B40912) was used as a secondary antibody and Alexa Fluor 488/594 Tyramide SuperBoost kit (PerkinElmer, Waltham, MA,NEL771B001KT, NEL774001KT) were used to develop the signal. Sections were counterstained with DAPI and imaged using a Leica DMI 3000B.

Single-cell RNA sequencing
Soft palate tissue (posterior third of the palatal region) was digested from E13.5, E14.5, and E15.5 controls and E14.5 Osr2-Cre;Runx2 fl/fl embryos by TrypLE express enzyme (Thermo Fisher Scientific, Waltham, MA) at 37˚C with shaking at 600 rpm for 20 min. Single-cell suspension was prepared according to the 10X Genomics sample preparation protocol. Seventeen thousand cells were loaded into the 10X Chromium system and prepared for single-cell library construction using the 10X Genomics Chromium single cell 3' v3 reagent kit. Sequencing was performed on the Novaseq 6000 platform (Illumina, San Diego, CA). Library quality control, sequence alignment, and read counts were analyzed using the CellRanger pipeline version 3.0.2. Raw read counts from each single cell in each sample were analyzed using Seurat R package (Stuart et al., 2019). Cell clusters and variably expressed genes in each cluster were identified by using Log Normalize, Find Variable Genes, Scale Data, and RunPCA functions. Seurat three package was used to combine the single-cell data from three stages as well as E14.5 control and Osr2-Cre;Runx2 fl/fl embryos to perform the integration analysis. Shared variances between different datasets were identified using the function FindIntegra-tionAnchors, then Seurat objects were processed using IntegrateData function. Scaledata, PCA, and UMAP visualization were then used for downstream analysis and visualization. Pseudotime trajectory analysis was done by Monocle three using Seurat 3 UMAP embedding to show cell fate restriction of CNC-derived soft palate mesenchymal cells across three development stages. Gene ontology and pathway analysis of enriched genes in different CNC-derived clusters was performed using Ingenuity Pathway Analysis (QIAGEN. Inc, Hilden, Germany).

ATAC-seq
Single-cell suspension was prepared from the soft palate of E13.5-E14.0 control mice as described above and processed to generate ATAC-seq libraries according to a published protocol (Buenrostro et al., 2015). Sequencing was performed on the NextSeq 500 platform (Illumina, San Diego, CA). ATAC-seq reads were aligned to the UCSC mm10 reference genome using BWA-MEM (Li, 2013). ATAC-seq peaks were called by MACS2 (Zhang et al., 2008). Peaks were annotated and known transcription factor binding motifs were analyzed in the ATAC-seq peaks by HOMER (Heinz et al., 2010).

RNA sequencing
Soft palate tissue was collected from control and Osr2-Cre;Runx2 fl/fl embryos at E14.5. mRNA was isolated using RNeasy Micro Kit (QIAGEN,Hilden,Germany,74404). Samples with RNA integrity number (RIN) >9.0 were used for cDNA library construction and sequencing by UCLA Technology Center for Genomics and Bioinformatics. Pair-end reads with 150 cycles sequencing were performed on Illumina NextSeq 500 platform. Sequence reads were trimmed and aligned using STAR (version 2.6.1d) using mm10 as the reference genome. Read counts were normalized using the upper quartile and differential expression was calculated using gene-specific analysis on the Partek Flow platform (Partek Inc, St. Louis, MO).

Statistical analysis
T-tests were performed for statistical analysis using GraphPad Prism 7. Statistical data are presented as mean ± SEM.