TGF-β signaling and Creb5 cooperatively regulate Fgf18 to control pharyngeal muscle development

The communication between myogenic cells and their surrounding connective tissues is indispensable for muscle morphogenesis. During late embryonic development in mice, myogenic progenitors migrate to discrete sites to form individual muscles. The detailed mechanism of this process remains unclear. Using mouse levator veli palatini (LVP) development as a model, we systematically investigated how a distinct connective tissue subpopulation, perimysial fibroblasts, communicates with myogenic cells to regulate mouse pharyngeal myogenesis. Using single-cell RNAseq data analysis, we identified that TGF-β signaling is a key regulator for the perimysial fibroblasts. Loss of TGF-β signaling in the neural crest-derived palatal mesenchyme leads to defects in perimysial fibroblasts and muscle malformation in the soft palate in Osr2Cre;Tgfbr1fl/fl mice. In particular, Creb5, a transcription factor expressed in the perimysial fibroblasts, cooperates with TGF-β signaling to activate expression of Fgf18. Moreover, Fgf18 supports pharyngeal muscle development in vivo and exogenous Fgf18 can partially rescue myogenic cell numbers in Osr2Cre;Tgfbr1fl/fl samples, illustrating that TGF-β-regulated Fgf18 signaling is required for LVP development. Collectively, our findings reveal the mechanism by which TGF-β signaling achieves its functional specificity in defining the perimysial-to-myogenic signals for pharyngeal myogenesis.


Introduction
Inductive tissue-tissue interaction constitutes a key mechanism for organ development and tissue homeostasis (Souilhol et al., 2016;Takahashi et al., 2013;Vainio and Müller, 1997). The interaction and communication of myogenic cells with the surrounding connective tissues are also indispensable for muscle morphogenesis. Although the detailed cellular and molecular aspects of this dialogue remain to be elucidated, the connective tissues have been shown to serve as a signal source and 'pre-pattern' for muscle differentiation and patterning during development (Michailovici et al., 2015;Nassari et al., 2017).
At the onset of pharyngeal muscle development, myogenic progenitors migrate from cranial paraxial mesoderm into the center of each pharyngeal arch to form the primary myogenic sites surrounded by cranial neural crest (CNC)-derived cells. Later, in coordination with the development of other CNC-derived craniofacial connective tissues, myogenic progenitor cells in these mesodermal cores migrate to the final myogenic sites while being segregated and patterned into discrete muscle masses (Noden and Francis-West, 2006;Sambasivan et al., 2011;Shih et al., 2008;Ziermann et al., 2018). A distinct regulatory network upstream of myogenesis has been identified during earlystage pharyngeal muscle development (Buckingham and Rigby, 2014;Sambasivan et al., 2011), but the regulatory mechanism that establishes the fine-tuned craniofacial muscle anlagen in their final myogenic sites remains unknown (Noden and Francis-West, 2006).
In the final stage of pharyngeal muscle formation, the first and second pharyngeal arch mesodermal cores respectively contribute to the muscles responsible for mastication and facial expression; the muscles in the oropharyngeal region crucial for breathing, swallowing, and speaking are derived predominantly from the fourth pharyngeal arch, including all the pharyngeal constrictors and soft palatal muscles except the tensor veli palatini (Frisdal and Trainor, 2014;Li et al., 2019;Michailovici et al., 2015;Shiba and Chhetri, 2019;Sugii et al., 2017). Notably, interactions between CNCderived and myogenic cells become more prominent as the CNC-derived cells start to contribute to connective tissues while individual muscles segregate and proceed to their final locations (Ziermann et al., 2018).
During this late pharyngeal muscle development in mice, most pharyngeal muscles of the head have already segregated and reached their ultimate locations between E11 and E13.5 (Noden and Francis-West, 2006). In contrast, the soft palatal muscle anlage development starts later, concurrent with the soft palatal shelf primordium formation (Han et al., 2021;Li et al., 2019). From E13.5 onwards, individual soft palatal muscles, particularly the levator veli palatini (LVP), can be clearly detected segregating from the rest of the pharyngeal muscles and populating the soft palate primordium, patterning and migrating towards the midline in parallel to the lateral-to-medial outgrowth of the connective tissues consisting of the CNC-derived palatal mesenchymal cells (Grimaldi et al., 2015;Han et al., 2021;Li et al., 2019). Therefore, soft palatal muscle development provides an optimal opportunity to study the initiation, segregation, and migration of muscle anlage formation during craniofacial development.
Clinically, the most common congenital craniofacial malformation is cleft lip with or without cleft palate (Vyas et al., 2020). When the soft palate is affected, disruption of the soft palatal shelves separating the oral and nasal cavities co-occurs with defects in soft palatal muscle formation and oropharyngeal function (Li et al., 2019;Carvajal Monroy et al., 2012). Furthermore, mouse models have demonstrated that loss of key regulators in the perimysial cells, a subset of the palatal mesenchymal connective tissues closely associated with the muscles, ultimately affects soft palatal muscle differentiation and patterning (Han et al., 2021;Sugii et al., 2017). These findings suggest that CNCderived perimysial cells possibly contribute to part of the embryonic myogenic niche for craniofacial muscle development at the later stages. Thus, soft palate development can serve as a useful model for studying how perimysial cells regulate late pharyngeal muscle anlage formation through tissue-tissue interactions.
In the present work, we investigated the mechanism that regulates perimysial-to-myogenic communication using the development of the pharyngeal muscle LVP as a model. Using unbiased single-cell RNAseq (scRNAseq) data analysis combined with mouse genetic approaches, we identified TGF-β signaling as a predominant and specific signaling activity that enables the perimysial cells to interact with the adjacent myogenic cells during late pharyngeal muscle development. Furthermore, we identified that a perimysial regulator Creb5 cooperates with TGF-β signaling to co-activate the expression of specific perimysial-to-myogenic signals such as Fgf18 to regulate pharyngeal myogenesis. Taken together, our findings reveal that TGF-β signaling and perimysial-specific regulators may cooperatively define perimysial-to-myogenic signals in regulating pharyngeal myogenesis.

Results
Unbiased screening using scRNAseq data analysis identifies TGF-β signaling as predominant and specific signaling for the perimysial cells within the soft palatal mesenchyme During late pharyngeal muscle development, myogenic cells of the soft palatal muscles, particularly those of the LVP, can be observed segregating from the existing myogenic site associated with middle pharyngeal constrictor to populate the soft palatal shelves and pattern alongside the CNCderived palatal mesenchymal cells (Han et al., 2021;Li et al., 2019). In parallel, during this stage, our recent scRNAseq analysis identified that the morphologically homologous palatal mesenchyme is also patterned into heterogeneous cellular domains associated with distinct anatomical locations and gene expression profiles (Han et al., 2021). We thus hypothesized that the perimysial population, anatomically adjacent to the myogenic cells, serves as a microenvironment or niche for embryonic pharyngeal muscle development. Furthermore, there must be a specific regulatory mechanism that differentiates the perimysial cells from the rest of the palatal mesenchymal cell populations so that the perimysial cells can provide specific signals to define the myogenic site.
Since signaling pathways have been reported to play critical roles in the patterning and fate determination of craniofacial mesenchymal cells (Mishina and Snider, 2014;Neubüser et al., 1997;Vincentz et al., 2016;Xu et al., 2019), we performed an unbiased signaling activity analysis of the cell populations of the soft palate, focusing on identifying key signaling pathways specifically activated in the perimysial cells. Using CellChat's incoming signaling analysis from scRNASeq data of E13.5-E15.5 soft palate, we identified inferred signaling pathway activities of different cell types based on enriched signaling interactions ( Figure 1A-B, Figure 1-figure supplement 1). Notably, the signaling activity of all the palatal mesenchymal cell types was categorized into the same group due to the shared expression of similar signaling pathways ( Figure 1B), consistent with their common connective tissue identity. Among this group of pathways, Ephrin-Eph (EphA/EphB) (Benson and Serrano, 2012;Xavier et al., 2016), non-canonical Wnt (ncWnt) (Reynolds et al., 2019), IGF (Marchant et al., 2020), TGF-β (Iwata et al., 2011), FGF (Nie et al., 2006), and Hedgehog signaling (HH) pathways (Everson et al., 2018) were previously identified to be required for the development of craniofacial mesenchymal tissues including palatal mesenchyme, so we reasoned that they are more likely to be functionally required by the perimysial cell fate determination, too. We thus next analyzed the intensities of activities of these pathways in the perimysial cells and found that the TGF-β signaling activity was inferred to be the highest ( Figure 1C). Moreover, a known TGF-β signaling downstream target, Tgfbi (transforming growth factor-beta-induced protein), was among the top 10 conserved marker genes for the perimysial cells across various developmental stages ( Figure 1D), suggesting that TGF-β signaling is not only a predominant signaling pathway but also specifically associated with regulating the perimysial cell fate and maintaining a perimysial-derived signal.
To validate the signaling activity analysis result from scRNAseq data, we next evaluated the TGF-β signaling activity in vivo during different stages of soft palatal development. Consistent with scRNAseq analysis, TGF-β signaling was present in the palatal shelf mesenchyme primordium during soft palatal development from E12.5 onwards (Figure 2A-H) and was most active in the cells surrounding the LVP myogenic cells, anatomically identified as perimysial cells , particularly during the early segregation and separation process at E13.5 ( Figure 2B, F, J, and N) and E14.5 ( Figure 2C, G, K, and O). Notably, TGF-β signaling was ubiquitous in the mesenchymal cells close to the pharyngeal wall at E12.5 prior to the arrival of myogenic cells (Figure 2A, E, I, and M). This specific activation of TGF-β signaling in the presence of perimysial cells suggested that TGF-β signaling may play a role in patterning perimysial cells to establish a myogenic niche in the specific anatomical location for the LVP. This postulation is supported by our previous work, which showed that loss of TGF-β in Wnt1-Cre;Tgfbr1 fl/fl mice in the early pre-migratory neural crest cells and their derivatives from E9.5 onward affect proliferation and differentiation of the tongue and other craniofacial muscles (Chai et al., 2000;Han et al., 2014).
E13.5-E15.5 soft palatal scRNAseq analysis led to the inference that the perimysial cells contain at least two distinct subpopulations with different differentiation statuses, namely earlier perimysial progenitors (Aldh1a2+) and more committed perimysial fibroblasts, which express markers other than Aldh1a2, such as Tbx15 (Han et al., 2021). We therefore sought to further specify the cellular identity    of the region with TGF-β signaling activity. After subclustering the perimysial cell population from the E13.5-E15.5 soft palatal scRNAseq data, we found these cells could be further divided into 3 clusters ( Figure 2-figure supplement 1A). Since Aldh1a2 expression was most enriched in cluster 2 ( Figure 2-figure supplement 1B), this cluster appeared to be associated with the perimysial progenitor cells, whereas clusters 0 and 1 appeared to be associated with the perimysial fibroblasts. Within the two perimysial fibroblast clusters, the previously reported marker Tbx15 was expressed by most cells in cluster 0, but not by the majority of cells in cluster 1 (Figure 2-figure supplement 1C). Instead, we identified that another perimysial marker, Smoc2, labeled both clusters, including the Tbx15-negative cells (Figure 2-figure supplement 1D). We confirmed that Smoc2  Loss of TGF-β leads to defects in soft palatal shelf mesenchyme and LVP in Osr2 Cre ;Tgfbr1 fl/fl mice To test whether this specific TGF-β signaling in the perimysial cells in later developmental stages also plays any functional role in regulating the fate of myogenic cells, we utilized Osr2 Cre , which affects the palatal mesenchyme at a later stage than Wnt1-Cre (Chen et al., 2009). Using Osr2 Cre ;Rosa26 LSL-tdTomato mice, we confirmed that Osr2 Cre can specifically label the palatal mesenchyme without affecting the myogenic cells during soft palatal muscle development ( Figure 3-figure supplement 1A-D). Thus, we generated Osr2 Cre ;Tgfbr1 fl/fl mice, in which TGF-β signaling is ablated in the palatal mesenchyme due to the loss of TGF-β Type 1 receptor Alk5 (Tgfbr1). Loss of TGF-β signaling in the palatal mesenchyme, including the perimysial cells, led to morphologically deformed palatal shelves and cleft palate in the LVP region in E18.5 Osr2 Cre ;Tgfbr1 fl/fl mice ( Figure 3B and D), compared with the intact control palatal shelves ( Figure 3A and C), indicating intrinsic defects in the palatal mesenchymal cells. In addition, while abundant MHC+ muscle fibers were present in the control palatal shelves ( Figure 3E, I, G, and K), in the defective palatal shelves of Osr2 Cre ;Tgfbr1 fl/fl mice, LVP formation was so severely affected that no muscle fibers could be detected at this stage ( Figure 3F, J, H, and L). This confirmed that TGF-β signaling is indispensable for the perimysial-derived signaling for myogenesis. We next evaluated the progression of the palate mesenchymal and myogenic cell defects in the Osr2 Cre ;Tgfbr1 fl/fl mice. We found that the palatal shelves of Osr2 Cre ;Tgfbr1 fl/fl mice started to appear morphologically different from those of the controls at around E14.5, suggesting that defects in the palatal mesenchymal cells start at this stage (        TGF-β signaling is required for perimysial fibroblasts in regulating pharyngeal muscle development To precisely identify the molecular changes of the perimysial population in the Osr2 Cre ;Tgfbr1 fl/fl mice, we performed scRNAseq analysis to compare cell-type-specific gene expression profiles of cell populations in E14.5 Osr2 Cre ;Tgfbr1 fl/fl and control soft palates. Using integration analysis based on shared variance ( Figure 4-figure supplement 1A), we identified similar clusters in the soft palates of control and Osr2 Cre ;Tgfbr1 fl/fl mice at this stage. We first distinguished the palatal mesenchymal cells (Runx2+/ Twist1+) from other non-mesenchymal cell types in the soft palate ( Figure 4-figure supplement 1B-L). Then, using more detailed markers we recently established from E13.5-E15.5 soft palatal scRNAseq analysis (Han et al., 2021), we further distinguished the perimysial population (Tbx15+/ Hic1+/Aldh1a2+) from the remaining palatal mesenchymal populations of the E14.5 Osr2 Cre ;Tgfbr1 fl/fl and control mice ( Figure 4-figure supplement 1M-N).
To more specifically evaluate the changes in the perimysial cells following the loss of TGF-β signaling, we subsetted the perimysial population from the mesenchymal clusters of the integrated Osr2 Cre ;Tgfbr1 fl/fl and control soft palate scRNAseq data at E14.5 ( Figure 4A). We first compared the expression of perimysial progenitor and perimysial fibroblast markers between the Osr2 Cre ;Tgfbr1 fl/fl and control perimysial cells at E14.5 ( Figure 4B-D), to identify whether one or both subpopulations of the perimysial cells were affected. Consistent with the overlapping of perimysial fibroblast populations in the region of active TGF-β signaling, in our scRNAseq data we found reduced expression of perimysial fibroblast markers Tbx15 and Smoc2 in the perimysial population in E14.5 Osr2 Cre ;Tgfbr1 fl/fl soft palates relative to controls ( Figure 4B and C), while the Aldh1a2+ perimysial progenitor subpopulation was not affected ( Figure 4D). More specifically, in vivo expression patterns further showed that Tbx15+ cells were located predominantly in between the LVP myogenic cells, while Smoc2 expression was present mostly around the periphery of the myogenic site ( Figure 4F and H). Consistent with the scRNAseq analysis, expression of Tbx15 and Smoc2 was affected in perimysial cells in Osr2 Cre ;Tgfbr1 fl/ fl mice at these respective sites ( Figure 4G and I). Also consistent with the scRNAseq analysis, the Aldh1a2+ cells located mainly adjacent to the middle pharyngeal constrictor in the control were not affected in the Osr2 Cre ;Tgfbr1 fl/fl mice ( Figure 4J and K). While our recent study showed that Runx2 is a regulator of Aldh1a2+ early perimysial progenitors (Han et al., 2021), here we identified TGF-β signaling as the main regulator for the fate determination of both populations of the more committed perimysial fibroblasts, suggesting a specific role of TGF-β signaling in late stage-pharyngeal muscle development.
Since TGF-β signaling has been reported to regulate proliferation of the hard palate mesenchyme (Ito et al., 2003), we further investigated whether the loss of Tbx15 and Smoc2 at E14.5 in the Osr2 Cre ;Tgfbr1 fl/fl mice was due to the failure of palatal mesenchymal cells to survive, differentiate, and/or proliferate. To evaluate the activation of these two markers during development, we analyzed their expression patterns at E13.5-E14.0 in control and Osr2 Cre ;Tgfbr1 fl/fl mice. We found that Tbx15 expression was not activated at early stages: Tbx15 was not detectable E13.5, and only starting to be activated at a low level by a few cells at E14.0 in both and Osr2 Cre ;Tgfbr1 fl/fl mice (Figure 4-figure supplement 2A-D). This low level of Tbx15 increased dramatically at E14.5 in the control ( Figure 4F) but failed to increase in a similar manner in Osr2 Cre ;Tgfbr1 fl/fl mice ( Figure 4G), suggesting that the reduced Tbx15 expression at E14.5 was due to a failure of activating the differentiation of Tbx15+ cells. In parallel, we found that Smoc2 was expressed extensively in the palatal mesenchyme at E13.5 at comparable levels in control and Osr2 Cre ;Tgfbr1 fl/fl mice, and its expression started to reduce in the  CellChat analysis of soft palatal scRNAseq data identifies putative perimysial-to-myogenic signaling molecules downstream of TGF-β signaling Since the perimysial fibroblast population is anatomically adjacent to the myogenic cells, they are most likely to be part of a microenvironment that supports myogenesis through signaling communication. To systemically identify potential signaling interactions between these two cell populations, we examined the E13.5-E15.5 integrated soft palatal scRNAseq data with the outgoing signaling analysis of the CellChat package (Han et al., 2021;Jin et al., 2021) to infer enriched intercellular signaling interactions at the single-cell level. In particular, we focused our analysis on signals sent from perimysial fibroblasts and received by myogenic cells. In doing so, we detected several enriched interactions including those between perimysial-derived signaling molecules in ncWnt, FGF, Notch, and BMP signaling pathways (Wnt5a, Fgf18, Dlk1, and Bmp4) and corresponding receptors expressed in the myogenic cells (Fdz4, Fgfr1/4, Notch3, and Bmpr1a/Acvr2b, respectively) ( Figure 5A). We also identified other interactions that may be associated with cell behavior through modulating extra-cellular matrices, such as Thrombospondin (Thbs3-Sdc1), Pleiotrophin (Ptn-Ncl), Nectins (Nectin3-Nectin1), Midkine (Mdk-Ncl), and Laminins (Lama4/b1/c1-Dag1) ( Figure 5A).
To identify perimysial-to-myogenic interactions that are functionally required for regulating myogenic fate, we decided to narrow our selection to the signaling molecules more specifically expressed by the perimysial fibroblasts. By checking the expression patterns of these signaling molecules in E13.5-E15.5 soft palatal scRNAseq data, we found that Lama4, Fgf18, and Dlk1 are the top three molecules more enriched in the perimysial fibroblasts than in other cells ( Figure 5B). This finding suggests they play a more specific role in regulating myogenesis. Among these three molecules, expression of Fgf18 and Lama4 was found to be reduced in the Osr2 Cre ;Tgfbr1 fl/fl perimysial cells, whereas Dlk1 expression appeared unaffected ( Figure 5C-F). This suggests that TGF-β signaling could also regulate the expression of perimysial-derived Lama4 and Fgf18 during muscle development. In vivo expression of these two genes further showed that Fgf18 expression is more restricted to the region of perimysial fibroblasts ( Figure 5H) than that of Lama4 ( Figure 5I), similar to the TGF-β signaling activity. Thus, we identified Fgf18 as the most likely perimysial-to-myogenic signal regulated by TGF-β signaling to support myogenesis.
To better understand the role of perimysial-myogenic signals during myogenesis, we further investigated the in vivo expression of Fgf18 in perimysial cells during LVP development in control and Osr2 Cre ;Tgfbr1 fl/fl mice. We found consistently enriched expression of Fgf18 in the location of perimysial fibroblasts most adjacent to the myogenic cells between E13.5-E16.5 ( Figure 5J        This strong correlation between the perimysial-fibroblast-derived Fgf18 and Myod1+ myogenic cells further supported the potential role of Fgf18 in myogenesis. We therefore next evaluated the expression of FGF receptors in the MyoD+ myogenic cells to test whether myogenic cells could receive an Fgf18 signal. Consistent with CellChat predictions, perimysial fibroblast-derived Fgf18 signals were more likely to communicate with Fgfr1 and/or Fgfr4 in the myogenic cells, which did not express Smad2/3 and Creb5 cooperatively regulate Fgf18 signaling for pharyngeal muscle development To investigate how TGF-β signaling serves as a key regulator of perimysial-to-myogenic communication, particularly Fgf18 signaling, we first explored whether TGF-β signaling mediators Smad2/3 can directly bind to the promoter region of Fgf18. Using JASPAR TF binding sites (TFBSs) predictions incorporated in the UCSC Genome Browser (Castro-Mondragon et al., 2022), we identified that several putative TFBSs for Smad2/3 are present in the promoter region 2 kb upstream of the transcription start site (TSS) ( Figure 6A) and focused on the binding motif at chr11:33098353-33098362, which had the highest predicted binding score, since it is most likely to be the binding site ( Figure 6A). Using a Cut and Run assay for Smad2/3, we identified significantly greater binding Smad2/3 to DNA fragments around this binding motif when compared to IgG control ( Figure 6B), suggesting direct binding of Smad2/3 to this predicted TFBS in the Fgf18 promoter region.
To further test the functional requirment of this TFBS, we repressed the transcription initiation of this TFBS using a CRISPRi system that sterically prevents the association between DNA motifs and transcription factors when the gRNA is targeting the promoter region (Larson et al., 2013), which Source data 1. Source data for Figure 5U.
Source data 2. Source data for Figure 5V.  have also been used for functional analysis of transcription binding sites . We designed a CRISPRi plasmid (Fgf18 CRISPRi Plasmid) whereby the gRNA guides the CRISPRi complex to target only this Smad2/3 binding site, as other binding sites are located distantly, hence blocking its transcriptional activity. Since gRNAs for CRISPRi are reported to have an optimal target effect within a range of 150 bp (MacLeod et al., 2019), we examined gRNA targeting sites from +150 bp to -150 bp around the TFBS and identifed a gRNA with the highest targeting efficency and lowest off-target rate, targeting very close to (47 bp downstream of) the TFBS. Next, we transfected the Fgf18 or control CRISPRi plasmid into the cultured soft palatal mesenchymal cells and compared the response of these cells to TGF-β1 treatment. Since cellular response to TGF-β ligands can be transient rather than accumulative (Sorre et al., 2014), we first identified the optimal time to analyze Fgf18 expression levels in the soft palatal mesenchymal cells following TGF-β1 treatment. Consistent with previous studies which found that pSmad2 activity starts to decline around 5 hr following TGF-β1 treatment (Inman et al., 2002), we found that 4 hr of TGF-β1 treatment led to a significant increase of Fgf18 expression, which reduced when analyzed 18 hr post-treatment ( Figure 6C). Thus, 4 hr post-treatment was used as the time point for evaluating the cellular response to the TGF-β1 treatment in the following experiments. Compared with the control CRISPRi plasmid, Fgf18 CRISPRi plamid significantly attenuated the increase in expresssion of Fgf18 following TGF-β1 treatment ( Figure 6D). This finding suggested that the binding to this TFBS in the Fgf18 promoter region is functionally required for the TGF-β-Smad2/3 signaling cascade to directly activate Fgf18 expression.
The biological effects of TGF-β signaling are contextual (Morikawa et al., 2016). One mechanism for cell-type-specific response to TGF-β signaling is through the cooperation with cell-type-specific master transcription factors, which enable TGF-β signaling to specifically affect the cell-type-specific genes bound by these master regulators (Mullen et al., 2011). To identify master regulators that help TGF-β signaling activate a perimysial-specific response in palatal mesenchymal cells, we performed Single-Cell Regulatory Network Interference and Clustering (SCENIC) analysis using E13.5-E15.5 calculations. Data is presented as mean ± SEM. (C) qPCR analysis of Fgf18 expression in E14.5 soft palatal mesenchymal cell culture after treatment with 5 ng/ml TGF-β1 or BSA compared after 4 or 18 hr. Note the increase of Fgf18 at 4 hr is higher than at 18 hr. *, p≤0.05; ***; p≤0.001. Statistical significance was assessed by unpaired t-test with two-tailed calculations. Data is presented as mean ± SEM. (D) qPCR analysis of Fgf18 expression of E14.5 soft palatal mesenchymal cell culture transfected with Fgf18 or Control CRISPRi plasmid followed by treatment with 5 ng/ml TGF-β1 or BSA for 4 hr. '+' or '-' under the plots indicates the presence or absence of the indicated treatment. **, p≤0.01; ***, p≤0.001. Statistical significance was assessed by ANOVA. Data is presented as mean ± Source data 1. Source data for Figure 6B.
Source data 2. Source data for Figure 6C.
Source data 3. Source data for Figure 6D.
Source data 4. Source data for Figure 6O.
Source data 5. Source data for Figure 6P.
Source data 6. Source data for Figure 6Q.     (Aibar et al., 2017;Han et al., 2021). We identified that gene regulatory networks (GRNs) mediated by Hic1, Creb5, Meox1, Etv5, Tbx15, Bach2, and Nr2f2 are most associated with perimysial fibroblasts (Figure 6-figure supplement 1). The expression patterns of these genes in scRNAseq data indicated that Hic1, Creb5, Meox1, Etv5, Tbx15, and Bach2 are more enriched and specific to the perimysial fibroblasts than other palatal mesenchymal populations, and are thus more likely to be associated with the perimysial fibroblast-specific cell behaviors ( Figure 6E). Next we explored whether these identfied perimysial factors could interact with TGF-β signaling to activate perimysial-specific TGF-β downstream genes such as Fgf18. Among the perimysial-associated factors, Creb5 has been reported to be required for TGF-β signaling to induce Prg4 expression during chondrocyte development . In vivo expression patterns further validated that, among these regulatory genes ( Figure 6I-N), Creb5 was most abundantly expressed in the perimysial fibroblasts ( Figure 6I), similar to that of Fgf18 and TGF-β signaling activity ( Figure 6G and H). Moreover, unlike that of Tbx15, Creb5 expression was present in the palatal mesenchyme as early as E13.5, similar to TGF-β signaling activity and Fgf18 expression ( Figure 6-figure supplement 2A). Since its expression did not change in the mutant from E13.5 to E14.5, it is likely to function as a partner rather than a downstream target of the TGF-β signaling ( Figure 6-figure supplement 2A-D). We thus hypothesized that Creb5 is likely one of the candidates that may function cooperatively with TGF-β signaling to activate perimysial-specific signaling genes such as Fgf18. To test the function of Creb5, we first evlauated its role in regulating Fgf18 expression using E14.5 soft palatal mesenchymal cell culture and found that reduction of Creb5 in soft palatal mesenchymal cells led to reduction of Fgf18 expression ( Figure 6O and P), indicating that Creb5 can regulate Fgf18 expression. To evaluate the cooperation of Creb5 and TGF-β signaling, we compared the increase of Fgf18 following TGF-β1 treatment with or without Creb5 reduction, and found that reducing Creb5 signficantly attenuated the increase of Fgf18 expression in soft palatal mesenchymal cells ( Figure 6Q), suggesting that TGF-β signaling requires cooperation with Creb5 to efficiently activate downstream targets like Fgf18 in the perimysial fibroblasts. Furthermore, we also found that TGF-β signaling activity (pSmad2), Creb5 and Fgf18 are also expressed in a similar pattern in the perimysial cells during the development of the other pharyngeal muscles, such as the first pharyngeal arch-derived masseter muscle (Figure 6figure supplement 3), confirming the possibility that this TGF-β signaling and Creb5 cooperation in regulating perimysial signals, which we identified using the soft palatal muscle model, can plausibly be extended to other craniofacial muscles.
Fgf18 reduction in perimysial cells leads to reduced myogenic cells in the LVP region of Osr2 Cre ;Fgf18 fl/fl mice while Fgf18 restoration partially rescues myogenic defects in Osr2 Cre ;Tgfbr1 fl/fl samples To test whether perimysial-derived Fgf18 is functionally required for myogenesis, we generated Osr2 Cre ;Fgf18 fl/fl mice. Unlike in the Osr2 Cre ;Tgfbr1 fl/fl mice, loss of Fgf18 in the Osr2 Cre ;Fgf18 fl/fl mice did not result in any obvious palatal shelf formation defects compared to the control mice at birth ( Figure 7A-D), consistent with recent studies demonstrating that Fgf18 is not required for the development of mesenchymal cells themselves in the hard palate (Xu et al., 2016;Yue et al., 2021). In contrast, the sizes of the soft palatal muscles of the Osr2 Cre ;Fgf18 fl/fl mice appeared to be smaller than those of controls, which was especially obvious in the posterior part of the LVP ( Figure 7G and  H). This was probably due to the difference in muscle fiber organization in the anterior and posterior parts of the LVP. The LVP consists of muscle fibers oriented both parallel and perpendicularly to the section plane in the anterior region ( Figure 7I) and predominantly perpendicularly in the posterior region ( Figure 7K). While all the fibers appear more sparse in the anterior and posterior LVP of the Osr2 Cre ;Fgf18 fl/fl mice ( Figure 7J and L), the perpendicular fibers are significantly reduced ( Figure 7J). This reduction of LVP muscle size was confirmed by reduced muscle cross-section area ( Figure 7M-Q) and muscle volume ( Figure 7R-V) in the Osr2 Cre ;Fgf18 fl/fl mice compared with the controls. To further analyze the progression of muscle defects, we assessed the soft palatal and myogenic phenotypes from E14.5 onwards (Figure 7-figure supplement 1). While the MyoD+ myogenic cells appeared comparable in the control and Osr2 Cre ;Fgf18 fl/fl LVP at E14.5 (Figure 7-figure supplement 1A-E), the myogenic cells in the Osr2 Cre ;Fgf18 fl/fl mice started to appear progressively reduced compared to the control, and this difference became more apparent at E16.5 (Figure 7-figure supplement 1F-J), confirming the role of Fgf18 in late-stage myogenesis. We also confirmed that the loss of functional   (Figure 7-figure supplement  1K-N). Taken together, these findings suggest that Fgf18 may not be required for determining the fate of palatal mesenchymal cells themselves, but rather may serve as a paracrine molecule for the myogenic cell development.
To further investigate how loss of Fgf18 led to reduce myogenic cells form E14.5 onwards, we analyzed change of proliferation and cell survival of myogenic cells in Osr2 Cre ;Fgf18 fl/fl mice compared with the controls. At this stage, the myogenic cells in Osr2 Cre ;Fgf18 fl/fl mice and controls are rarely apoptotic (Figure 7-figure supplement 2A and D), but display proliferative activity in both mutants and controls (Figure 7-figure supplement 2B and E). Therefore, we focused on analyzing the proliferative status of the myogenic cells. Since Myf5 expression was mostly associated with proliferating cells with enriched Fgfr4 expression, we next compared the proliferative status of the Myf5+ myogenic progenitors specifically in Osr2 Cre ;Fgf18 fl/fl and control mice. We could see a significant reduction of proliferation of Myf5+ cells in Osr2 Cre ;Fgf18 fl/fl mice (Figure 7-figure supplement 2C, F, and G), suggesting that loss of Fgf18 led to reduced myogenic progenitor proliferation in the Osr2 Cre ;Fgf18 fl/ fl mice, which ultimately led to smaller muscle mass. Furthermore, to confirm that Fgf18 could directly work on myogenic cells, we treated C2C12 mouse myogenic cells, which we maintained as undifferentiated myoblasts in growth medium, with FGF18 and found that the proliferation rate of the C2C12 cells significantly increased in growth medium as early as 1 day later (Figure 7-figure supplement  2H-J), eventually leading to significantly increased C12C2 cell numbers at 3 days post-treatment (Figure 7-figure supplement 2K), supporting the role of Fgf18 in directly regulating myogenic progenitor proliferation in vivo.
To test whether Fgf18 treatment could restore some of the severe muscle formation defects following the loss of TGF-β signaling in Osr2 Cre ;Tgfbr1 fl/fl mice, we optimized an organ culture system (Alfaqeeh and Tucker, 2013;Humpel, 2015) to culture 300 μm-thick slices of embryonic head tissues containing the LVP region. This enabled the maintenance of the three-dimensional structure of the pharyngeal region and allowed for efficient nutrient penetration and bead implantation ( Figure 7W). Using this culture system, we found that FGF18 bead treatment in the soft palatal shelves led to increased MyoD+ cells compared with BSA beads in Osr2 Cre ;Tgfbr1 fl/fl slice cultures ( Figure 7X and Y). This suggests that Fgf18 is one of the key regulators of perimysial-to-myogenic signaling affected in the LVP in the Osr2 Cre ;Tgfbr1 fl/fl mice. Taken together, our studies identified that TGF-β signaling interacts with perimysial regulator Creb5 to specify the perimysial-to-myogenic signaling, such as Fgf18, in individual pharyngeal muscle development.
for quantification in Q. *, p≤0.05. N=3. Statistical significance was assessed by unpaired t-test with two-tailed calculations. Data is presented as mean ± SEM. (R-V) CT scanning and quantitative analysis of the muscle volume of control and Osr2 Cre ;Tgfbr1 fl/fl LVP at P0.5. A representative reconstructed control LVP from CT scanning is indicated in yellow (R and S) and an Osr2 Cre ;Tgfbr1 fl/fl reconstructed LVP is indicated in teal (T and U). **, p≤0.01. N=4. Statistical significance was assessed by unpaired t-test with two-tailed calculations. Data is presented as mean ± SEM. (W) A 300 μm coronal slice of the LVP region at E14 from Osr2 Cre ;Fgf18 fl/fl mouse for slice culture following bead implantation. Arrows point to BSA-or FGF18-treated bead. (X-Y) Immunofluorescence of MyoD (red) in the coronal section of LVP region from Osr2 Cre ;Fgf18 fl/fl mouse cultured for 3 days with BSA bead (X) and FGF18 bead (Y). White circles indicate the location of the BSA beads in X and FGF18 beads in Y. N=3. Fgf18 fl/fl or Fgf18 fl/+ littermates were used as controls for Osr2 Cre ;Fgf18 fl/fl mice. Scale bars in A, E, I, M, R, S, and W indicate 500 μm in in A-D, 100 μm in E-H, 50 μm in I-L, 400 μm in M-P, 400 μm in R and T, 100 μm in S and U, and 100 μm in W-Y, respectively.
The online version of this article includes the following source data and figure supplement(s) for figure 7: Source data 1. Source data for Figure 7Q.
Source data 2. Source data for Figure 7V.

Discussion
In the later stage of pharyngeal muscle development, myogenic progenitors migrate into discrete myogenic sites to form individual muscle anlagen anatomically resembling their adult counterparts (Noden and Francis-West, 2006;Sambasivan et al., 2011;Shih et al., 2008;Ziermann et al., 2018). How these fine-tuned individual muscles form the right morphology at the right location in order to perform their physiological functions remains largely unknown. Here, using the pharyngeal muscle LVP as a model, our study reveals that an important aspect of this mechanism is the establishment of a distinct pro-myogenic perimysial subdomain adjacent to the myogenic site. The cells in this domain have unique identities and regulatory mechanisms distinct from the rest of connective tissue cell populations within the palatal mesenchyme, enabling perimysial cells to uniquely provide various promyogenic signals and define specific myogenic sites to support region-specific myogenesis. Using unbiased screening with scRNAseq analysis combined with mouse genetic approaches, we identified TGF-β signaling as a predominant and specific regulator for perimysial cell fate determination during this stage, and perimysial transcriptional factor Creb5 assists TGF-β signaling to achieve functional specificity in supporting pharyngeal myogenesis via pro-myogenic signals such as Fgf18.
In early myogenesis, CNC cells have also been shown to induce myogenic differentiation by secreting both BMP and Wnt inhibitors to antagonize the dorsal neural tube-derived BMP and Wnt signaling molecules that repress craniofacial skeletal muscle formation (Tzahor et al., 2003). In the later developmental stages, when the CNC-derived palatal mesenchymal cells are more differentiated and patterned into distinguishable functional domains, we have the optimal opportunity to explore a specific regulatory mechanism of myogenesis. We found that during the LVP development, the late perimysial population (perimysial fibroblasts) with the most active TGF-β signaling is also distributed according to a muscle-specific pattern in both the fourth pharyngeal arch-derived LVP and first pharyngeal arch-derived masseter muscle; and this TGF-β signaling function is required for the formation of the muscle analgen of both the LVP and masseter (Han et al., 2014). Similarly, in limb muscle development, a group of pre-patterned Tcf4+ limb mesodermal cells also pre-determines the basic pattern of the muscles (Kardon et al., 2003). This indicates a potentially universal mechanism in which a pre-defined perimysial domain distinct from the rest of the connective tissues is required to establish specific myogenic sites that allow for proper muscle analge specification in late developmental stages.
By investigating perimysial-to-myogenic communication in this study, we have identified promyogenic signaling from the neighboring late perimysial cells (perimysial fibroblasts). It thus appears that these embryonic myogenic progenitors may also require 'embryonic niches' to support their contribution to muscle development at specific myogenic sites. In adults, muscle stem cells -satellite cells -reside in specific niches, and their ability to potentiate muscle repair and regeneration is also supported by signals from these niches (Andersen et al., 2013;Relaix et al., 2021). By comparing the 'adult muscle stem cell niche' with the 'embryonic niche', more molecular and cellular similarities can be further identified. For example, several perimysial-derived signaling pathways for the embryonic myogenic progenitors are also required for satellite cell fate regulation after birth (Bigas and Espinosa, 2016;Farin et al., 2016;Stantzou et al., 2017). In addition, while Hic1+ mesenchymal progenitors in the adult niche coordinate various aspects of skeletal muscle regeneration (Scott et al., 2019), they are also inferred to be regulators of the perimysial cells of the embryonic niche during development. These similarities suggest potentially conserved mechanisms between embryonic muscle development and adult muscle repair/regeneration, consistent with the clinical observation that patients with disrupted soft palatal muscle development also exhibit impaired differentiation and a reduced number of satellite cells (Carvajal Monroy et al., 2012). Thus, the mechanism we identified during development could also be potentially useful for promoting muscle growth after birth.
Craniofacial anomalies such as DiGeorge syndrome, as well as other syndromic and non-syndromic forms of cleft lip and palate, often affect pharyngeal muscle formation and lead to difficulties in eating, facial expression, speaking, and swallowing (Kelly et al., 2004;Kernahan et al., 1984;Carvajal Monroy et al., 2012;Scambler, 2000). While the identification of a detailed Pitx2-Tbx1-Msc-Tcf21 transcription factor regulatory network for the pharyngeal myogenic cells may influence the intrinsic regulatory network within the myogenic cells to promote pharyngeal muscle development (Buckingham and Rigby, 2014;Sambasivan et al., 2011), manipulation through paracrine signaling is a more likely mechanism. In this study, we have identified a variety of potential pro-myogenic signaling molecules from perimysial-to-myogenic cell interaction analysis. A 'chemokine/cytokine cocktail' may be needed to restore the muscle detects in our cleft soft palate animal model (Osr2 Cre ;Tgfbr1 fl/fl ). As a proof of concept, we found partial improvement of the muscle defects after Fgf18 restoration and thus confirmed the potential of restoring pharyngeal muscle defects with the combination of promyogenic factors we identified in the study. Moreover, some of these molecules, including Fgf18 and Dlk1, have also been suggested to be involved in the development and regeneration of various muscles throughout the body including in the limb, diaphragm, and tongue (Andersen et al., 2013;Ito et al., 2018;Mok et al., 2014;Yue et al., 2021), extending their potential application to defects in other muscles in the body.
In summary, our study has been able to decipher the molecular and cellular composition of the embryonic myogenic niche for the development of pharyngeal muscles and potentially other muscles in the body. Our work will contribute to a better understanding of the fine-tuned regulatory network of late-stage muscle morphogenesis and lead to the development of novel treatment strategies for regenerating/repairing muscle defects, particularly in infants with birth defects.

Animal studies
The Osr2 Cre (gift from Rulang Jiang, Cincinnati Children's Hospital; Chen et al., 2009), Rosa26 LSL-tdTomato (JAX#007905, Jackson Laboratory; Madisen et al., 2010), Tgbfr1 fl/fl (Alk5 fl/fl ) (Dudas et al., 2006;Larsson et al., 2001) mice. All mice were genotyped using genotyping primers as previously reported. The embryonic samples and newborn pups were collected and used for analysis without consideration of sex. For in vivo BrdU labelling, mice were injected intraperitoneally with 10 mg/ml BrdU in PBS, and 100 mg/kg BrdU was administered to each mouse. The mice were then euthanized for proliferation analysis 2 hr after BrdU administration. All animal handling followed federal regulation and was performed with the approval of the Institutional Animal Care and Use Committee (IACUC) at the University of Southern California documented under protocol numbers 9320 and 11765.

Tissue processing
Embryonic and newborn mouse heads were dissected and fixed in 10% formalin (HT501128, Milli-poreSigma) overnight at room temperature following decalcification in 10% EDTA depending on the stage. For paraffin sections, samples were processed in serially ascending concentrations of ethanol solution at room temperature followed by xylene and paraffin wax at 60 °C, then embedded in paraffin wax and sectioned at 8 μm using a microtome (RM2255, Leica). Deparaffinized sections were stained with Hematoxylin and Eosin (H&E). For cryosections, samples were dehydrated in 15% sucrose/PBS solution followed by 30% sucrose/50% Tissue-Tek OCT compound (4583, Sakura). Samples were embedded in the OCT compound and frozen on a block of dry ice. Embedded samples were then cryosectioned at 8 μmusing a cryostat (CM13050S, Leica).

Immunostaining
Sections were antigen-retrieved for 10 min in preheated antigen unmasking solution (H-3300, Vector Laboratories) followed by 10 min of incubation with 1% Triton X in PBS (T8787, Sigma Aldrich

Single-cell RNAseq
Soft palate tissue (dissected from the posterior third of the palatal region) from E14.5 control and Osr2 Cre ;Tgfbr1 fl/fl embryos was digested and dissociated into single-cell suspension with TrypLE Express enzyme (12605010, Thermo Fisher Scientific) at 37 °C with shaking at 600 rpm for 15 min using a Thermomixer (2231000269, Thermo Fisher Scientific). Single-cell suspension was loaded into the 10x Genomics Chromium system and prepared for single-cell library construction using the 10x Genomics Chromium single-cell 3' v2 reagent kit (PN-120267, 10x Genomics) according to the manufacturer's protocol. Sequencing was performed on Novaseq 6000 (Illumina). Library quality control, sequence alignment, and read counts were analyzed using Cell Ranger 4.0.0. Two control samples were combined using the function cellranger aggr so that the combined control sample had comparable cell numbers to the sample from Osr2 Cre ;Tgfbr1 fl/fl mice. One combined control and one Osr2 Cre ;Tgfbr1 fl/fl sample were used for analysis. For each sample, raw read counts from every single cell were analyzed to identify cell clusters and variably expressed genes in each cluster using the Seurat R package (Hao et al., 2021) as previously described (Han et al., 2021). Seurat was also used to combine the E14.5 control and Osr2 Cre ;Tgfbr1 fl/fl embryos to perform integration analysis as previously described (Han et al., 2021). RunPCA and RunUMAP visualization were used for downstream analysis and visualization. The integrated Seurat object combining E13.5-E15.5 control soft palates generated in our previous study (Han et al., 2021) was used to analyze gene regulatory network inference and ligand-receptor interactions. Gene regulatory network inference was performed using the R package SCENIC (Aibar et al., 2017). Transcription factors for each cell population in the palate were identified using GENIE3 and compiled into regulons, then subjected to cis-regulatory motif analysis. Regulon activity was then scored using AUCell. CellChat (Jin et al., 2021) was used to identify the potential ligand-receptor interactions between cell populations in the soft palate. Pre-processing functions (identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData) and core functions (computeCommunProb, computeCommunProbPathway, and aggregateNet with standard parameters) were applied along with other functions (netVisual_circle, netVisual_bubble, netA-nalysis_signalingRole_heatmap netAnalysis_dot, and netAnalysis_river) to determine the senders and receivers.

CRISPRi plasmid generation
gRNAs targeting close to the Smad2/3 binding sites (+150 bp to -150 bp) were identified using IDT guide RNA design tool. A gRNA with the sequence TTTC AGTT GAGT CACC CCAC was selected due to its targeting the proximity (47 bp downstream, 33098286-33098306) to the binding site, high on-target score, and low off-target probability. This gRNA sequence was then cloned to pCas-Guide-Puro-CRISPRi plasmid (GE100083, Origene Technologies) by Origene Technologies to generate a Fgf18 CRISPRi Plasmid. pCas-Guide-Puro-CRISPRi-Scramble plasmid (GE100084, Origene Technologies) was used as control.

C2C12 cell culture
C2C12 mouse myoblast cell line (CRL-1772, ATCC) was obtained from ATCC and cultured in growth medium composed of DMEM (11965092, Thermo Fisher Scientific)/10% Fetal bovine serum (12662029, Thermo Fisher Scientific)/1% Penicillin-Streptomycin (15140148, Thermo Fisher Scientific) and maintained at low density to avoid differentiation. After 1-2 passages, C2C12 cells were seeded at 3000 cells/cm 2 and cultured for 1 day in growth medium prior to the treatment. For proliferation assay, 500 ng/ml FGF18 (100-28, Peprotech) or an equal volume of BSA (RB02, R&D Systems) were added to C2C12 cells in growth medium for 1-3 days. For BrdU labelling, 10 µM BrdU was added to the medium for 2 hr. Cells were fixed in 4% paraformaldehyde (PFA) for 15 min at RT.

Imaging and image analysis
The immunofluorescent and RNAScope in situ hybridization fluorescent signals were captured using a Leica DMI3000 B research microscope. The brightfield images were captured using Keyence BZ-X710. Quantitative image analysis was performed by Image J. The same thresholds were set up for measuring the same group of samples. Adjacent sections in the LVP regions were measured and average values across sections were used for analysis. 'Analyze Particles' and 'Measure' functions were used to measure the number of cells and percentage of area of staining out of the whole palate area or whole image area.

MicroCT analysis
MicroCT scans were performed at the University of Southern California Molecular Imaging Center using a SCANCO µCT50 device. The embryonic heads were scanned with the X-ray source at 70 kVp and 114 µA to generate images at a resolution of 10 μm. Morphometric analysis was performed using the AVIZO 9.1 software package (Visualization Sciences Group). Four biological replicates were performed.

Statistical analysis
Sample sizes were extrapolated from sample sizes of previously published studies. No randomization or blinding were performed. No samples were excluded from the analysis. For qPCR analysis, 3 technical replicates were used for measuring each sample. GraphPad Prism was used for statistical analysis. All bar graphs display mean ± SEM. Un-paired two-tailed t-test was applied to assess statistical significance between two groups of samples. One-way ANOVA was used for analysis of more than two groups of samples. The chosen level of significance for all statistical tests in this study was p≤0.05.

Ethics
All animal handling followed federal regulation and was performed with the approval of the Institutional Animal Care and Use Committee (IACUC) at the University of Southern California documented under protocol numbers 9320 and 11765. Additional files

Supplementary files • MDAR checklist
Data availability The E14.5 control and Osr2 Cre ;Tgfbr1 fl/fl soft palatal single-cell RNAseq data generated in this study have been deposited in Gene Expression Omnibus (GEO) under accession code GSE203035. The single-cell RNAseq data for E13.5-E15.5 soft palate (Han et al., 2021) used for analysis in this study have been deposited in GEO under accession code GSE155928.
The following dataset was generated: