Single cell transcriptomics and developmental trajectories of murine cranial neural crest cell fate determination and cell cycle progression

Cranial neural crest (NC) cells migrate long distances to populate the future craniofacial regions and give rise to various tissues, including facial cartilage, bones, connective tissues, and cranial nerves. However, the mechanism that drives the fate determination of cranial NC cells remains unclear. Using single-cell RNA sequencing combined genetic fate mapping, we reconstructed developmental trajectories of cranial NC cells, and traced their differentiation in mouse embryos. We identified four major cranial NC cell lineages at different status: pre-epithelial-mesenchymal transition, early migration, NC-derived mesenchymal cells, and neural lineage cells from embryonic days 9.5 to 12.5. During migration, the first cell fate determination separates cranial sensory ganglia, the second generates mesenchymal progenitors, and the third separates other neural lineage cells. We then focused on the early facial prominences that appear to be built by undifferentiated, fast-dividing NC cells that possess similar transcriptomic landscapes, which could be the drive for the facial developmental robustness. The post-migratory cranial NC cells exit the cell cycle around embryonic day 11.5 after facial shaping is completed and initiates further fate determination and differentiation processes. Our results demonstrate the transcriptomic landscapes during dynamic cell fate determination and cell cycle progression of cranial NC lineage cells and also suggest that the transcriptomic regulation of the balance between proliferation and differentiation of the post-migratory cranial NC cells can be a key for building up unique facial structures in vertebrates.


Introduction 37
The craniofacial complex is one of the most diversified anatomical parts among vertebrates. 38 Many studies have been focused on external factors like feeding behaviors that drive the 39 selection of species-specific patterns of the craniofacial complex (Langenbach and van Eijden,

A single-cell graph of NC developmental process in the murine embryos 91
We combined lineage tracing and scRNA-seq to address the spatiotemporal dynamics and 92 investigate the cell fate decisions involved in NC cell differentiation after they migrate to their 93 target regions (Fig. 1A). To genetically label NC cells and their progenies, we crossed Wnt1-Cre 94 mice with the Rosa26-tdTomato/eGFP (Rosa26-mT/mG) reporter mouse line. The Wnt1-Cre-95 mediated recombination has been shown to occur in cranial and cardiac neural crest cells, 96 midbrain, and developing neural tube (Tavares and Clouthier, 2015). Therefore, it can label 97 these tissues and their progeny cells with green fluorescence. Since the cranial NC cells migrate 98 to their target region around E8.5 and the generation of the facial prominences is completed by 99 E12.5 (Ji et al., 2020), we collected eGFP positive cells from E9.5, E10.5, E11.5, and E12.5 100 embryos using fluorescence-activated cell sorting (Fig. 1B, C). Therefore, our data covered 101 most of the circuital stages for forming, fusion, and merging craniofacial prominences (Ji et al., 102 2019). Then we performed 10X single-cell RNA sequencing using live eGFP positive cells. The 103 sequencing detected roughly 3,382, 3,380, 3,012, and 2,465 genes per cell from E9.5, E10.5, 104 E11.5, and E12.5 embryos, respectively (Fig. S1). 105 In addition to the dorsal spinal cord, where the premigratory NC cells are located, the 106 Wnt1-Cre is also expressed at the midlines of the midbrain, the caudal diencephalon, and the Hbb-y, Hba-x, , which are small clusters containing only eGFP negative cells, were also 116 removed from further analysis. 117 After removing the non-NC cells from our dataset, a total of 19,640 single-cell 118 transcriptomes passed quality control for further analysis, including 3,620, 4,284, 4,302, and 119 7,434 cells from E9.5, E10.5, E11.5, and E12.5 embryos, respectively. Using Uniform Manifold 120 Approximation and Projection (UMAP), we analyzed the transcriptional heterogeneity of NC-121 derived cells. When the data from all four time points were analyzed together, we found that the 122 NC cells clustered into four major cell populations. 123 124

Identification of the major NC-derived cell types 125
We identified the major NC-derived cell types based on the expression of known marker genes. 126 The major NC-derived cell populations were identified as pre-EMT, which express Wnt1 and 127  Table S1). Notably, a small cluster representing the cardiac 131 neural crest cells that migrate to the outflow tract were identified based on their expression of 132 Tbx20 and Acta2 (Singh et al., 2005). These results argue that our data captured a majority of 133 known NC-derived cell types. 134 Previous studies showed that both pre-EMT and early migratory NC cells maintained 135 multipotency in mice (Baggiolini et al., 2015). One important characteristic of multipotent stem 136 cells is their ability to self-renew. Therefore, we analyzed the cell cycle states across the major 137 cell populations using the Seurat scRNA-seq analysis pipeline (Fig. 2C, Fig. S3). Indeed, 64.2% 138 of the pre-EMT and 86.7% of the early migratory neural crest cells are in the S or G2/M phase of the cell cycle, suggesting they are fast-dividing cells consistent with the hypothesis that they 140 are multipotent progenitor cells. In contrast, only 9.3% of neural linage cells are at the S or 141 G2/M phase, indicating that the neural linage is no longer actively dividing. Interestingly, our 142 results showed that 80.5% of the mesenchymal cells are at the S or G2/M phase, indicating that 143 most mesenchymal cells are also proliferating, the frequency of which was dynamic across the 144 sampled time points (Fig. 2D). The frequency of pre-EMT and early migratory NC cells declined 145 over time, while mesenchymal and neural lineage cells increased dramatically after E11.5. 146 These results suggest that our data represent the developmental process that NC progenitor 147 cells lose multipotency and differentiate into various fates. 148 149

Reconstruction of developmental trajectories of NC cells 150
One of the benefits of scRNA-seq is that it allows mapping dynamic differentiation process by 151 densely sampling cells at different stages, in our case E9.5, E10.5, E11.5, and E12.5. 152 Therefore, these NC-derived cells sampled at different time points can be used to create a 153 continually NC developmental trajectory. To reconstruct the progression that cells undergo 154 during their differentiation from pre-EMT NC progenitor cells to fate determined neural or 155 mesenchymal cells, we performed pseudotime analysis on all of the eGFP positive cells (Fig.  156 3A, B). The resulting lineage tree demonstrates the transcriptional changes associated with cell 157 fate splits. After emigrating from the dorsal neural tube, the fate determination that separate 158 neural and mesenchymal cell lineages happen twice (Fig. 3C). The first differentiation separates

Characterization of the mesenchymal cells at facial prominences 166
Once committed to a mesenchymal fate, cranial NC cells generate most of the bone and 167 cartilage in the craniofacial regions. However, it remains unclear how post-migratory NC-derived 168 mesenchymal cells are induced to differentiate into the specific craniofacial skeletal structures 169 with the correct size and shape. To reveal the mechanisms that drive the fate determination of 170 craniofacial mesenchymal cells, we re-clustered the 9,875 NC-derived mesenchymal cells 171 collected from E9.5, E10.5, E11.5, and E12.5 embryos (Fig. 4A, B, Table S2). Using cluster-172 specific marker genes and wholemount in situ hybridization (WISH), we identified fifteen clusters 173 associated with NC cells that migrate into the frontonasal prominences and the maxillary 174 prominence of the first pharyngeal arch (Fig. 4C-E). For example, the aristaless-like homeobox 175 Alx genes, Alx1, Alx3, and Alx4, are widely expressed in the midfacial complex at E10.5, and 176 cells that highly express these genes are enriched in cluster m1 (mesenchymal 1), suggesting 177 these clusters are from the midfacial prominences. By E11.5, Alx1 and Alx4 continue to be 178 widely expressed in all three paired midfacial primordia, the lateral nasal prominence (LNP), the 179 medial nasal prominence (MNP), and the maxillary prominence (MxP), suggesting that clusters 180 m1, 3, 7, and 11 represent the mesenchymal cells at midfacial prominences. The expression of 181 Alx3 is limited to MNP, indicating that clusters m3 and 11 are the MNP. In contrast, Barx1 is 182 widely expressed in the mandibular prominence of the first pharyngeal arch, the second, third, 183 fourth, and sixth pharyngeal arches at E10.5, suggesting clusters m2, 4, 5, 6, 8, 9, 10 are from 184 the pharyngeal regions. At E11.5, a small part of ventral MxP also begins to express Barx1, and 185 those cells are found in cluster m3. Together, using cluster-specific markers, we identified m1, 186 3, 7, and 11 as the mesenchymal cell populations of the midfacial prominences. Interestingly, 187 cells from early developmental stages (E9.5 and E10.5) are enriched in cluster m1, while cells 188 from late developmental stages (E11.5 and E12.5) are scattered in all four clusters, suggesting that the transcriptional heterogeneity for cells from E11.5 and E12.5 is more complicated 190 compared to E9.5 and E10.5 (Fig. 4B). 191

192
The differentiation of NC-derived craniofacial mesenchymal cells starts at a relatively late 193

developmental stage 194
We next performed a pseudotime analysis of cells from clusters m1, 3, 7, and 11 using Monocle 195 2 to reconstruct the development processes (Trapnell et al., 2014). Monocle 2 performed an 196 unsupervised analysis to order the cells and reconstructed a tree-like trajectory, beginning with 197 E9.5 and E10.5 cells and ending with E11.5 and E12.5 cells (Fig. 5A, B). Notably, cells in the 198 "root" branch are mostly from E9.5 and E10.5 embryos. By contrast, cells from E11.5 and E12.5 199 embryos are spread out in the other four states, indicating that the differentiation of craniofacial 200 mesenchymal cells is likely to initiate at a later developmental stage (Fig. 5C). based on the number of genes that were detected per cell (Dorrity, 2020). To test whether the 205 differentiation of craniofacial mesenchymal cells starts after E11.5, we measured the 206 developmental progression of craniofacial mesenchymal cells (Fig. 5F, G). The result shows 207 that the transcriptional complexity decreases in E11.5 and E12.5 embryos, leading to a higher 208 developmental progression score, while E9.5 and E10.5 embryos remain at early developmental 209 stages (Fig. 5F, G). In addition, we identified that more than 10,000 genes that exhibit temporal 210 expression patterns, which fall into two distinct gene groups (Fig. 5D, Table S3). One group of 211 genes are highly expressed in the early development stages (E9.5 and E10.5). The other group 212 of genes is highly expressed at the late development stages (E11.5 and E12.5) and includes genes like Msx1 and Col1a1 that are critical for craniofacial development (Fig. 5E) Therefore, although NC-derived mesenchymal cells populate at the craniofacial complex and 227 give rise to different prominences as early as E9.5, they remain at an undifferentiated, 228 homogeny stage. The differentiation of NC-derived mesenchymal cells into different bones in 229 the facial region does not initiate until E11.5, suggesting the fate of cranial NC cells might not be 230 intrinsically programmed but is acquired from the environment. 231 232

The initiation of the ossification is coupled to cell cycle progression 233
To reveal the mechanism that triggers the differentiation of cranial NC-derived mesenchymal 234 cells, we performed branch-dependent expression analysis at the first branch point (Fig. 6A, 235 Table S4). The results showed that one group (Group 1 in Fig. 6A) of genes is highly expressed 236 in the root branch, and the expression of these genes significantly decreased in both branches, 237 especially in cells at fate 1. Gene ontology (GO) enrichment analysis suggests that these genes 238 are involved in cell cycle regulation, chromosome segregation, and microtubule cytoskeleton 239 organization. In addition, many genes that inhibit the differentiation process are also 240 downregulated in both branches (Fig. 6B). For example, the inhibitor of DNA binding and cell 241 differentiation protein (Id1) has been shown to be expressed in embryonic and somatic stem 242 Cre dataset. Moreover, the expression level of all genes, the patterns of marker genes were 264 also similar between Wnt1-Cre and Sox10-Cre datasets (Fig. 7). Therefore, our analysis 265 revealed highly reproducible NC-derived cell populations associated with facial prominences, 266 which allows us to combine Wnt1-Cre and Sox10-Cre datasets for further analysis. 267 The facial mesenchymal cells from E11.5 embryos were populated into four clusters. . This is consistent with our results that, at E11.5, LNP and 278 MNP cells (m0, m1, and m2) were grouped into three clusters at different differential stages 279 while most of the MxP cells (m3) are still at an early developmental stage. 280 To further reveal the possible mechanisms that drive the differentiation of LNP and MNP 281 cells, we analyzed the differentially expressed gene expression between m2 and m0 (Table S6) 282 and m2 and m1 (Table S7). As a result, both Twist1 and Id2, known to be required for cell formation, such as Sox9, Col2a1, Itm2a, Igfbp5, and Col9a1, are expressed at a higher level in 286 m0 than in m2, suggesting cells in m0 are chondroprogenitors (Fig. 8F). Although cells in m1 do not have many highly expressed genes compared to m2, Igfbp5 and Itm2a were found to also 288 highly express in m1. Igfbp5 has been shown to stimulate bone cell growth (Miyakoshi et al., 289 2001). Studies in mice also showed that Itm2a was involved in osteogenic differentiation 290 (Tuckermann et al., 2000). These results indicate that m1 might represent a transition stage 291 between mesenchymal stem cells and chondroprogenitors. Moreover, we also found that 292 Wnt5a, Lef1, and Crabp1 were highly expressed in m2. In contrast, a WNT antagonist gene 293 Dkk2 was highly expressed in m0, suggesting Wnt and retinoic acid signaling might be essential 294 for maintaining the self-renewal of mesenchymal stem cells. showed that the NC-derived mesenchymal stem cell differentiation is also cell cycle-dependent. 329 The molecular mechanisms that underline this phenomenon still need further study. One 330 possible mechanism is the epigenetic landscape at the developmental genes that might change 331 in the G1 phase. In pluripotent stem cells, most developmental genes are H3K4 and H3K27 332 trimethylated near their transcription start sites to be silenced by a polycomb-dependent The embryos were collected at E10.5 and E11.5 and fixed in 4% PFA overnight at 4 ˚C. As   The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. 452  showing the developmental progression score for each cell type. Wilcoxon rank sum test was 734 used to compare cells from each cluster with m2. *** = p< e -14 , ** = p< e -3 (F, G) Volcano plots 735 showing differentially expressed genes between m2 and m0 (F) or m2 and m1 (G). Highly 736 expressed genes (avg_logFC > 0.5, P adj < 0.05) are colored in blue (in m2) or red (in m0 or 737 m1) dots, respectively. 738  Table S1. Top marker genes for major NC-derived cell types. 751 Table S2. Top marker genes for NC-derived mesenchymal clusters. 752 "p_val" and "p_val_adj" represent the probability for a gene to mark the cluster. avg_logFC, log2 753 of the ratio of the average expression level of a gene for all cells within the cluster compare to 754 all the other cells. "pct.1" and "pct.2" are the percent of cells in or out the cluster that express 755 the gene. 756 Table S3. Top gene lists in groups 1 and 2 that exhibit temporal expression patterns, 757 related to figure 5. 758