Research ArticleNEUROSCIENCE

Spatiotemporal cellular movement and fate decisions during first pharyngeal arch morphogenesis

See allHide authors and affiliations

Science Advances  16 Dec 2020:
Vol. 6, no. 51, eabb0119
DOI: 10.1126/sciadv.abb0119


Cranial neural crest (CNC) cells contribute to different cell types during embryonic development. It is unknown whether postmigratory CNC cells undergo dynamic cellular movement and how the process of cell fate decision occurs within the first pharyngeal arch (FPA). Our investigations demonstrate notable heterogeneity within the CNC cells, refine the patterning domains, and identify progenitor cells within the FPA. These progenitor cells undergo fate bifurcation that separates them into common progenitors and mesenchymal cells, which are characterized by Cdk1 and Spry2/Notch2 expression, respectively. The common progenitors undergo further bifurcations to restrict them into osteogenic/odontogenic and chondrogenic/fibroblast lineages. Disruption of a patterning domain leads to specific mandible and tooth defects, validating the binary cell fate restriction process. Different from the compartment model of mandibular morphogenesis, our data redefine heterogeneous cellular domains within the FPA, reveal dynamic cellular movement in time, and describe a sequential series of binary cell fate decision-making process.


The heterogeneous structures that form the face together serve crucial physiological and sociological functions in human life. From birth, faces are essential to individuals’ identity and serve as powerful indicators of our emotions and health status (1). The multiple distinctive cellular origins of the various facial tissues and sensory organs are vital to their anatomy and function, yet poorly understood (2). Numerous congenital craniofacial abnormalities affect the form and function of the face, such as cleft lip/palate and micrognathia (3). A comprehensive understanding of facial morphogenesis at the cellular and molecular levels is fundamental to the explanation and, ultimately, prevention of these developmental defects, which have serious physiological and psychosocial consequences and can substantially compromise the quality of life.

During facial morphogenesis, cranial neural crest (CNC) cells contribute to a variety of tissue types, including craniofacial skeletal structures, chondrocytes, odontoblasts, dental pulp cells, connective tissue, melanocytes, neurons, and glia of the peripheral nervous system (4, 5). Different sets of transcriptional networks specify the early (induction and migration), intermediate (intra-arch patterning), and late (terminal differentiation) stages of CNC cell development. Specifically, a hierarchical gene network composed of sequentially expressed transcriptional factors including Brn3c, Lhx5, Sox8, Tfap2b, and Ets1 has been identified to play an important role in regulating early CNC induction and migration (6). During the intermediate stage of CNC development, Soldatov and colleagues (7) recently performed single-cell RNA sequencing analysis and demonstrated the sequential binary fate restrictions of CNC cells during epithelial-to-mesenchymal transition and early migration. Coexpression and competition of genes drive the alternative fate programs toward sensory neuron, autonomic neuron, and mesenchyme (7). Shortly after CNC migration into the first pharyngeal arch, patterning domains are established within the mandibular process to set the primary identities of CNC mesenchymal cells. A new set of specifiers, such as homeobox genes Alx1, Barx1, Lhx8, and Gsc, which have been proposed to be induced by epithelial signals, divide the mandible into four patterning domains along the proximal-distal and oral-aboral axes (810). Multiple signaling networks have been identified as responsible for patterning domain establishment. Through epithelial-mesenchymal interactions, signaling molecules such as endothelin1, Bmp4, Fgf8, and Shh induce the expression of domain specifiers in the mesenchyme (9, 11, 12). Recently, this patterning process was tested by single-cell RNA sequencing analysis based on gene expression profiling, and it was revealed that Bmp and Shh signaling are involved in regulating patterning along the oral-aboral axis (12). However, despite advances in identifying specifiers during early CNC cell differentiation, we have limited understanding of whether the postmigratory CNC cells continue to move within the first pharyngeal arch and whether their fates are solely restricted by the ectodermal signaling, as proposed in the compartment/unit condensation model, and subsequently determined through cellular differentiation during craniofacial development (13, 14). From a developmental perspective, it is fascinating to investigate the dynamic processes through which postmigratory CNC cells commit to different fates and subsequently give rise to a variety of structures within the craniofacial complex. The expression of master regulators of chondrogenesis, osteogenesis, and odontogenesis, such as Sox9, Runx2, and Dspp, respectively, in the mandibular primordium indicates the beginning of the third stage, the terminal differentiation of postmigratory CNC cells (15, 16). Furthermore, the link between the domains established within the first pharyngeal arch in the patterning stage and different tissue types formed afterward is also largely missing (17, 18). Last, the locations of progenitor cells for the differentiated cell types and their migratory paths by which they reach the final destinations also remain unknown.

In this study, we investigated cellular heterogeneity during early mandible development and refined and mapped patterning domains with distinct differentiation potential back into the first pharyngeal arch in fine detail. We also transcriptionally profiled single cells from critical developmental stages, closely tracked CNC cells within the first pharyngeal arch, and identified a series of bifurcating fate decisions as they formed different lineages. Our findings reveal that the most proximal group of postmigratory CNC cells undergoes dynamic intrapharyngeal arch movement, following routes from the proximal toward the distal, aboral, and oral domains, and contributes to multiple lineages within the first pharyngeal arch. Last, we used single-cell analysis to identify how impairment of the patterning process affects CNC cell fate decision and leads to specific mandibular dysmorphology, which provides novel insight into the cellular dynamics of mandibular morphogenesis in normal and abnormal development.


Refining the patterning domains of the embryonic day 10.5 mouse first pharyngeal arch

To understand the cellular heterogeneity of postmigratory CNC cells and their roles in early mandible development, we isolated the mandibular primordium from mouse embryos at embryonic day 10.5 (E10.5) and sequenced individual cell transcriptomes (Fig. 1A). Through unsupervised clustering and marker analysis, we identified eight different cell types within the E10.5 first pharyngeal arch, including CNC-derived mesenchymal cells (Dlx5+ and Prrx1+), oral epithelium (Fgf8+ and Vgll2+), pharyngeal cleft (Wnt6+ and Tfap2b+), pharyngeal endoderm (Shh+ and Nkx2-3+), mesoderm-derived cells (Tbx1+ and Myf5+), endothelial cells (Cdh5+ and Cldn5+), erythroid cells (Hba-a1+ and Hbb-bt+), and glial cells (Plp1+ and Sox10+) (Fig. 1B) (1923). RNAscope analysis of known marker genes in mandibular sections confirmed the clusters derived from different germ layers (fig. S1A).

Fig. 1 Patterning domains in the first pharyngeal arch.

(A) Scanning electron microscopy scan of mouse embryo at E10.5 shows the general scheme of the intrapharyngeal arch patterning domains (top); histology of the first pharyngeal arch at E10.5 shows the homogeneous cellular morphology (bottom). PA1, first pharyngeal w of E10.5 mandibular primordium identified 20 major clusters representing eight different cell types. (C to F) FeaturePlots of patterning genes representing different domains. Lhx6 is expressed in the oral domain, Gsc is expressed in the aboral domain, Nbl1 is expressed in the proximal domain, and Hand1 is expressed in the distal domain. (G to J) RNAscope analysis of patterning genes representing different domains. (K) DotPlot of signature genes of mesenchymal clusters. (L) Schematic of cluster mapping within the first pharyngeal arch. The dotted circle indicates the mesodermal core. For example, clusters 2, 12, 6, and 7 as shown with Hand1, Alx3, and Lix1 expression are in the distal domain of the first pharyngeal arch.

In the E10.5 samples, patterning domains along the proximal-distal and oral-aboral axes were well separated, which was confirmed using RNAscope in situ analysis of genes known to be differentially expressed in each domain, including Lhx6, Gsc, Nbl1, and Hand1 (Fig. 1, G to J). We noticed that the clusters of cells that formed the oral and aboral as well as the proximal and distal domains were located in regions that mimicked their in vivo anatomical locations (Fig. 1, C to F). We also reclustered CNC-derived mesenchymal cells with a lower principal components analysis (PCA) dimension factor to identify the major differences among these cells. We observed five clusters. Clusters 0 through 3 represented the four patterning domains along the proximal-distal and oral-aboral axes, suggesting a significant influence of the patterning process on CNC-derived mesenchymal cells (fig. S1, B and C). Cluster 4, which was located in the center, had no apparent patterning gene marker expression, revealing an unidentified cell population residing either scattered throughout the tissue or in the central region of the E10.5 first pharyngeal arch (fig. S1, B and C).

In E10.5 samples, unsupervised clustering predicted 20 cell clusters, with 13 of them representing the CNC-derived mesenchyme, suggesting considerable complexity and heterogeneity among the postmigratory CNC cell population within the first pharyngeal arch (Fig. 1B). To facilitate the understanding of this complex patterning process, we performed extensive in situ analysis and mapped the 13 cell clusters representing CNC-derived mesenchyme back into their in vivo locations (Fig. 1L and fig. S2). On the basis of established patterning domain markers for mice (Lhx6 and Etv4 for oral, Gsc for aboral, Hand1 and Alx1 for distal, and Barx1 for proximal), we fit our 13 newly identified clusters within the CNC-derived mesenchyme into these four patterning domains. Specifically, clusters 1, 3, 4, and 15 represent the proximal domain; clusters 5 and 6 represent the oral domain; clusters 0, 7, 8, 9, and 11 represent the aboral domain; and clusters 2 and 12 represent the distal domain. In the zebrafish first pharyngeal arch, previous studies have named three specific patterning domains: dorsal, intermediate, and ventral (24, 25). Using partition-based graph abstraction (PAGA) analysis and Louvain clustering, we also observed a clear separation among these regions based on marker gene expression (Pou3f3 for dorsal, Dlx3 for intermediate, and Hand1 for ventral) (fig. S1, E and F) (26). The traditional patterning domains in the first arch give us the fundamental positional identity of CNC cells. However, the 13 fine-grained patterning domains we identified will allow us to reveal more nuanced differences and investigate the distinct differentiation potentials and future contributions of these heterogeneous groups of postmigratory CNC cells at a much higher resolution.

Single-cell analysis of lineage specification of CNC cells within the mandibular primordium at E12.5

As early as E12.5, the establishment of Meckel’s cartilage (MC) and bone primordium in the mandible suggests that CNC-derived mesenchymal cells have already acquired their fates and begun the process of differentiating toward distinct lineages (27). Clustering of the transcriptomes of samples at E12.5 and E14.5 revealed the expansion of a variety of cell types, such as epidermal, CNC-derived mesenchymal, and endodermal cells, over developmental time (Fig. 2A and fig. S3). An array of terminally differentiated tissue types derived from CNC cells emerged during this period, such as osteoblasts (Runx2+ and Sp7+), chondroblasts (Sox9+ and Col2a1+), odontoblasts (Pax9+ and Barx1+), and tendon and ligament precursors (Scx+ and Tnmd+) (Fig. 2A and fig. S3).

Fig. 2 Lineage specification of CNC cells within the mandibular primordium at E12.5.

(A) Cell annotation on UMAP plot of E12.5 mandibular primordium indicates different cell types based on known differentiation marker genes expression. (B) Pseudotemporal analysis of mesenchymal population showing the multilineage differentiation process of postmigratory CNC cells within the first pharyngeal arch. (C) Selected marker gene expression of each cell lineage in green. (D) Pseudotime trajectory of chondrogenic lineage. Top left: Pseudotime trajectory colored by cell clusters; bottom left: pseudotime trajectory colored by timeline, and heatmap showing transcriptional dynamics during differentiation. Specific genes are listed on the right to show their expression dynamics. (E and H) FeaturePlot of genes in the intermediate stage of chondrogenic differentiation at E10.5. (F, G, I, and J) Immunofluorescence (F and G) and RNAscope (I and J) detection of indicated protein or mRNA transcripts for genes in the intermediate stage (cluster II) at E10.5 and E12.5, respectively. M, Meckel’s cartilage.

To better demonstrate the multilineage differentiation of CNC cells, we isolated mesenchymal cell clusters at E12.5 and performed pseudotemporal analysis using Monocle 2 (28). The pseudotemporal differentiation analysis predicted that the early mesenchymal progenitors generated distinct trajectories to osteocyte, chondrocyte, and fibroblast lineages (Fig. 2B). Each of these trajectories was identifiable by known markers corresponding to the respective cell types (Fig. 2C). Consistent with Monocle, PAGA analysis (26) also showed a similar differentiation trajectory of E12.5 mandibular mesenchyme (fig. S5). We further examined each lineage trajectory individually. On the basis of the dynamics of their differentiation patterns, the genes associated with each trajectory were divided into three clusters (numbered I, II, and III) (Fig. 2D). Genes in cluster I were strongly expressed in the early stages of differentiation process. Postmigratory CNC cells in the mandibular arch strongly expressed genes in cluster I, such as Twist2, Epha4, and Osr2. Genes in cluster III were highly expressed in the terminally differentiated cells, so these genes indicated the characteristics of each progenitor cell type, including chondrogenic cells (Sox9 and Col2a1), osteoblasts (Runx2 and Mef2c), and fibroblasts (Dlk1 and Shox2) (Fig. 2D and fig. S4, A and B). Genes in cluster II were highly expressed in the intermediate stage of differentiation. We identified several first-arch aboral/distal domain patterning genes expressed in the intermediate stage of chondrogenesis, such as Gsc and Maf (Fig. 2D). In addition, Igf1 was also highly expressed in the aboral domain within the first pharyngeal arch at E10.5. These data suggest that the cells in the aboral/distal domain of the first arch are primed for mandibular chondrogenic differentiation, which is consistent with previous findings (29).

We further analyzed the contribution of cells located in the aboral/distal domain (clusters 7 to 9 and 12) at E10.5 to the future mandibular development. At E10.5, Gsc was highly expressed in the aboral/distal domain of the first pharyngeal arch, whereas Foxp2 expression was found in two distinct domains: the aboral/distal domain where Gsc expression was seen, and the aboral/proximal domain (Fig. 2, F and I). At E12.5, both Gsc and Foxp2 were specifically expressed in the perichondrium of MC, suggesting a possible regulatory role of these genes in mandibular chondrogenesis (Fig. 2, G and J). Gsc regulation of chondrogenesis has been reported previously. In Gsc−/− mice, defects of MC and the middle ear structure, which is derived from MC, were observed (30). We further investigated the contribution of Foxp2+ cells to chondrogenesis during mandible development using Foxp2Cre;tdTomato mice. As expected, at E15.5, Foxp2-derived cells contributed to both the perichondrium and mature chondrocytes (fig. S6, A and B), suggesting a potential role of Foxp2 in chondrogenesis (31). Overall, the pseudotime analysis faithfully reflected the multilineage differentiation process of postmigratory CNC cells within the first pharyngeal arch and also offered new insights into the potential molecular mechanisms that control the commitment of each cell lineage.

Integrative analysis of mouse mandibular development reveals sequential fate decisions of postmigratory CNC cells

To investigate the complex developmental process of mandibular morphogenesis over time, we performed integrated analysis of single-cell RNA sequencing data from three different developmental stages (E10.5, E12.5, and E14.5) using Seurat 3 (32). We identified 17 clusters that could be manually identified with 17 major cell types through their expression of highly specific marker genes (Fig. 3A; erythroid, endothelial, and myeloid clusters were removed for visualization purposes). On the basis of the developmental stage at which each cluster emerges, we can see the developmental progress of different cell types that give rise to each cell lineage through terminal differentiation (Fig. 3, A and C). For example, in the epithelium, early cells from E10.5 in the central region branched in three different directions, representing the oral epithelium (Pdgfa+), skin epithelium (Cxcl14+), and dental epithelium (Pitx2+) at later developmental stages (Fig. 3, A and C). The differentiation trajectory of CNC-derived mesenchymal cells predicts a different developmental progression pattern from that of the epithelium. Mesenchymal cells at E10.5 were tightly clustered on one side with two streams of cells branching out toward the later-stage clusters on the other side, connecting to cells from E12.5 and E14.5, among which we identified the osteogenic/odontogenic lineage, chondrogenic lineage, dermal fibroblast, and perimysial fibroblast lineages (Fig. 3, A and C).

Fig. 3 Integrative analysis of mouse mandibular development reveals sequential fate decisions of postmigratory CNC cells.

(A) Seurat 3 UMAP plot integration of mandibular primordium at E10.5, E12.5, and E14.5. (B) Marker gene expression in the integration analysis showing terminally differentiated lineages. (C) Seurat 3 UMAP plot integration of mandibular primordium at E10.5, E12.5, and E14.5 based on the developmental stages. The red-boxed area is the epithelial compartment, and the blue-boxed area shows the initial fate bifurcation of the postmigratory CNCs. (D) Monocle 3 pseudotime trajectory of CNC-derived mesenchymal cells showing the sequential fate decisions of postmigratory CNC cells. Top: Pseudotime trajectory colored by stage; bottom left: pseudotime trajectory colored by timeline; bottom right: schematic drawing representing the hierarchy of differentiation trajectory.

The multilineage differentiation ability of CNC cells has been well characterized in previous investigations, showing that during their epithelial-mesenchymal transition and early migration, they lose multipotency following a structured pattern of lineage restrictions (7). However, the means by which postmigratory CNC mesenchymal cells make their fate decisions and whether their lineage commitment occurs simultaneously in a compartmentalized model or in a dynamic sequential manner remain unclear. We used Monocle 3 pseudotime analysis to computationally model the cell fate decisions of CNC-derived mesenchyme. Our results demonstrated that CNC-derived progenitors acquire their fates through an orderly sequential bifurcation process. Specifically, the progenitors first bifurcate into two different groups: the Notch2+ stromal cell lineage and the common progenitors of the odontogenic, osteogenic, chondrogenic, and fibroblast lineages. The second bifurcation separates the odonto-/osteoprogenitors from the chondrogenic and fibroblast lineages (Fig. 3D).

The integration analysis predicts that the first fate choice of postmigratory CNC cells happens at E10.5. The commitment to each fate is marked by expression of Spry2, a gene encoding an Fgf signaling antagonist, in one branch, and Cdk1, a cell cycling gene, in the other (Fig. 4A). We identified that some cells strongly express both marker genes, suggesting the common origin of these two lineages (Fig. 4B). Consistent with this, in situ analysis showed Spry2 is highly expressed in the oral domain of the first arch at E10.5, whereas strong Cdk1 expression was detected in the proximal region with partially colocalized expression of Spry2 (Fig. 4, D to F), suggesting a possible progenitor role of postmigratory CNC cells in the proximal region of the first pharyngeal arch.

Fig. 4 Fate decision of postmigratory CNC cells at E10.5.

(A) FeaturePlot of markers representing the first bifurcated lineages toward common progenitors (Cdk1+) and stromal cells (Spry2+). (B) Cdk1 and Spry2 show largely branch-specific expression patterns, but are coexpressed in some cells at E10.5. X and Y axes indicate the expression level of each gene in the cells. (C) Heatmap showing expression of branch-specific genes in cells assigned to either lineage. (D to F) RNAscope detection of indicated mRNA transcripts showing lineage bifurcation of the postmigratory CNC cells starting from the proximal region of the first pharyngeal arch toward the oral and aboral directions. Arrow points to Spry2 expression in the oral domain, and arrowhead points to the Cdk1 expression in the proximal region. Boxed area highlights the coexpression of Spry2 and Cdk1. (G) E10.5 clusters in Seurat 3 integration embedding showing different differentiation potential of cells within the first pharyngeal arch mesenchyme. (H) RNA velocity analysis projected onto the t-SNE plot. The arrows indicate the predicted future cell states. (I) Heatmap of whole-mount immunofluorescence of Ki67 expression in E10.5 embryo. Arrow points to the strong expression of Ki67 in the proximal region of the first pharyngeal arch. (J to L) Immunofluorescence of Ki67 expression in E10.5 embryo section. Red-boxed area shows the proximal region of the first pharyngeal arch, and yellow-boxed area shows the distal region of the first pharyngeal arch. Arrows point to Ki67-positive cells.

Dynamic intra-arch cellular movement contributes to postmigratory CNC cell fate decision

We integrated the clusters gleaned from the E10.5 samples into the combined E10.5 to E14.5 datasets and mapped each cluster in the E10.5 mesenchyme to the differentiation trajectories defined by Seurat 3 integration analysis. This approach allowed us to investigate the differentiation potential and heterogeneity of CNC cells in the first pharyngeal arch at a much higher resolution than the traditional patterning domains. To our surprise, each cluster in the CNC-derived mesenchyme showed a unique pattern, collectively revealing a well-defined hierarchy of postmigratory CNC cell movement and differentiation. Cells in cluster 3 in the proximal domain represented the most undifferentiated population in the first pharyngeal arch mesenchyme, which has also been demonstrated in zebrafish (33), while cells in clusters 5 and 11 formed bifurcated streams very close to cluster 3, suggesting that these cells were just acquiring their cell fate and committing to different lineages. Cells in cluster 5 in the oral domain were strongly associated with the Notch2+ lineage, while the cells in cluster 11 served as common progenitors that further separated into trajectories toward clusters 6 and 7, which later on contribute to osteogenic/odontogenic, chondrogenic, and fibroblast lineages (Fig. 4G). However, this prediction of the differentiation potential of E10.5 cells is derived from gene expression similarities; genuine lineage tracing will be necessary for further analysis.

Using RNA velocity analysis of the E10.5 CNC-derived mesenchyme, we also validated this postmigratory CNC maturation direction from proximal to other domains (34). The velocity field projected onto the t-stochastic neighbor embedding (t-SNE) plot showed strong directional flow originating from clusters 1, 3, and 4, representing cells in the proximal domain moving in two distinct directions: cluster 12/2, representing cells in the aboral/distal domain, and cluster 5, representing cells in the oral domain (Fig. 4H). RNA velocity analysis based on scVelo (35) predicted the same directional flow from the proximal domain to the aboral/distal domains and oral domain (fig. S8, A and B). Consistent with this, pseudotime analysis of E10.5 mesenchyme also predicted that cells in the proximal domain would give rise to two lineage trajectories: one toward the oral domain, and the other toward the aboral/distal domains (fig. S8, C and D). Taking these findings together, our data suggest that besides anatomical difference, a cellular hierarchy of differentiation exists within the first pharyngeal arch mesenchyme.

To validate this discovery in vivo, we first examined the proliferation activity of cells in the first pharyngeal arch at E10.5. The heatmap based on whole-mount immunofluorescence of Ki67 expression indicated that cells in the proximal domain had significantly increased proliferation relative to those in other domains, suggesting that the cells in this region may be less differentiated (Fig. 4, I to L). This observation was further confirmed by phospho-histone H3 staining and quantification analysis (fig. S7). We next identified that Gbx2 is specifically expressed in the proximal domain of the first pharyngeal arch at E10.5 (Fig. 5, A and B). To investigate the contributions of Gbx2+ cells in the proximal domain, we performed a cell lineage tracing study using Gbx2CreER;tdTomato mice. We first harvested Gbx2CreER;tdTomato mice samples 10 hours after tamoxifen induction, showing that only the cells in the proximal domain of the first pharyngeal arch were labeled (Fig. 5, C and D). Twenty-four hours after tamoxifen induction, Gbx2-derived cells were still largely residing in the proximal region of the first pharyngeal arch. However, these Gbx2+ cells started to split into two separate streams that moved in the oral and aboral directions around the mesoderm core, confirming the prediction made by the integration analysis (Fig. 5, E and F). After 4 days of tracing, Gbx2-derived cells contributed to all the major lineages in the mandible, including osteoblasts, chondroblasts, dental mesenchyme, and connective tissue (Fig. 5, I to S, and fig. S9). Furthermore, after 6 days of tracing, the majority of the mandible was labeled with tdTomato+ cells, with the proximal region labeled more abundantly than the distal region, suggesting a directional movement of cells in the first pharyngeal arch (Fig. 5, G and H).

Fig. 5 Dynamic intra-arch cellular movement contributes to postmigratory CNC cell fate decision.

(A) FeaturePlot of Gbx2 expression at E10.5. (B) RNAscope detection of Gbx2 expression at E10.5 in the proximal region of the first pharyngeal arch. (C to G) Whole-mount views of lineage tracing of Gbx2CreER;tdTomato mice at different time points, as indicated. (C, E, and G) Sagittal views. (D and F) Frontal views. Arrows point to tdTomato-positive cells. (H) Section of E15.5 mandible of Gbx2CreER;tdTomato mouse. (I to M) Immunostaining of Sox9 expression on section of Gbx2CreER;tdTomato at E12. Boxed area in (L) is enlarged in (M). Arrows in (M) point to chondroblasts that are derived from Gbx2+ cells. FITC, fluorescein isothiocyanate. (N to S) Immunostaining of Runx2 expression on section of Gbx2CreER;tdTomato at E13.5. Yellow and green boxes in (Q) are enlarged in (R) and (S), respectively. Arrows in (S) point to osteoblasts that are derived from Gbx2+ cells. (T to X) Mandibular arch ex vivo culture with DiI labeling on day 0 and day 2. Boxed areas in (T) to (W) are enlarged in (T′) to (W′), respectively. Arrows in (U′) point to migrated CNC cells. (X) Quantification of the DiI-positive area of the ex vivo culture after 2 days. Scale bars, 500 μm.

To better monitor the dynamic movement of the CNC-derived mesenchymal cells in the first pharyngeal arch, we used an independent method, DiI (a lipophilic stain), to trace cells in ex vivo first pharyngeal arch culture. After 2 days, labeled cells underwent vigorous movement from the proximal region in the distal direction, avoiding the mesodermal core in the center (Fig. 5, T to U′). In contrast, in the distal domain, the cells labeled with DiI mainly remained in the local region without much cellular movement throughout the culture period (Fig. 5, V to W′). However, the movement of the cells in the proximal region slowed down after 3 days of culture, and after 6 days of culture, DiI-labeled cells contributed to the perichondrium of the MC (fig. S11, A and C). To better track the movement of the labeled cells in the proximal domain, we also performed time-lapse imaging of DiI labeling of the proximal mandibular arch. To optimize culture conditions, pictures were manually taken at 4-hour intervals after DiI injection. In the video generated from these images, we visualized cellular movement from the proximal region toward the distal (ventral) direction (movie S1). Taking these findings together, our data suggest previously unappreciated movement of postmigratory CNC cells within the first pharyngeal arch from the proximal to the oral and aboral/distal domains, eventually differentiating into odontogenic, chondrogenic, osteogenic, and fibroblast lineages. This dynamic cellular movement likely occurs in parallel with local cell expansion during mandibular morphogenesis. Overall, our study highlights how postmigratory CNC cells undergo cell differentiation in vivo through orderly migration and sequential cell fate decisions.

Disruption of CNC cell patterning leads to impaired lineage commitment

Next, to test whether disruption of the early patterning domain would affect CNC cell fate commitment, we took advantage of a well-established mutant model, Wnt1Cre;Alk5fl/fl mice, in which the transforming growth factor–β (TGF-β) type I receptor is ablated in all the CNC cells. These mice have compromised canonical TGF-β signaling, which is known to disrupt first pharyngeal arch patterning (36). At E10.5, Alk5 mutant mice showed increased apoptosis in the proximal domain of the first pharyngeal arch, compromised progenitor population, and, eventually, defective mandible development. Using Seurat 3 to integrate single-cell sequencing data from control (4705 cells) and Wnt1Cre;Alk5fl/fl mice (4690 cells) at E12.5, we discovered that there was an increased number of cells in cluster 2 (Dlk1+) and a reduced number of cells in cluster 7 (Runx2+) in the Alk5 mutant sample (Fig. 6A). Violin plotting visualized the down-regulation of Runx2 and up-regulation of Dlk1 expression. Specifically, Dlk1 was ectopically expressed in cluster 7 in the Alk5 mutant, whereas in the control, the Dlk1 expression was minimal (Fig. 6C). Quantification of cell number percentage based on three biological replicates showed that there was a decrease in Runx2-expressing osteoblasts and an increased number of Dlk1-expressing fibroblasts in the Alk5 mutant mice (Fig. 6B). These data suggest a potential CNC cell fate switch from osteogenic lineage to fibroblast lineage. To test this hypothesis, we performed pseudotime analysis on both samples, and consistent with our observations, the osteogenic lineage trajectory was missing in the mutant, while the chondrogenic and fibroblast trajectories remained (Fig. 6C). To validate these data, we analyzed Runx2 and Dlk1 expression in vivo. In the control, Runx2 was expressed on the lateral side of MC where the bone primordium was forming. Dlk1 was highly expressed in dermal fibroblasts underneath the skin as well as perimysial fibroblasts. In the Alk5 mutant sample, Runx2 expression at E12.5 was markedly reduced and Dlk1 expression expanded into the bone-forming region, consistent with the computational analysis (Fig. 6, D to K). By combining these approaches, we were able to monitor how the CNC cell domains and their migratory contributions are regulated by TGF-β signaling in normal and abnormally developing models.

Fig. 6 Disruption of CNC cell patterning leads to impaired lineage commitment.

(A) UMAP visualization of canonical correlation analysis of control and Wnt1Cre;Alk5fl/fl mice mandible primordia. Arrows point to the clusters representing fibroblasts (2) and osteoblasts (7). (B) Quantification of cell number percentage of Runx2+ osteoblasts and Dlk1+ fibroblasts in control and Wnt1Cre;Alk5fl/fl mice mandible primordia at E12.5. (C) Violin plots of Runx2 and Dlk1 expression in each cluster. (D) Pseudotime trajectories of mesenchymal cells of control and Wnt1Cre;Alk5fl/fl mice mandible primordia. The colors represent the cell clusters in t-SNE embedding. (E to L) Immunofluorescence of Runx2 and Dlk1 expression in control and Wnt1Cre;Alk5fl/fl mice at E12.5. Boxed areas in (E), (G), (I), and (K) are enlarged in (F), (H), (J), and (L), respectively. * in (H) indicates the reduced expression of Runx2. Dotted line in (J) outlines the mandible bone primordium. Arrows in (L) point to the ectopic expression of Dlk1 in mandible bone primordium in Wnt1Cre;Alk5fl/fl mice. M, Meckel’s cartilage; T, tongue.


The heterogeneity of postmigratory CNC cells has been well established through extensive studies on first pharyngeal arch patterning (9, 37). The anatomical position of the postmigratory CNC cells within the first pharyngeal arch dictates their differentiation potential, possibly due to the presence of local epithelial inductive signals (8, 38). The general consensus has been that the cells in the first pharyngeal arch are compartmentalized, undergo unit expansion with very limited movement, and only contribute to the structures in their immediate surroundings in both fish and mammals (3, 25, 39). Our study revealed three crucial properties of the postmigratory CNC population: (i) diverse heterogeneity represented by multiple cell clusters within the first pharyngeal arch, (ii) subpopulations with different movement ability and trajectories, and (iii) variation in the cells’ differentiation potential as they move from proximal toward the distal region of the first pharyngeal arch. Computational analysis predicted that the cells in the proximal domain are the least differentiated compared with those in other domains and that these cells choose their fate through a cellular movement process of bifurcation toward either the oral domain or the aboral/distal domain. This differentiation hierarchy within the first arch mesenchyme has not been previously realized (40). The Notch2 and Spry2 expression–driven lineage in the oral domain does not contribute to the odontogenic lineage. Although odontogenic-associated genes such as Pax9 and Lhx6 have similar expression patterns to those of Notch2 and Spry2, they belong to different lineages, as shown in our integration analysis. These data suggest that the oral region contains a mixed mesenchymal cell population with different functions and fates during first pharyngeal arch morphogenesis. This observation also suggests that the fate decisions of postmigratory CNC cells may not entirely depend on local inductive signaling. In the later stages, these cells gradually lose their location identities and spread out through the whole mandible (fig. S12, A to D). We also identified that multiple chromatin remodeling genes, such as Cbx3, Kdm6b, and Kmt2d, have highly enriched expression in the Notch2/Spry2+ cell population. Cell populations with similar gene expression profiles have also been identified in other developing organs, such as the frontonasal process and nephron (41, 42), suggesting these cells could be either transient supporting cells during organogenesis or a precursor cell population that will later differentiate into different cell types.

The traditional compartment/developmental unit model based on gene expression divides the first pharyngeal arch into several anatomical domains that are fairly localized (13, 14). These compartments have been suggested to be responsible for the local development of each region of the future mandible. Our cell proliferation analysis of the mouse E10.5 first pharyngeal arch also supports the notion that the local expansion of different domains may serve as a critical driving force guiding mandibular morphogenesis. However, this model has certain limitations in explaining the tremendous and dynamic shape changes during mandibular morphogenesis. Studies have shown that CNC cells may continue to move in the peripheral region of the pharyngeal arch, remaining mostly undifferentiated. These cells will interact with paraxial mesoderm-derived, ectodermal, or endodermal cells and differentiate during mandibular arch morphogenesis (17, 43). We show here that Gbx2-labeled cells in the very restricted proximal region retain some migration ability and undergo dynamic yet directional movement toward the distal region, contributing to multiple differentiated cell types. The proximal region of the mandibular primordium has higher mobility and proliferation activity than the distal region and may serve as a driving force to guide mandibular morphogenesis, possibly due to cells in other domains being more committed to their lineages, while cells in the proximal domain continued to maintain undifferentiated status, providing for the growth of the mandible. In the later stages, the condyle and angle of the mandible located in the proximal region also serve as growth centers that are responsible for mandibular bone elongation in postnatal growth (44). The critical role of cellular movement in embryonic morphogenesis and other biological processes has long been established, such as in convergent extension during axis elongation and in epithelial-mesenchymal transition requiring the loss of cell-cell contact (45). The mobility of the cells in the proximal region of the first pharyngeal arch is also associated with higher proliferation activity, suggesting that this population plays a precursor role. Consistent with this, osteogenesis initiates from the distal region of the mandibular primordium (46). Further studies are necessary to determine the precise movements of neural crest cells from proximal to oral, aboral, or distal domains within the first pharyngeal arch. Uncovering this cellular hierarchy of differentiation potential and cell mobility builds on the compartment model and further improves our understanding of the organized postmigratory CNC cell fate restrictions and the dynamic mandibular morphogenesis process.

Pseudotime analysis of E12.5 mandible primordium has revealed dynamic changes in gene expression during the differentiation of each cell type. Notably, genes in cluster II (Fig. 2D), which are highly expressed in the intermediate stage of each differentiation process, may represent candidate genes that control the differentiation of each progenitor cell type. For example, in the chondrogenesis trajectory, genes in the insulin growth factor (IGF) signaling pathway such as Igf1 and Igfbp5 are enriched in the intermediate stage, suggesting the involvement of IGF signaling in the initiation of chondrogenic differentiation (Fig. 2D). This is consistent with the finding that Igf1 induces chondrogenesis in mesenchymal stromal cell culture (47).

Ablation of TGF-β type I receptor in CNC cells leads to partial loss of the precursor population in the proximal domain in the first pharyngeal arch at E10.5 due to increased apoptosis (34). Consistent with this, Alk5 mutant mice show severe disruptions of both tooth initiation and mandibular bone formation. Significantly, although apoptosis only affects the very proximal region of the first arch, marked bone loss is found in both the proximal and aboral regions of the mandible of Alk5 mutant mice, suggesting the potential cellular contribution of cells from the proximal domain to the aboral region (36). We also observed a fate change of postmigratory CNC cells from osteogenic lineage to fibroblast lineage in Alk5 mutant mandibles. Our integration analysis suggested that, after the initial bifurcation, a second binary fate decision separates the common progenitors into osteo-/odontogenic and fibroblast/chondrogenic lineages. At E12.5, MC is still formed in the Alk5 mutant, whereas bone and tooth formation are both defective, suggesting that loss of TGF-β signaling mainly affects the second binary fate decision, not the first.

Our single-cell RNA sequencing analysis of the developing mandible provides a resource for studies of facial morphogenesis and various craniofacial defects associated with micrognathia. The dynamic movement and differentiation hierarchy of postmigratory CNC cells in the first pharyngeal arch may also be shared by other craniofacial structures that are derived from CNC cells, and highlight the dynamic nature and heterogeneity of cells in building the craniofacial complex.


Animal studies

The Gbx2-CreER knock-in (JAX#022135, the Jackson laboratory) (48), tdTomato conditional reporter (JAX#007905, the Jackson laboratory) (49), Wnt1-Cre transgenic line, and Alk5 floxed alleles (36, 50) mouse lines have all been described previously. The Institutional Animal Care and Use Committee at the University of Southern California reviewed and approved all animal work performed in this study. All work adhered to institutional guidelines. Timed matings were set up to recover embryos at the appropriate stage (E9.5 to E16.5). The sexes of the embryos were not determined.

Tamoxifen administration

Tamoxifen (T5648, Sigma-Aldrich) was dissolved in corn oil (C8267, Sigma-Aldrich) at a concentration of 20 mg/ml and injected intraperitoneally at a dosage of 0.8 mg/10 g body weight.

RNAscope in situ hybridization

Embryos harvested at predetermined stages were dissected and fixed in fresh 4% paraformaldehyde overnight at 4°C. Samples were passed through a sucrose series, embedded in OCT (optimal cutting temperature) compound (Tissue-Tek), and sectioned on a cryostat at 8 μm. Staining was performed using either an RNAscope 2.5 HD Reagent Kit-RED assay (322350, Advanced Cell Diagnostics) or an RNAscope Fluorescent Multiplex Detection Reagent kit (320851, Advanced Cell Diagnostics) according to the manufacturer’s instructions. Fgf8 (313411), Tbx1 (481911), Dlx5 (478151), Lhx6 (422791), Gsc (504821), Nbl1 (454541), Hand1 (429651), Alx1 (403161), Ebf1 (433411), Gbx2 (314351), Maf (412951), Osr2 (517481), Pitx1 (587451), Lix1 (552041-C2), Cgnl1 (587441), Spry2 (425061), and Cdk1 (476081-C2) probes were designed and synthesized by Advanced Cell Diagnostics.


Sections were immersed in preheated antigen unmasking solution (Vector, H-3300) in an Electric steam cooker (976l, Cell Marque, Sigma-Aldrich) for 10 min, followed by cooling at room temperature for 30 min. Sections were then incubated with blocking reagent (PerkinElmer, FP1012) for 1 hour and then primary antibody overnight at 4°C. After three washes in phosphate-buffered saline, sections were incubated with Alexa-conjugated secondary antibody (Invitrogen). Sections were counterstained with 4′,6-diamidino-2-phenylindole (DAPI; D9542, Sigma-Aldrich). Images were captured using a fluorescence microscope (Leica DMI 3000B) with filter settings for DAPI/FITC (fluorescein isothiocyanate)/TRITC (tetramethyl rhodamine isothiocyanate).

Immunostaining was carried out using the following antibodies: Runx2 (1:100, Cell Signaling, #12556), Dlk1 (1:100, Abcam, ab16667), Foxp2 (1:50, Abcam, ab9361), Ki67 (1:100, Abcam, ab28538), and Sox9 (1:1000, Abcam, 185230). Alexa Fluor 568 and Alexa Fluor 488 (1:200, Invitrogen) were used for detection.

First pharyngeal arch explant culture

The first pharyngeal arch of E10.5 embryos was dissected out and cultured in vitro in a Trowell culture system, which has been described before (51). In brief, the dissected first arch was cultured on the membrane filter supported by metal grids in BGJb media with 10% fetal bovine serum (FBS). The first arch was cultured up to 6 days, and the culture medium was changed every other day.

DiI labeling

The technique of DiI injection has been described before (52). DiI (Molecular Probes CellTracker CM-DiI, C-7000) was dissolved in 100% ethanol at a concentration of 1 μg/μl. The DiI stock was further diluted with BGJb medium containing 10% FBS at 1:1 ratio. Micropipettes were made using 1- to 5-μl Drummond Calibrated Pipettes (St. Louis, MI). The pipette tip was carefully placed at the desired position of the first arch, and DiI was injected using gentle air pressure. The position of DiI injection was validated by using whole mount or sections of the explant.


Data analysis of Fig. 5W (n = 7) was performed using GraphPad Prism 8, expressed as dot plot with mean point and error bars. The statistical significance was calculated using two-sided Student t tests.

Single-cell RNA sequencing

Isolation of cells and sequencing. Single-cell transcriptomes were obtained from the digestion of mandibles of embryos at different stages (E10.5, E12.5, and E14.5). In brief, a whole mandible was placed in trypLE enzyme (Gibco) on a thermomixer (Eppendorf) at 37°C for 15 to 30 min, depending on the stage of the embryo to release cells from the tissue. For mutant samples, fast genotyping was performed, while embryos were dissected out at designated time points and submerged in BGJb medium with 10% FBS on ice. Embryos were collected, and 1 to 2 mm of embryonic sac was harvested from each sample. Two hundred microliters of digestion mix (2 μl of protease kinase + 198 μl of tail buffer) was added to each sample, which was then digested in the ThermoMixer for 5 min (55°C + 2000 rpm). The digestion mix was then heat inactivated on an 85°C heat block for 15 min followed by regular polymerase chain reaction (PCR) amplification and gel running process. Then, tissue was dissociated into single cells according to the PCR results. Dissociation quality and cell viability were tested with a Countess II automated cell counter (Thermo Fisher Scientific) before proceeding to GEM (gel bead-in emulsion) generation. All our samples had more than 90% viable cells. For each sample, we targeted 10,000 cells. The number of actually sequenced cells was 3058 cells for E10.5, 4705 cells for E12.5, and 6788 cells for E14.5. On average, we had a sequencing depth of 80,000 read pairs per cell. At each stage, we used two biological replicates for sequencing analysis; for control and mutant single-cell RNA sequencing, three replicates were used for sequencing. Quality control, mapping, and count table assembly of the library were performed using the CellRanger pipeline version 3.0.2.

Identifying variable genes and dimensionality reduction. Raw read counts from the cells in each stage were analyzed using the Seurat R package (53). Standard Seurat object was generated, followed by filtering out low gene expression cells (cutoff, 200 genes per cell). Cell cycle regression normalization, ScaleData, and RunPCA functions (15 statistically significant principal components, cutoff of P = 0.05) were performed for the t-SNE dimension reduction and final visualization of the clustering.

Pseudotemporal reconstruction of lineages. We used the Monocle 2 algorithm (28) to reconstruct the multilineage differentiation pathways across the CNC-derived mesenchymal cells. The differential expressed gene list across different lineages was generated by M3DropDifferentialExpression function, and then standard reduceDimension and orderCells function were performed to estimate the alignment of cells along the trajectories. The estimation of the root of the trajectory was based on the cluster identities and marker gene analysis. We selected cells that were projected onto specific lineages on the basis of both the cell clusters they belong to and the “state” they were projected onto. The analysis of each trajectory was performed with standard protocol with default parameters. The unbranched paths were analyzed to identify genes that varied across pseudotime.

Integrative analysis of samples at different stages. Seurat 3 was used to combine the single-cell data of three stages and perform integration analysis. Standard workflow was used to first identify “anchors” among different datasets with the function FindIntegrationAnchors, and then Seurat objects were returned by passing these anchors to the IntegrateData function. Scaledata, PCA, and UMAP (uniform manifold approximation and projection) visualization were then performed for the downstream analysis and visualization. Monocle 3 pseudotime trajectory analysis was performed with Seurat 3 UMAP embedding to cell fate restriction hierarchy of postmigratory CNC cells across the different developmental stages. The identities of the cells at the first bifurcation point were extracted to analyze their Cdk1 and Spry2 expression level and generate the heatmap of their transcriptional profiles.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank B. Samuels for the critical comments on this manuscript. We also thank A. Smith and G. Crump for comments. Funding: This study was supported by NIH R01 DE012711, U01 DE024421, and U01 DE028729 to Y.C. Author contributions: Y.Y. and Y.C. designed the study, organized experimental work, and wrote the manuscript. Y.-h.E.L. performed computational analysis of single-cell data. T.-V.H., K.G., and A.W. generated whole-mount cell proliferation heatmap. X.H., J.F., J.H., and J.J. gave feedback on experimental aspects and implemented the data interpretation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All single-cell RNA-seq datasets are available through the FaceBase data repository (Record ID 1-DTK2, doi: 10.25550/1-DTK2; All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article