Integrated Transcriptome and Network Analysis Reveals Spatiotemporal Dynamics of Calvarial Suturogenesis.

Craniofacial abnormalities often involve sutures, the growth centers of the skull. To characterize the organization and processes governing their development, we profile the murine frontal suture, a model for sutural growth and fusion, at the tissue- and single-cell level on embryonic days (E)16.5 and E18.5. For the wild-type suture, bulk RNA sequencing (RNA-seq) analysis identifies mesenchyme-, osteogenic front-, and stage-enriched genes and biological processes, as well as alternative splicing events modifying the extracellular matrix. Single-cell RNA-seq analysis distinguishes multiple subpopulations, of which five define a mesenchyme-osteoblast differentiation trajectory and show variation along the anteroposterior axis. Similar analyses of in vivo mouse models of impaired frontal suturogenesis in Saethre-Chotzen and Apert syndromes, Twist1+/- and Fgfr2+/S252W, demonstrate distinct transcriptional changes involving angiogenesis and ribogenesis, respectively. Co-expression network analysis reveals gene expression modules from which we validate key driver genes regulating osteoblast differentiation. Our study provides a global approach to gain insights into suturogenesis.

INTRODUCTION

Craniofacial sutures are fibrous joints between adjacent skull bones. Within the developing suture, non-ossifying suture mesenchyme (SM) separates the edges of adjacent bones called the osteogenic fronts (OFs). Growth at the OFs of the calvarial sutures occurs through intramembranous ossification, by which mesenchymal osteoprogenitors proliferate and differentiate directly into osteoblasts (OBs). OBs secrete an extracellular collagen matrix called osteoid that is mineralized to form bone. During development, the bones remain separated to allow growth of the skull in coordination with the brain and to accommodate physical forces affecting the skull, even through adulthood (Morriss-Kay and Wilkie, 2005; Richtsmeier and Flaherty, 2013).

The frontal suture (FS), also known as the metopic suture of the forehead in humans, lies between the frontal bones at the anterior midline of the calvarium. Uniquely, the metopic suture undergoes fusion during the third to ninth months of life in humans, while other calvarial sutures do not fuse until adulthood (Manzanares et al., 1988). In mice during the first month of life, posterior fusion of the intramembranous frontal bones occurs by a secondary process of endochondral ossification, by which cartilage anlagen first form within the SM and are replaced by OBs that bridge the frontal bones. In contrast, the anterior frontal and other calvarial sutures never fuse (Sahar et al., 2005). Thus, the FS serves as a model for both suture formation and physiologic calvarial fusion.

Metopic suture defects are a significant source of human pathology such as craniosynostosis (CS), the premature fusion of the suture, and other craniofacial dysostoses. These conditions can require significant invasive surgical intervention in infants, particularly with CS, to prevent increased intracranial pressure leading to secondary neurological deficits (Twigg and Wilkie, 2015). Metopic CS is the second most common single-suture CS with a prevalence of about 1 per 5,000 live births, increasing in recent decades (Cornelissen et al., 2016), and 75% of cases are non-syndromic (Lajeunie et al., 1998). The genetic etiology of metopic CS is unknown for the majority of cases, but rare mutations have been found in ASXL1, COLEC11, ERF, FGFR1, FREM1, GLI3, IHH, LRP5, KRAS, MEGF8, MSX2, RUNX2, TGFBR1, and TGFBR2 (Heuzé et al., 2014; Twigg and Wilkie, 2015), implicating a wide variety of molecular mechanisms and cellular processes. Conversely, the metopic suture is pathologically wider in syndromes such as cleidocranial dysplasia, craniofrontonasal syndrome, and other frontonasal dysplasias (Hennekam et al., 2010). Genes mutated in these phenotypes include ALX3, ALX4, EFNB1, MSX2, and RUNX2.

Mice provide an excellent model for the study of suturogenesis, and numerous mutant mouse models reproduce the phenotypes of human mutations (Holmes, 2012; Lee et al., 2019b). Twist1 and Fgfr2 are important in regulating the balance between maintenance of SM and osteogenic differentiation. TWIST1 proteins inhibit or promote Fgfr2 expression in the SM or OFs, respectively, depending on their level of heterodimerization with other basic helix-loop-helix transcription factors or homodimerization (Connerney et al., 2006, 2008). Fibroblast growth factor (FGF) signaling promotes osteoprogenitor proliferation and differentiation in the OFs (Iseki et al., 1999). In Saethre-Chotzen syndrome, caused by loss-of-function TWIST1 mutations (el Ghouzzi et al., 1997; Howard et al., 1997), newborns can present with wide metopic sutures (Thompson et al., 1984; Young and Swift, 1985), and Twist1 haploinsufficiency causes a wide suture defect in neonatal mice (Ishii et al., 2003). This frontal defect persists in later development with delayed and less robust bone formation in the posterior frontal fusion (Hermann et al., 2012; Behr et al., 2011) and decreased repair of surgically induced frontal bone defects (Hermann et al., 2012). In Apert syndrome, caused by activating FGFR2 mutations (Park et al., 1995; Wilkie et al., 1995), newborns also present with wide metopic sutures that fuse after being filled in with ectopic bone (Faro et al., 2006), and a wide suture is found in an Apert syndrome mouse model (Wang et al., 2005).

Understanding FS development requires a detailed transcriptome map of the spatiotemporal organization of the suture. We used laser capture microdissection (LCM) and bulk RNA sequencing (RNA-seq) of the SM and OF regions of the FS at embryonic days (E)16.5 and E18.5 from wild-type (WT) mice to generate a comprehensive atlas of genes involved in normal suturogenesis. Distinct gene expression signatures between these regions identified functional specializations such as cell communication and signaling in the SM and proliferation and ossification in OFs. Differential gene splicing highlighted the importance of post-transcriptional regulation for modulating the composition of the extracellular matrix (ECM). Single-cell RNA-seq (scRNA-seq) of dissected sutures also at E16.5 and E18.5 identified mesenchymal and osteogenic cell subpopulations that were spatially arranged along a differentiation trajectory of osteogenesis and differed along the anteroposterior (AP) axis of the suture. We examined in vivo changes to the transcriptome and cell subpopulation structure in mutant FSs from Twist1 and Fgfr2 mice. Transcriptional changes affecting angiogenesis and ribogenesis distinguished the two mutants, respectively, while the cell subpopulation structure was not significantly altered. Co-expression network analysis of the SM and OFs further characterized the transcriptional organization of these regions and identified a mesenchymal gene expression module that included Twist1 and several key driver genes involved in OB differentiation.

RESULTS

Comprehensive RNA-Seq Defines Distinct Transcriptional Profiles of SM and OFs

To create a comprehensive atlas of gene expression within the FS, we performed bulk RNA-seq profiling of the SM and OFs of the FS from WT C57BL/6J mice. These regions were isolated by LCM at E16.5, when OFs are widely separated, and E18.5, when OFs are more closely opposed and sutures are more morphologically distinct (Figure 1A). We first characterized expression in the SM and OFs and found that across both stages, there were a combined 4,282 differentially expressed genes (DEGs) out of 12,947 detected genes (Figure 1B). Of these, 2,139 were more highly expressed in the SM (false discovery rate [FDR] ≤ 0.01; Figure 1C), and 2,141 were more highly expressed in the OF (FDR ≤ 0.01; Figure 1D) at one or both stages. Additionally, the expression of two genes (Usp13 and 4930546H06Rik) was higher in the SM at E16.5 but higher in the OF at E18.5. In both regions, approximately half of the DEGs were upregulated at both stages (Figures 1C and 1D). Gene Ontology (GO) analysis revealed significant biological processes (BPs) specific to each region. At both stages, SM gene expression was enriched for anatomical structure development, multicellular organismal processes, cell communication, and signaling (Figures 1C and S1A). At both stages, OF gene expression was enriched for ossification and cell division processes (Figures 1D and S1A), while sterol biosynthetic processes were prominent at E16.5 (Figure 1D).

Figure 1.

Regional and Temporal Gene Expression Differences between SM and OFs of the FS

(A) LCM of the FS at E18.5. Section planes through frontal bones (F) and suture (S) at anterior, middle, and posterior suture levels are indicated as black lines on an ALPL-stained mouse head. SM and OF regions before and after LCM are shown on cresyl violet-stained sections. Bottom schematic indicates position of unmineralized osteoid and mineralized bone. Scale bar: 200 μm.

(B) Hierarchical clustering of average expression changes for 4,282 genes (rows) with significantly increased (red) or decreased (blue) expression (FDR ≤ 0.01) in the OF compared to SM at E16.5 and E18.5 (columns). Color bar indicates the average log2 fold-change (FC).

(C) Overlaps (top) and GO enrichment (bottom) of SM DEGs at E16.5 and E18.5 (p % 0.05, dashed line).

(D) Overlaps (top) and GO enrichment (bottom) of OF DEGs at E16.5 and E18.5 (p ≤ 0.05, dashed line).

(E) Heatmap of average expression changes between regions and stages (left) and GO BP categories (right) for 12 SM or OF genes selected for ISH validation. Significant expression changes are asterisked. Associated GO BP categories are indicated in green. Gene names are at left.

(F) RNA ISH (blue) for the DEGs from (E) in the E18.5 suture, counterstained with nuclear fast red. Arrows indicate bone edges visualized by overlaid autofluorescence (yellow green). Bottom row shows bone autofluorescence (left, green) corresponding to mineralized bone stained with von Kossa (middle) as seen in merged images of adjacent sections (right). Scale bar: 200 μm.

(G) Hierarchical clustering of average expression changes for 2,480 genes with significantly increased (red) or decreased (blue) expression at E18.5 compared to E16.5 in the SM and OFs (FDR ≤ 0.01).

(H) Overlaps (top) and GO enrichment (bottom) of E16.5 DEGs by region (p ≤ 0.05, dashed line).

(I) Overlaps (top) and GO enrichment (bottom) of E18.5 DEGs by region (p ≤ 0.05, dashed line).

(J) Relative expression (2−ΔCT) by qRT-PCR for the indicated stage-specific DEGs by region. Data are shown as individual measurements from biological replicates with mean ± SD. *p % 0.05 (Student’s t test).

See also Figure S1 and Table S1.

Few SM-associated genes have been described to date, but known genes identified in our data included Fgf9, Fgf18, Il11ra1, Msx1, Msx2, and Thbs1 (Connerney et al., 2006; Kim et al., 1998; Nieminen et al., 2011; Ohbayashi et al., 2002; Opperman, 2000). Likewise, known OF-associated genes identified in our data included Bmp2, Dlx5, Fgfr1, Fgfr2, Fgfr3, Id1, Ihh, and Ptch1 (Holleville et al., 2003; Klopocki et al., 2011; Opperman, 2000). OFs were enriched for genes known to be expressed from early to late OB differentiation, including Alpl, Runx2, Sp7, Ibsp, Bglap, and Dmp1 (Table S1). Among the top DEGs in SM (FDR % 0.01, >5 fold-change) were many genes with no defined function in sutures, including those with specialized roles in muscle (Myh8) and myofibroblast (Acta2) cytoskeletal contractility, tendon (Tnmd), cytoskeletal interactions (Cpxm2), Wingless-related integration site (WNT) signaling (Sfrp2), and a marker of adult WNT-regulated stem cells (Lgr5) (Figures 1E and S1A). The most significant OF-associated genes included those encoding a cross-linking ECM collagen (Col22a1), a nuclear envelope lamin-interacting protein (Mlip), a bone morphogenetic protein 8a (Bmp8a), a sphingolipid metabolic protein required for cell growth (Sgms2), and proteins involved in WNT signaling (Dkk1, Tcf7) (Figures 1E and S1A). In each case, RNA in situ hybridization (ISH) at E18.5 confirmed their region-specific expression patterns (Figure 1F). The OF genes showed a similar expression along the AP axis of the suture, but some SM-enriched genes were variably expressed along this axis. Myh8 and Acta2 expression was detected only in the anterior SM, while Tnmd showed expression throughout the anterior SM but became restricted to stratified domains within the posterior suture (Figure S1C).

The SM-enriched genes suggested interesting SM functions. Smooth muscle actin encoded by Acta2 is a marker of myofibroblasts and is critical for their cytoskeletal contractile function. The expression of myofibroblast markers such as Col4a1, Fn1, Itgb5, and Postn (Stempien-Otero et al., 2016) was significantly enriched in SM at E16.5 and E18.5 (Table S1). Additionally, the expression of tendon markers such as Scx, Mkx, Col3a1, Col5a1, Col12a1, Col14a1, Fmod, Lum, and Dcn (Gaut and Duprez, 2016) was significantly enriched in SM at E16.5 and E18.5 (Table S1). The presence of these genes suggests a role for the SM as a mechanoresponsive connective tissue. Lgr5 is a potential stem cell marker and was expressed in scattered cells within the SM (Figure 1F). In the postnatal suture, Axin2, Ctsk, Gli1, and Prrx1 expression reporter constructs have been found to mark suture stem cells (Debnath et al., 2018; Maruyama et al., 2016; Wilk et al., 2017; Zhao et al., 2015). In our data, Prrx1 and another potential stem cell marker, Lrig1, were enriched in SM at E16.5 and E18.5 (Table S1). Axin2 was enriched in SM at E18.5, while Gli1 was enriched in OFs at both stages (Table S1). Suture stem cells have been identified only in the postnatal suture, and the timing of their formation is not known (Holmes et al., 2020a).

Transcriptional Profiles of SM and OFs Differ between Developmental Stages

We next characterized DEGs by stage and found that across both regions, there were a combined 2,480 DEGs (FDR ≤0.01; Figure 1G). Of these, 1,243 were more highly expressed at E16.5 (Figure 1H), and 1,234 were more highly expressed at E18.5 (Figure 1I). Additionally, the expression of three genes (Nell1, Rnft2, Slc24a4) was higher at E16.5 in the SM, but higher at E18.5 in the OFs. At both stages, approximately half of the DEGs were more highly expressed in SM, a quarter were more highly expressed in OF, and a quarter were more highly expressed in both regions (Figures 1H and 1I). GO BP analysis showed that at E16.5, both SM and OF gene expression was enriched for nucleosome assembly, cell cycle, and DNA replication (Figures 1H and S1B). At E18.5, SM gene expression was enriched for ECM organization and cell surface receptor signaling, while OF gene expression was enriched for regulation of transcription, inactivation of MAPK activity, and ossification (Figure 1I). At different stages, both regions showed altered expression of immune response genes. These transcriptional changes are consistent with a decrease in proliferation and an increase in maturation in both the SM and OFs from E16.5 to E18.5. We validated several DEGs in the major GO categories by qRTPCR, confirming that Capn6 and Igfbp3 were upregulated in SM and that Tnn was upregulated in OFs at E18.5 compared to E16.5 (Figure 1J).

Differential Splicing between SM and OFs Predominantly Impacts ECM Genes

We used an annotation-free approach to examine differential splicing (DS) between regions and stages of the FS, based on changes in intron usage of aligned RNA-seq reads. We identified 261 DS intron clusters in 238 genes (Figure 2A), of which 52% contained at least one previously unidentified exon. The most common type of DS event was exon skipping, followed by alternative 5′ exon usage presumably driven by alternative promoters, and alternative 5′ and 3′ splice site usage (Figure 2B). The numbers of DS genes identified between the SM and OFs, and between E16.5 and E18.5, were 136 and 146, respectively (Figure 2C). Of these, 44 were significantly enriched (FDR q <0.05) for mRNA splicing regulation and the ECM (Figure 2C; Table S2), suggesting that DS events mainly modulate the functions of splicing factors themselves and of ECM proteins.

Figure 2.

Differential Splicing in the FS

(A) Heatmap of the maximum percent spliced-in (PSI) changes for 261 intron clusters with significant (FDR ≤ 0.05) SM and OF regional and E16.5 and E18.5 temporal differential splicing (DS) identified by LeafCutter.

(B) Proportional breakdown of event types for regional and temporal DS intron clusters (A5SS, alternative 5′ [donor] splice site; A3SS, alternative 3′ [acceptor] splice site).

(C) Overlaps (top) and top GO enrichment categories (bottom; BP, biological process; MF, molecular function; CC, cellular component) for 238 genes with DS intron clusters between regions (orange) or stages (green).

(D) Overview of significant DS intron clusters in Col5a1 associated with increased expression of exon E64A and decrease of the mutually exclusive exon 64B in SM compared to OFs at E16.5 and E18.5 (TSPN, thrombospondin N-terminal-like domain; COLFI, fibrillar collagen C-terminal domain; Laminin_G, laminin G domain; Collagen, collagen triple helix repeat [20 copies] domain). Increased intron usage in SM or OF is indicated in purple and green, respectively. FDR-corrected p values are indicated for each cluster for E16.5 (turquoise) and E18.5 (salmon). Details of the DS event indicating the PSI for each intron in SM and OF and the changes in PSI (DPSI) between regions for both stages. The p values for DS are indicated above each intron cluster.

(E) Relative expression (2−ΔCT) by qRT-PCR confirming PSI changes at Col5a1 exons E64A and E64B (clusters 3107 and 3108) in SM and OF at E18.5. Data are shown as individual measurements from biological replicates with mean ± SD.

(F) Similar to (D), but for the Postn gene (FAS1, fasciclin-like domain; TGFBI/POSTN, transforming growth factor-beta-induced protein/periostin domain). The p values for DS are indicated above each intron cluster.

(G) Similar to (E), but confirming PSI changes at Postn exon E17 (cluster 2561; E17I, included; E17S, skipped).

See also Figure S2 and Table S2.

Many of the DS events had potential functional effects on protein isoforms, with 46.3% of events altering predicted protein domains. Notable examples included Col5a1 and Postn. The tendon marker Col5a1 encodes the alpha 1 chain of type V collagen, and the SM showed increased expression of isoforms containing exon 64A with a concomitant decrease in exon 64B at E16.5 and E18.5, potentially altering a C-terminal fibrillar collagen domain that may be involved in substrate binding (Mitchell et al., 2012) (Figure 2D). Postn encodes periostin, a secreted ECM protein that is a key regulator of cell behavior and ECM organization (Bonnet et al., 2016). Postn isoforms in the E16.5 and E18.5 SM showed increased inclusion of exon 17 that affects the TGFBI/POSTN (transforming growth factor-beta-induced protein/periostin) domain (Figure 2F). DS events (percent spliced-in) in each gene were validated at E18.5 by qRT-PCR (Figures 2E and 2G). Other examples of DS events include genes encoding proteins with known roles in mineralization and bone development such as Col6a3, Fn1, and Macf1 (Bentmann et al., 2010; Izu et al., 2016; Yin et al., 2018), as well as transcription regulation, such as Zfp207 (Figure S2). Taken together, our data suggest that alternative splicing plays a key role in the regional and temporal organization of ECM proteins in the developing FS.

Single-Cell Profiling Identifies Two SM- and Three OF-Related Subpopulations

While the sequencing depth of bulk samples provides a comprehensive regional representation of gene expression and transcript diversity, scRNA-seq can identify different cell types within the suture. We profiled single cells of dissected FSs in duplicates at E16.5 (1,450 cells) and E18.5 (5,182 cells) to a mean depth of 28,062 unique molecular identifiers (UMIs) per cell and a median of 4,517 genes per cell (22,240 genes in total). Unsupervised graph clustering of all scRNA-seq profiles identified 15 cell clusters (Figures 3A and S3A; Table S3). We preserved differences between proliferating and non-proliferating cells, as these are tightly linked to embryonic development and cell differentiation (Mayer et al., 2018). The identified clusters included various cell types of the hematopoietic lineage, such as dendritic cells (1.6%), mast cells (4.2%), monocytes (3.0%), macrophages (14.1%), and osteoclasts (1.4%). Capillary endothelial cells (3.8%) and pericytes (1.3%) were also present (Figure 3A). All cell types were found at both stages (Figure 3B) but at different relative proportions. The increase in macrophages at E18.5 (Figure 3B) was consistent with an increased expression of genes involved in immune processes at E18.5 in the bulk data (Figure S1B). Gene expression differences between stages within each cell type further contributed to the stage-specific DEG profiles we observed in the bulk data (Figures S3B–S3D). Therefore, changes in SM and OF gene expression during suture development result from a combination of changes in cell composition and expression changes within each cell type.

Figure 3.

scRNA-Seq Analysis of the FS at E18.5

(A) Uniform Manifold Approximation and Projection (UMAP) plot of cell-type clusters detected by unsupervised graph clustering of 6,632 cells from both replicates at E16.5 and E18.5. Related clusters of FS1–5, OB, HD, and DM cell types are outlined with a dashed line. Gene signatures used to identify cell types are in Table S3.

(B) Aggregate proportion of cell types in each stage. n = 1,450 cells for E16.5 and 5,182 cells for E18.5.

(C) Cell-type enrichment among genes with SM- or OF-associated expression at E16.5 or E18.5 in bulk RNA-seq data (Figure 1B). Significantly enriched cell types are asterisked (FDR ≤ 0.05).

(D) Significant GO BP categories of FS1–FS5, OB, HD, and DM expression signatures (p % 0.05, Fisher’s exact test).

(E) Average expression of the top five differentially expressed marker genes (boxed) ranked by FC (FDR ≤ 0.01, lnFC ≥ 0.1) for FS1–FS5, OB, HD, and DM. Genes selected for smFISH in (F) are colored according to their cluster membership. The complete list of cluster markers is in Table S3.

(F) smFISH of FS1–FS5, OB, HD, and DM for the indicated genes in the E18.5 suture. OFs and adjacent bone are indicated by dashed lines, where identical shapes are used for images from the same section. The merged image at lower left is a composite of channels from two adjacent sections. The merged image at lower right is of channels from a single section. Expression is summarized in the schematic at bottom, with FS5 omitted. Scale bar: 100 μm.

(G) Graph abstraction of the relationships among FS1–FS5, OB, HD, and DM. Cells are colored according to their cluster membership (left) or by pseudotime (right). A principal graph was fitted (solid and dashed lines) and rooted on the central FS1 subpopulation (circle) to determine the position of cells along an arbitrary pseudotime scale. See also Figure S3 and Table S3.

The majority of cells (72.1% at E16.5 and 70.2% at E18.5) were aggregated in a constellation of related clusters (Figure 3A) that included OBs, remnants of overlying hypodermis (HD) and underlying dura mater (DM) that adhered to the suture during dissection, and five FS clusters (FS1–FS5) that could not be matched to specific cell types. Several of these unidentified clusters were enriched for either SM- (FS1, FS2) or OF-enriched genes (FS3, FS4) identified in our bulk RNA-seq data (Figure 3C), suggesting that they comprised a spectrum of mesenchymal and osteoprogenitor cells in the developing SM and OFs. This was supported by GO analysis and comparisons of the FS1–FS5 subpopulations to the Mouse Cell Atlas (Han et al., 2018) (Table S3). FS1 and FS2 were enriched for gene signatures of mesenchymal stem cells (FS1) and/or neonatal calvaria stroma (FS1, FS2) (Han et al., 2018), as well as BPs related to ECM organization and regulation of extracellular signaling (FS1) or mesenchyme migration and tissue development (FS2) (Figure 3D). FS2 was notable for the expression of genes involved in cell contractility such as Acta2, Tagln, Actg2, and Myl9 and tendon markers such as Tnmd, Scx, and Col12a1 (Figures 3E and S3A; Table S3). The OF-associated subpopulations bore signatures of neonatal calvaria (pre)OB/bone/cartilage (FS3) and embryonic mesenchyme (FS4) (Han et al., 2018). The FS4 subpopulation was enriched for cell proliferation markers (Figure 3D; Table S3). Although a proliferation signature was also detected in other cell types, the majority of proliferating cells were within the FS4 subpopulation (Figures 3A and S3E).

Mesenchymal and OF Subpopulations Are Spatially Arranged along a Differentiation Trajectory

To further characterize the cell types of the FS, we investigated the spatial organization of the FS1–FS5, OB, HD, and DM cell populations at E18.5 using multiplexed single-molecule fluorescence ISH (smFISH). Probes were chosen to target at least one of the top five marker genes with high expression in each subpopulation and specific or significantly elevated expression levels relative to the other subpopulations (Figure 3E; Table S3). Visualization of HD (Clec3b) above and DM (Cxcl12) below the SM and frontal bones confirmed their identities as the hypodermis and the dura mater, respectively (Figure 3F). Ibsp expression, enriched in the OB subpopulation, was strongest at the bone edges and, together with the HD and DM locations, provided landmarks for the organization of the FS1–FS5 cell types. Among the SM-enriched clusters, Col8a1 expression, enriched in the FS1 subpopulation, was present throughout the SM from anterior to posterior (Figures 3F and S3F), while Acta2 expression, enriched in the FS2 subpopulation, occupied the central domain of anterior SM, as previously identified using our bulk RNA-seq data (Figures 3F and S3F). Acta2 expression was strongest in the upper half of the anterior SM (Figure 3F). Genes for both of the OF-enriched clusters, FS3 (Npnt) and FS4 (Top2a), showed the strongest expression at the edge of the frontal bones, marking the presumptive OFs (Figure 3F). Top2a expression was interspersed within the Npnt expression domain, consistent with previous observations of increased osteoprogenitor proliferation in the OFs (Iseki et al., 1999). Interestingly, FS3 also was detected in the anterior SM, adjacent to FS2, with stronger expression of Npnt in the lower half of the SM (Figure 3F). This FS3 mesenchymal expression was posterior to the Wormian bone typically found within the anterior FS of C57BL/6J mice (Figure 4A) (Zimmerman et al., 2019) and was not strongly expressed in more posterior SM (Figure S3F). Thus, the AP axis of the FS is characterized by the presence of regionally distinct cell subpopulations. Finally, Icam1 expression, enriched in the FS5 subpopulation, was detected diffusely throughout much of the suture (Figure 3F). FS5 was enriched for genes known to be induced as a result of tissue dissociation (Adam et al., 2017; van den Brink et al., 2017; Denisenko et al., 2019; O’Flanagan et al., 2019) (Figure S3G). As such, FS5 may not represent a specific endogenous cell subpopulation.

Figure 4.

Transcriptional Profiles of the FS in Twist1+/− and Fgfr2+/S252W Mouse Models of Suture Dysgenesis

(A) Morphology of E18.5 WT, Twist1, and Fgfr2 FSs at anterior, middle, and posterior locations (see Figure 1A for locations). Frontal bones, including OFs, are stained for ALPL (red) and nuclei (blue). Arrowheads indicate edge of OF. Wormian bones in anterior sections are asterisked. Scale bar: 100 μm.

(B) FS widths at E18.5 for WT, Twist1, and Fgfr2. Data are shown as mean ± SD. *p < 0.05 (Student’s t test).

(C) Hierarchical clustering of average expression changes of genes with significantly increased (red) or decreased (blue) expression in Twist1 or Fgfr2 FSs compared to WT in bulk RNA-seq data of SM and OF regions. Expression changes are shown as the log2 FC in expression between mutant and WT genotypes for each region and stage. Corresponding FDR-corrected p values are shown on the left as a heatmap, where shades of yellow represent significant changes.

(D) Cell-type enrichment of mutant expression signatures in bulk RNA-seq data. Significantly enriched cell types are asterisked (FDR ≤ 0.05).

(E) Venn diagrams showing the overlaps of significant DEGs between each mutant compared to WT in bulk and scRNA-seq datasets with corresponding hypergeometric p values indicated below each diagram.

(F) Similar to (E), but showing the overlaps of significant DEGs between bulk and scRNA-seq datasets for each mutant compared to WT.

(G) GO categories significantly enriched (p % 0.05; green) for DEGs between Twist1 and WT in bulk RNA-seq at E16.5 or E18.5 in SM or OF (top row; bulk), and per cell type at E18.5 (bottom panel of rows; single-cell). Categories are ordered from left to right by decreasing significance of enrichment among DEGs in the bulk dataset.

(H) Similar to (G) but for DEGs between Fgfr2 and WT.

(I) qRT-PCR for the indicated Twist1 and Fgfr2 DEGs in SM at E18.5. Data are shown as individual measurements from biological replicates with mean ± SD. *p % 0.05 (Student’s t test).

See also Figure S4 and Tables S1 and S4.

We next inferred an approximate graph abstraction to determine the relationships between the FS1–FS5 and OB clusters and delineate potential differentiation trajectories. The trajectory graph contained a main branch that showed a linear progression along pseudotime through the FS1, FS2, FS3, FS4, and OB subpopulations (Figure 3G), consistent with an osteogenic differentiation model from a mesenchymal to OB fate, and their spatial organization in the FS (Figure 3F). A side branch extended from FS1 to the HD, while FS5 and the DM were arranged together in another branch connected to FS2 (Figure 3G, dashed lines). It is unlikely that these additional branches represent alternative differentiation trajectories for SM cells, as the HD and DM are determined early in craniofacial development (Dasgupta and Jeong, 2019; Yoshida et al., 2008). Alternatively, it could reflect their common origin from the mesenchymal neural crest lineage that forms these tissues at this location (Yoshida et al., 2008).

Distinct Mesenchymal Expression Profiles Underlie Enlarged FSs of In Vivo Models of Suture Dysgenesis

Our definition of gene expression signatures and cell-type organization in the WT suture provides a basis for analyzing the transcriptome and identifying associated functional processes in mouse models of FS dysgenesis. The perinatal FS is abnormally wide in Saethre-Chotzen syndrome (Twist1) and Apert syndrome (Fgfr2) mouse models, similar to the human conditions. This widening was previously observed qualitatively by alizarin red staining of mineralized bone (Ishii et al., 2003; Wang et al., 2005). Therefore, it is unclear whether width defects result from disturbances to earlier processes of osteoprogenitor differentiation at the OFs or to later processes of OB mineralization. We measured the FS width in WT and mutant mice between the OFs visualized by staining for alkaline phosphatase (ALPL) enzymatic activity at E18.5. As the width of mutant sutures varied along the AP axis, measurements were taken within the anterior, middle, and posterior sutures (Figures 1A and 4A). The Twist1 FS was significantly wider compared to the WT, anteriorly and posteriorly (Figure 4B). This extends the phenotype of this mutant, as widening previously was described only in the posterior suture at P4 (Ishii et al., 2003). The Fgfr2 FS was significantly wider compared to the WT at all three locations (Figure 4B). These results show that earlier processes regulating suture differentiation are affected in both models.

To characterize the molecular causes of these defects, we first compared bulk gene expression in the SM and OFs isolated by LCM of both mouse models to WT mice. In total, 115 genes and 127 genes at E16.5 and/or E18.5 were differentially expressed in the Twist1 and Fgfr2 sutures compared to WT, respectively (FDR % 0.05). The vast majority of these changes (88.6% and 76.3%) occurred in the SM (Figures 4C, S4A, and S4B), and DEGs for both models were enriched for genes expressed in the SM-associated FS1 and FS2 subpopulations (Figure 4D), suggesting that it is the primary region impacted in each model. Despite these general similarities, only six genes (AW551984, C1qtnf3, Col2a1, Csmd3, Eya1, and Lhfpl2) had altered expression in both models (Figures 4E and S4; Table S1). Interestingly, homozygous mutants of Eya1 have severely impaired calvarial bone development (Xu et al., 1999), and Eya1 was downregulated approximately 50% at E18.5 in the SM (and OFs with borderline significance) of Twist1 mice and in the OFs of Fgfr2 mice compared to WT (Figure S4; Table S1). Each model was characterized by distinct changes in genes associated with the capillary endothelial cell population in Twist1 and with FS3, FS4, and HD cell populations in Fgfr2 mice (Figure 4D). We did not find significant differences in splicing between mutant and WT SM or OFs.

In the Twist1 suture, most SM changes occurred at E16.5 (Figures 4C and S4A; Table S1). GO analysis showed that DEGs were most significantly enriched for categories including multicellular organismal processes, vasculature development, anatomical structure development, and plasma membrane (Figure 4G). DEGs included Ctnna2, a catenin involved in cell-cell interactions mediated by adherens junctions, and the adult stem cell marker Lgr5, a top SM-associated gene in the WT suture (Figure S4A; Table S1), which were both upregulated in Twist1 SM compared to WT. The expression changes of Ctnna2 and Lgr5 were confirmed by qRT-PCR (Figure 4I).

In the Fgfr2 suture, the majority of SM changes occurred at E18.5 (Figures 4C and S4B; Table S1). Upregulated genes included a downstream target of FGF signaling, the transcription factor Etv5, which showed an almost 2-fold increase in Fgfr2 SM at E16.5 compared to WT (Table S1). ETV5 and TWIST1 proteins interact in the limb bud to repress Shh expression and regulate skeletal patterning (Zhang et al., 2010). Ptchd4, which encodes another Hedgehog (HH) pathway repressor, was upregulated in Fgfr2 SM compared to WT at E18.5 (Table S1). Downregulated genes in the Fgfr2 SM at E18.5, compared to WT, included 10 large and four small ribosomal protein genes. GO analysis showed that DEGs were most significantly enriched for categories related to the ribosome, amide biosynthesis, and extracellular region (Figure 4H). The latter category included Spon1, a known ECM regulator of endochondral ossification, whose expression change was confirmed by qRT-PCR (Figure 4I).

To assess changes at the cell-type level, we performed scRNA-seq profiling of a Twist1 and Fgfr2 FS at E18.5, when the mutant phenotype was first apparent (Figure 4A). All cell types identified in WT sutures (Figure 3A) were present at comparable levels in both mutant sutures (Figure S4C), indicating that there was no major reorganization of cell types in either mutant. Across all cell subpopulations, there were 1,239 DEGs in the Twist1 suture and 162 DEGs in the Fgfr2 suture compared to WT (Figure 4F; Table S4). The vast majority of these expression changes were less than 2-fold. There was significant overlap between DEGs identified in the bulk and scRNA-seq datasets for each mutant (Figure 4F). Between the Twist1 datasets, common genes included Twist1 and endothelial cell markers Emcn and Cdh5. Between the Fgfr2 datasets, common genes included the ribosomal protein genes Rpl35a, Rpl36al, and Rps26. The gene expression differences we observed specific to either bulk or scRNA-seq datasets likely reflect a difference in the ability of each approach to detect distinct aspects of transcriptional change. The former allows detection of more transcripts because of greater sequencing depth, but it may not detect expression differences specific to subpopulations, while the latter has greater resolution at the cell-type level but detects fewer transcripts. Despite these differences, there was strong concordance between the top BPs impacted by gene expression changes in the bulk and scRNA-seq datasets, with GO categories for cell-type DEGs overlapping with those of the bulk data across a range of cell types in both the Twist1 and Fgfr2 sutures (Figures 4G and 4H).

Gene Co-expression Network Analysis Identifies a Mesenchymal Module Involved in Suturogenesis

We next performed multiscale embedded gene co-expression network analysis (MEGENA) of our bulk WT and mutant RNA-seq data to identify modules of co-expressed genes that could signify specific suture processes and regulatory programs, especially those that may functionally differ between WT and mutants. MEGENA infers gene co-expression networks by embedding significant gene-gene correlations onto a topological sphere and identifies highly co-expressed gene modules at multiple compactness scales (Song and Zhang, 2015). Modules identified by MEGENA form a module hierarchy tree, where larger modules are progressively organized into branches of smaller submodules. Given the extent of gene expression differences between SM and OF (Figure 1), each region was analyzed separately, yielding an SM network with 286 modules and an OF network with 359 modules (Table S5). Finally, we leveraged our scRNA-seq data to superimpose cell-type-level information onto the network topology.

We focused on modules strongly and selectively enriched for suture-related processes (Figures 5A and S5A). In the OF network, 12 modules in six branches were related to developmental and/or ECM processes (Figure S5A; Table S5). Of these, the branch ending in the module OF107 was also highly enriched for genes specifically expressed in the OF (Figure S5B; Table S5) and the OB cluster (Figures S5C and S5D; Table S5), suggesting that they reflected key developmental processes related to reorganization of the ECM for osteogenesis. Module OF107 contained 103 genes, of which 15 were identified as hubs by key driver analysis (Zhang et al., 2013), which included important known regulators of mineralization such as Enpp6, Phex, and Spp1 and collagens such as Col1a1, Col1a2, and Col11a1 (Figure S5E).

Figure 5.

Co-expression Network Analysis of the SM

(A) Sunburst plot of the MEGENA co-expression network hierarchy of the SM based on aggregated bulk expression data from WT, Twist1, and Fgfr2 mice. The network hierarchy is represented as concentric rings where the center ring corresponds to the root modules, which further subdivide into increasingly smaller child modules in the outward rings. Modules are colored according to the most significantly enriched GO BP category (FDR % 0.05). The network hierarchy branch for module SM286 is outlined and indicated by an arrow.

(B) Similar to (A), but with modules colored according to bulk differential gene expression between region, stage, and/or genotype.

(C) Similar to (A), but with modules colored according to cell-type enrichment.

(D) Cell-type enrichment of the module SM286 genes. Significantly enriched cell types are asterisked (FDR ≤ 0.05).

(E) Co-expression network module SM286. Genes with significantly increased expression in SM or OF are shown in red or blue, respectively. Genes with altered expression in Twist1 and Fgfr2 (F) or tested in in vitro OB differentiation assays (G) are shown in bold.

(F) Heatmap of 24 genes in module SM286 with significant expression differences (FDR ≤ 0.05) between Twist1 or Fgfr2 and WT, in SM or OF regions at E16.5 or E18.5 (bulk), or per cell type at E18.5 (single-cell). Expression changes are shown as the log2 FC in expression between mutant and WT. Significant DEGs are asterisked.

(G) MC3T3-E1 cells overexpressing mCherry (control), Lgr5, or Cpxm2 were differentiated for 28 days in OB differentiation media and stained with alizarin red (left) or von Kossa (right) to identify mineralized modules.

(H) Expression of Lgr5 and Cpxm2 relative to Actb (×100) in MC3T3-E1 cells transfected with lentivirus expressing mCherry, Lgr5, or Cpxm2 (F) as indicated, determined by qRT-PCR. Lgr5 and Cpxm2 were overexpressed compared to mCherry control levels (*: Lgr5, p = 2.2 3 10 4; Cpxm2, p = 6.5 3 10 9; Student’s t test).

See also Figure S5 and Table S5.

The SM network showed a similar organization of basic cellular and suture-related processes, including modules enriched for DEGs in Twist1 and/or Fgfr2 mouse models (Figures 5A and 5B). Multiple lines of evidence pointed to a central role for the module SM286 in mesenchymal development. The parent modules of SM286 were enriched for GO categories of ECM processes (Figure 5A; Table S5). Comprising 205 genes (Figure 5E), SM286 was enriched for genes expressed in WT SM at E16.5 and E18.5 (fold enrichment [FE] = 4.39, adjusted p value = 1.94E 41; Figure 5B; Table S5) and FS1, FS2, and DM cell types (Figures 5C and 5D). The module also contained Twist1, along with 23 other genes that had regional or cell-type-associated expression changes in the Twist1 suture compared to WT (Figure 5F). Several of these have been previously implicated in bone mineralization or collagen organization (Dkk2, Thbs2) (Kyriakides et al., 1998; Li et al., 2005), OB differentiation (Satb2, Igfbp3, Dkk3) (Aslan et al., 2006; Britanova et al., 2006; Li et al., 2013), cartilage development (Itga10) (Bengtsson et al., 2005), and suture development (Antxr1) (Cullen et al., 2009). Key driver analysis identified 46 potential key regulators of the module SM286 (Figure 5E), which included the top SM-associated genes Cpxm2, Lgr5, and Sfrp2 (Figures 1E and 1F). Notably, Lgr5 was among the most significantly upregulated genes in Twist1 compared to WT SM (Figure S4A) and was also specifically induced in Twist1 FS1 cells (Figure 5F).

Finally, we sought to explore the role of the module SM286 in suture development. Two key driver genes in SM286, Igfbp3 and Dkk3, have been shown previously to inhibit osteogenic differentiation in vitro (Aslan et al., 2006; Li et al., 2013). Similarly, overexpression of SM- and FS1/FS2-associated key drivers Lgr5 and Cpxm2 strongly decreased mineralization compared to a homologous mCherry fluorescent reporter control in an in vitro OB differentiation and mineralization assay (Figures 5G and 5H). Altogether, these data provide strong evidence of an important role for module SM286 in controlling the differentiation process of mesenchymal cells in FS1 and FS2 subpopulations to OBs.

DISCUSSION

Understanding the cellular and molecular dynamics of craniofacial sutures in health and disease requires a comprehensive knowledge of their transcriptional organization. Using a complementary approach of combining bulk and scRNA-seq analyses, we revealed multiple aspects of FS organization and potential function. Significant expression of genes known to be involved in mechanoresponsive and adult stem cell processes and the organization of the ECM were identified. Transcriptional profiling of individual cells further revealed a complex suture environment comprising SM and OF cell subpopulations arranged along a differentiation axis from mesenchyme to OB. The presence of the FS2 subpopulation within the central region of the anterior SM indicates additional variation at the cellular level along the AP axis of the FS. Bulk and scRNA-seq analysis of in vivo models of FS dysgenesis showed that widening of the Twist1 and Fgfr2 SM is associated with changes in mesenchymal gene expression particular to each model, impacting angiogenesis and ribogenesis, respectively. Finally, integrated co-expression network analysis of WT and mutant data identified a mesenchymal gene expression module associated with FS1 and FS2 subpopulations that is implicated in controlling differentiation to OBs. The module contains Twist1, as well as the key driver genes Lgr5 and Cpxm2, for which we describe a role in regulating OB differentiation.

Mechanical tension positively regulates OB differentiation and postnatal suture development (Herring, 2008). A hypothesis of cranial development proposes that tension produced by the expanding brain promotes osteogenic precursor differentiation at the OFs (Moss, 1954; Pritchard et al., 1956; Richtsmeier and Flaherty, 2013). However, the role of the SM in this process is generally not considered. The expression of many genes typical of specialized contractile tissues such as myofibroblasts and tendon, evident in both the bulk and single-cell expression data, suggests that the embryonic frontal SM is a mechanoresponsive connective tissue. Acta2 expression is induced in the adult sagittal suture by tension (Takeshita et al., 2017), and its embryonic expression and that of other tension-responsive genes may be similarly induced. The SM or adjacent tissues also may respond to tension by the induction of extracellular osteogenic signals. Igf1 expression, enriched in HD, is potentially one such signal. Tension induces the expression of Igf1 and Igf1r, which is enriched in FS1, in the postnatal sagittal suture (Hirukawa et al., 2005; Takeshita et al., 2017). FS1 is also enriched for expression of Igfbp3 and Igfbp4, each of which can inhibit osteogenesis (Eguchi et al., 2018; Li et al., 2013; Zhang et al., 2003). Increased IGF1 pathway expression has been found in some cases of non-syndromic CS (Stamper et al., 2012), and primary OBs from such cases show increased contractility and reduced cell migration in vitro (Al-Rekabi et al., 2016). Signals induced by tension forces in the SM during calvarial growth then may play a role in defining the differentiation trajectory suggested by our scRNA-seq analysis. In the FS2 subpopulation that localizes to the anterior SM, there is enrichment of many of the myofibroblast and tendon marker genes; therefore, the detection of and response to mechanical forces may vary along the AP axis. This variation may contribute to the absence of fusion anteriorly, as compared to the occurrence of fusion in the posterior FS.

In the sutures, the precise origin of embryonic osteoprogenitor cells remains to be determined (Ishii et al., 2015). Postnatally, reporter constructs expressed by Axin2, Ctsk, Gli1, and Prrx1 have been found to mark suture stem cells within the SM (Debnath et al., 2018; Maruyama et al., 2016; Wilk et al., 2017; Zhao et al., 2015). The relationship between embryonic osteoprogenitors and postnatal stem cells in the suture, and the timing of postnatal stem cell formation, is unknown. We identified Lgr5 expression in isolated cells within the SM. Lgr5 is a marker of stem cells in a variety of tissues, including the intestine, stomach, hair follicles, mammary gland, and ovaries (Leung et al., 2018). Lgr5 therefore is a potential marker of an early suture stem cell. Lgr5− mice have no described skeletal phenotype (Morita et al., 2004). Other potential stem cell markers that we identified in the SM at both stages, and that were enriched in FS1, include Lrig1 and Prrx1.

Differential splicing of splicing factor genes themselves and of ECM genes between SM and OFs was a major finding of our bulk expression analysis. Spliceosomopathies caused by defects in splicing factors include craniofacial phenotypes, indicating the sensitivity of craniofacial development to these processes (Merkuri and Fish, 2019). The tendon marker Col5a1 encodes a fibrillar collagen that regulates ECM maturation and collagen fibril diameter in the skin, corneas, and tendons (Paladin et al., 2015). Mutations of Col5a1 cause Ehlers-Danlos syndrome, a collagenopathy with joint hypermobility, skin hyperelasticity, and tissue fragility with atrophic scarring (Malfait and De Paepe, 2014). Alternative exons 64A and 64B may modify the specificity of the collagen chain selectivity recognition domain encoded by exon 65 (Mitchell et al., 2012). Interestingly, exon 64B is used most frequently in actively proliferating tissue (Mitchell et al., 2012), and in our data, exon 64B was enriched in the more proliferative OFs compared to SM. Postn encodes an ECM protein involved in collagen fibrillogenesis and integrin-mediated signaling regulating cell adhesion and motility (Bonnet et al., 2016). The presence of the alternatively spliced exon 17 appears to promote a positive role for POSTN in angiogenesis (Nakama et al., 2016), an important process in intramembranous bone formation (Percival and Richtsmeier, 2013). The functional significance of differential splicing between SM and OF on the ECM environment, particularly for mechanical properties of the SM, modulation of extracellular signaling, cell adhesion, and motility, remains to be explored.

In the Twist1 and Fgfr2 FSs, scRNA-seq analysis showed that cell identities themselves apparently were not perturbed in either mutant line. Most bulk transcriptional changes occurred in the SM, suggesting that changes to SM processes affecting the timing or rate of osteoprogenitor differentiation resulted in the wider separation of the OFs. The largely non-overlapping identities of DEGs and associated GO terms in the bulk and scRNA-seq data between Twist1 and Fgfr2 suggests distinct etiologies of dysgenesis, although the expression of Eya1, which is required for proper intramembranous ossification in the calvaria, was downregulated in both models. Many of the genes misregulated in Twist1 play a role in vascular development. xTwist is required for normal embryonic vascular development in Xenopus (Rodrigues et al., 2008), TWIST1 promotes vascularization in tumor xenograft models (Mironchik et al., 2005), and Twist1 null mice die at E11 with vascular defects (Chen and Behringer, 1995). Decreased Twist1 expression therefore may adversely affect vascular development, impairing bone formation (Percival and Richtsmeier, 2013). Lgr5, a potential stem cell marker co-expressed with the Twist1 in the module SM286, was notably upregulated in Twist1 SM and FS1 cells, and overexpression of Lgr5 in MC3T3-E1 cells inhibited mineralization.

In the Fgfr2 SM at E18.5, 14 large or small ribosomal protein genes were downregulated compared to WT in the bulk RNA-seq data. Skeletal development is known to be sensitive to disturbances in ribosome biogenesis in the nucleolus, which is tightly regulated during OB differentiation and influenced by nucleolar localization of FGFR2 (Neben et al., 2017a, 2017b). Dominant FGFR2 mutations causing bent bone dysplasia syndrome increase the nucleolar occupancy of FGFR2, which upregulates rRNA transcription, resulting in increased osteoprogenitor proliferation and reduced OB differentiation (Neben et al., 2014). Potentially, reduced ribosome biogenesis may therefore restrict protein production, reducing advancement of the Fgfr2 OFs.

Co-expression network analysis of all bulk expression profiles identified a large number of modules and their associated key drivers that recapitulated known suture biology. Many of the key SM genes and expression changes in the mouse models, in particular Twist1, converged in the module SM286 in the mesenchymal network. This module was associated with ECM processes and included Twist1 as well as Frem1, mutations of which cause metopic CS in humans and mice (Vissers et al., 2011). Expression of two key regulators of this module, Lgr5 and Cpxm2, strongly decreased OB mineralization in MC3T3-E1 cells, suggesting an important role of this module early in the differentiation trajectory. Thus, co-expression network analysis was shown to be effective in identifying genes of functional significance in suture biology.

Altogether, our study provides detailed insights into the spatiotemporal organization of the FS and is a key resource for modeling human suturogenesis and dysostoses. The data resource presented here will facilitate a broad range of methodological approaches to study suture biology, from a more directed definition of spatial relationships between cell types and mutation of genes of interest to genome-wide assessments of gene networks that can be informed by other “-omics” datasets integrated within a multi-scale analysis of suture biology. These strategies will better inform approaches to understand the complexity of and to identify candidate genes for human skeletal dysplasias such as CS and other craniofacial dysostoses.

STAR★METHODS

RESOURCE AVAILABILITY

Lead Contact

Comments and/or requests for resources and/or reagents should be directed to and will be fulfilled by the Lead Contact, Harm van Bakel (harm.vanbakel@mssm.edu).

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

Bulk and single-cell RNA-seq data reported in this paper are deposited with FaceBase (https://www.facebase.org/) with accession numbers FaceBase: FB00000805, FaceBase: FB00000904, and FaceBase: FB00001013 (Brinkley et al., 2016; Holmes et al., 2020b).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mice

Mouse procedures were in compliance with animal welfare guidelines mandated by the Institutional Animal Care and Use Committee (IACUC) of the Icahn School of Medicine at Mount Sinai. Timed matings of WT male and female C57BL/6J mice (The Jackson Laboratory, 000664), male Twist1 (Chen and Behringer, 1995) and female C57BL/6J, or male Fgfr2 (Wang et al., 2005) and female homozygous EIIACre (Lakso et al., 1996) mice, were performed to obtain embryos at E16.5 and E18.5 as required for LCM or scRNA-seq analysis. Twist1, Fgfr2, and female EIIACre mice have been maintained on the C57BL/6J background for over 20 generations in the Jabs laboratory. For LCM, dissected embryos were decapitated, heads were placed in Tissue-Tek Optimal Cutting Temperature (OCT) compound (Sakura, 4583) in labeled square plastic molds (Thermo Fisher Scientific; Peel-A-Way Disposable Embedding Molds, S-22, 50–465-347) in a standard orientation, rapidly frozen in a methyl-2-butanol/dry ice bath, and stored at −80°C. To obtain E16.5 and E18.5 embryos for single-cell preparation, timed matings of WT male and female C57BL/6J mice were made. Genotyping was performed by PCR of DNA extracted from tail biopsies. Twist1 genotypes were detected with neomycinspecific primers and control Tcrd-specific primers as described previously (Duchon et al., 2011). Fgfr2 genotypes were detected with primers listed in Table S6. Sex genotypes were identified as described previously (Bean et al., 2001; Wang et al., 2010).

Cell lines

HEK293T/17 cells (ATCC, CRL-11268, female) and MC3T3-E1 Subclone 4 cells (ATCC, CRL-2593, sex undetermined) were maintained at 37°C with 5% CO2. These lines were purchased directly from the ATCC.

METHOD DETAILS

Laser capture microdissection

To exclude the confounding effects of having a mixed sex population in RNA-seq libraries, only male embryos were used for LCM. Embryonic heads were sectioned in the coronal plane using an Avantik QS11 cryostat at 12 μm onto PEN-Membrane 2.0 μm slides (Leica, 11505158). The start of the frontal suture was identified by inspection of sections using a microscope, and the entire suture was sectioned serially, with 8–12 sections per slide, on a total of approximately 12–18 slides. Each PEN slide was dried, placed on dry ice and transferred to −80°C for storage prior to LCM. As a “guide slide,” the first and subsequently every tenth section was collected on a Superfrost glass slide for ALPL and DAPI staining for orientation of the suture and selection of PEN slides.

For LCM, a minimum of one slide was selected from each of three locations along the AP axis of the frontal suture: anterior, within the midregion, and posterior of the eye. Slides were thawed briefly and then washed twice for 30 s each in 45 mL of 75% ethanol (Thermo Fisher Scientific, BP2818) with agitation in 50 mL Falcon tubes to remove OCT compound. Slides were then stained with 0.1% cresyl violet in 50% ethanol for 30 s, washed for 30 s in 45 mL of 75% ethanol, then dehydrated by passage for 30 s each through 45 mL of 95% ethanol, 100% ethanol, and xylene (Sigma-Aldrich, 214736) in 50 mL Falcon tubes. Slides were air-dried before proceeding to LCM.

LCM was performed using a Leica LMD6500 as described by the manufacturer. For each stage and genotype five individual embryos were used to produce five biological replicates of both the SM and frontal bone OFs. OF samples from the homologous frontal bones within each suture were combined as one sample. Tissue was collected in 50 μL of Arcturus Picopure extraction buffer (Thermo Fisher Scientific, KIT0204) in the cap of a 0.5 mL Eppendorf collection tube, immediately frozen on dry ice, and stored at −80°C. RNA was purified using the Arcturus Picopure kit using the Macrocap procedure, with on-column DNase digestion (QIAGEN, 79254), as described by the manufacturer, and eluted with 20 μL of elution buffer. Purified RNA was stored at −80°C. RNA concentration and purity (RNA Integrity Number, or RIN) was determined using an Agilent Bioanalyzer 2100 and the RNA Pico kit as described by the manufacturer (Motch Perrine et al., 2019; Wu et al., 2019).

Bulk RNA-seq library preparation and sequencing

Library preparation of RNA derived by LCM was performed by the Gene Expression Core Facility at the Cincinnati Children’s Hospital Medical Center (Cincinnati, OH). Briefly, an initial amplification step from 4–12 ng of RNA per sample was performed with the Ovation RNA-seq System v2 (NuGEN, 7102–32) to create double-stranded cDNA. Concentrations were measured using the Qubit dsDNA BR assay (Thermo Fisher Scientific, Q32850). The cDNA size distribution was determined by using a DNA 1000 Chip. To create libraries for each sample for sequencing with the Illumina protocol, the Nextera XT DNA Sample Preparation Kit (Illumina, 15032354 and 15032355) was used to create DNA library templates from the double-stranded cDNA as described by the manufacturer. Concentrations were measured using the Qubit dsDNA HS assay (Thermo Fisher Scientific, Q32851). Library size for each sample was measured using the Agilent HS DNA chip (5067–4626).

For each genotype a total of 20 libraries (five biological replicates of SM and OF at E16.5 and 18.5) were generated. WT and Twist1 library sets were synthesized and sequenced as a single pool of 40 samples. The Fgfr2;EIIACre library set was initially synthesized and sequenced as a separate pool of 20 samples. To allow for the correction of potential batch effects, RNA-seq libraries for three biological replicates for each region at each stage of WT and Fgfr2;EIIACre library sets were resynthe-sized together, multiplexed, and sequenced as a single pool of 24 samples in a single HiSeq run.

Illumina sequencing was conducted by the Genetic Resources Core Facility, Johns Hopkins Institute of Genetic Medicine (Baltimore, MD). DNA sequencing was performed on an Illumina HiSeq 2500 instrument using standard protocols for paired-end 100 bp sequencing. Average yield was ~15 Gb of raw sequencing data per lane, or ~300 million reads per lane. As per Illumina’s recommendation, 5% PhiX was added to each lane as a control.

Bulk RNA-seq analysis

Illumina data were processed through Illumina’s Real-Time Analysis (RTA) software generating base calls and corresponding base call quality scores. CIDRSeqSuite 7.1.0 was used to convert compressed bcl files into compressed fastq files.

After adaptor removal with cutadapt (Martin, 2011) and base quality trimming to remove 3′ read sequences if more than 20 bases with Q R 20 were present, paired-end reads were mapped to the murine mm10 reference genome using STAR (Dobin et al., 2013) and gene count summaries were generated using featureCounts (Liao et al., 2014). Raw fragment (i.e., paired-end read) counts were then combined into a numeric matrix, with genes in rows and experiments in columns, and used as input for differential gene expression analysis with the Bioconductor Limma package (Ritchie et al., 2015) after multiple filtering steps to remove low-expressed genes. First, gene counts were converted to FPKM (fragments per kb per million reads) using the RSEM package (Li and Dewey, 2011) with default settings in strand specific mode, and only genes with expression levels above 1 FPKM in at least 50% of samples were retained for further analysis. Additional filtering removed genes with less than 50 total reads across all samples or less than 200 nucleotides in length. Finally, normalization factors were computed on the filtered data matrix using the weighted trimmed mean of M-values (TMM) method, followed by voom (Law et al., 2014) mean-variance transformation in preparation for Limma linear modeling, where we fitted the next generation sequencing batch as a covariate with fixed effect. Data were fitted to a design matrix containing all sample groups and pairwise comparisons were performed between sample groups (i.e., suture region, developmental stage, and genotype). eBayes adjusted p values were corrected for multiple testing using the Benjamin-Hochberg (BH) method and used to select genes with significant expression differences (q ≤ 0.01).

Analysis of local splicing differences

Local splicing analysis was performed in an annotation-free approach using LeafCutter (Li et al., 2018). We used LeafCutter to first find clusters of variably spliced introns across all samples and then to identify differential splicing between suture regions and developmental time points jointly modeling intron clusters using the Dirichlet-Multinomial generalized linear model (GLM) (Li et al., 2018). We used LeafCutter to find intron clusters as follows: overlapping introns (i.e., exon junction-spanning reads) were clustered and filtered to keep intron clusters supported by at least 25 split reads across all samples, retaining introns of up to 100 kb and accounting for at least 1% of the total number of reads in the entire cluster. This yielded 10,412 clusters encompassing 27,024 introns in 5,853 genes that were used for further analysis. This intron count file was then used in the differential splicing (DS) analysis. DS intron clusters were identified in pairwise analyses between sample groups (OF compared to SM and E16.5 compared to E18.5). After discarding introns that were not supported by at least one read in 5 or more samples, clusters were analyzed for DS if at least 3 samples in each comparison group had an overall coverage of 20 or more reads. P values were corrected for multiple testing using the BH method and used to select clusters with significant splicing differences (FDR ≤ 0.05).

Gene ontology enrichment analyses

Gene ontology (GO) biological process (BP), molecular function (MF), and/or cellular component (CC) enrichment analyses were performed using the gProfileR R v0.6.4 package (Reimand et al., 2007). The background gene set was restricted to frontal suture-expressed genes (defined as genes with expression levels above 1 FPKM in at least 50% of samples). An ordered query was used, ranking genes by log 2 fold-change for differential gene expression analyses. P values were corrected using the g:SCS algorithm to account for multiple comparisons.

scRNA-seq library preparation

Frontal sutures, encompassing the SM, OFs, and narrow strips of adjacent bone, were dissected from C57BL/6J mice. Ectocranial and endocranial membranes were removed as much as possible while avoiding loss of suture tissue, and sutures of both sexes were combined for processing. For the WT analysis, two libraries at E16.5 and two libraries at E18.5 were created. The E16.5 libraries consisted of separate pools of 13 and 19 sutures and were each obtained on separate days. The E18.5 libraries consisted of separate pools of 7 and 10 sutures and were each obtained on separate days. Dissections were performed under a stereomicroscope in multiple changes of cold phosphate-buffered saline (PBS). For each of the four libraries suture strips were combined into a single pool and digested at 37°C in αMEM (GIBCO, 32571–036) with 0.2% collagenase type II (Worthington, LS004176), 0.2% dispase II (Sigma-Aldrich, 4942078001), and 1U/μl DNase (QIAGEN, 79254). E16.5 libraries were digested for up to 100 minutes with brief vortexing every 10 minutes. Cell suspensions were strained through a 40 μm filter (Falcon, 352340), pelleted at 400 g for 7 minutes, and resuspended in red blood cell lysis solution (Miltenyi Biotec, 130–094-183). Cells were washed and resuspended in PBS/1% BSA. E18.5 libraries were serially digested for 10–15 minutes with fresh changes of digest solution for up to 90 minutes. Successive fractions were pooled on ice with the addition of FBS to 2%. Cell suspensions were strained through a 40 μm filter and processed as described above, with the omission of red blood cell lysis. Libraries from each mutant line at E18.5 were prepared as described for E18.5 WT libraries, with an additional WT library prepared from the Fgfr2 line to increase statistical power. For the Twist1 line, 10 Twist1 sutures were pooled on the basis of the mutant calvarial phenotype, and genotypes were subsequently confirmed by PCR. For the Fgfr2 line, 4 WT littermates were pooled; 3 Fgfr2 littermates were pooled on the basis of the mutant calvarial phenotype; and their genotypes were subsequently confirmed by PCR except for one WT sample, which was found to be a partial recombinant without a mutant calvarial phenotype. scRNA-seq 3′ expression libraries were prepared on a Chromium instrument (10X Genomics, model GCG-SR-1) using the Chromium Single Cell Gene Expression kit (Version 3) by the Technology Development Facility at the Icahn School of Medicine at Mount Sinai. Sequencing was conducted by the Genetic Resources Core Facility, Johns Hopkins Institute of Genetic Medicine. Libraries were first run on an Illumina MiSeq at 26×8×98 and analyzed to confirm the number of captured cells and assess capture efficiency prior to sequencing. The initial four WT libraries were then sequenced as a single pool on an Illumina NovaSeq S1 using standard Illumina primers, where read 1 was the UMI, read 2 was the library index, and read 3 was the transcript. The mutant and WT littermate libraries were similarly sequenced as a second pool.

scRNA-seq analysis

Demultiplexing of cellular barcodes, obtaining 3′ read profiles for individual cells, read alignment to the murine mm10 genome reference, and collapsing of unique molecular identifiers (UMIs) were done using the 10X Genomics ‘Cellranger’ software version 3. In order to account for differences in library sizes, all samples were downsampled to the same sequencing depth and combined with the Cellranger’s aggr function. Filtered gene and barcode count matrices were used for scRNA-seq analysis with a custom R pipeline centered around the Seurat package v3.1.1 (Butler et al., 2018; Stuart et al., 2019), in addition to other publicly available R toolkits (see below). To remove potential artifacts due to low quality cells, only cells with > 1000 genes detected, and genes detected in > 5 cells with a minimum of 10 UMI counts per gene were included. High expression levels of mitochondrial genes and the number of detected genes were first used to remove outlier cells using multivariate outlier detection with the mvoutlier R package v2.0.8 (Filzmoser and Gschwandtner, 2015; Filzmoser et al., 2005). The UMI count matrix was then normalized using global scaling log normalization followed by variance stabilizing transformation as implemented in Seurat. In order to control for unwanted sources of variation, the levels of mitochondrial gene expression and number of UMI per cell were regressed out from the filtered and normalized data using linear regression. We analyzed the distribution of cell cycle genes in the data by scoring the expression of canonical cell cycle markers (Kowalczyk et al., 2015) and classifying each cell into S, G1, or G2/M cell cycle phase. We then ran a principal component analysis (PCA) to determine the contribution of cell cycle genes to the top PCs. From this analysis we determined that 37% of genes in the first five top PCs were cell cycle related genes; thus, we proceeded to adjust for the effects of cell cycle. As embryonic stages are expected to contain both proliferating and non-proliferating cells of different cell types, we treated the difference between cell cycle scores as a confounder variable and regressed it out by linear regression (Mayer et al., 2018). This removes differences in cell cycle phase within proliferating cells, while preserving differences between proliferating and non-proliferating cells.

Cell type identification and comparisons between conditions

PCA was performed on highly variable genes identified as those with the highest standardized variance as implemented in Seurat (Stuart et al., 2019). We then determined significant PCs by the Jackstraw method (Buja and Eyuboglu, 1992). As > 20 PCs were significant, we selected the top PCs based on the highest inflection point on the standard deviation of the PCs (elbow plot) with the extremum distance estimator method (Christopoulos, 2016). The top 18 PCs were then used to perform unsupervised shared nearest neighbor graph-based clustering (k parameter = 100) as implemented in Seurat. Clustering was performed at high resolution (r = 1.2). These parameters resulted in the maximum number of clusters, for which each cluster contained distinctive gene expression signatures after differential gene expression testing (see below). We performed UMAP dimensionality reduction (McInnes et al., 2018), in order to visualize the final clustering results overlaid on a UMAP plot. All samples aligned well and no stage-, genotype-, library-, or replicate-specific clusters were observed. Differential gene expression among cell clusters, stages, and genotypes was assessed using the generalized linear model method MAST (Finak et al., 2015), with an FDR ≤ 0.05 and lnFC expression threshold ≤ 0.25. In order to define specific cluster signatures, all possible pairwise comparisons among clusters were tested and the maximum FDR and minimum lnFC values were used to determine significant markers (Haber et al., 2017). To identify clusters we used the top differentially expressed genes as well as known markers of potential suture tissue types (mesenchyme, osteoblasts, hypodermis, dura mater, osteoclasts, capillary endothelial cells, pericytes, and hematopoietic lineages; Table S3). To aid cell type and cell function classifications we performed ranked gene enrichment analysis with the R extension of gProfileR (FDR ≤ 0.05) (Reimand et al., 2007). In addition, we queried the top markers of each cell cluster against indexes built for the Mouse Cell Atlas (Han et al., 2018) for bone marrow, brain, embryonic mesenchyme, embryonic stem cells, mesenchymal stem cells, and neonatal calvaria datasets using the search index Scfind (Lee et al., 2019a). We identified the hypodermis and dura mater clusters by reference to recent single cell studies of these tissues (DeSisto et al., 2019; Driskell et al., 2013; Philippeos et al., 2018; Sennett et al., 2015). Finally, to validate cell type clusters, we then scored each cell for each cell type signature module, obtained from differential gene expression analysis and known markers, as described above for cell cycle analysis and using the AddModuleScore function in Seurat. Cell type specific clusters were manually curated and annotated based on the obtained enrichment scores and their distribution on the UMAP plot. The curated clusters were then tested for a final round of differential gene expression analysis and gene enrichment.

After the initial clustering analysis we noticed that a small number of cells exhibited outlier expression of hemoglobin and erythroid cell-related genes; however none of the clusters was identified as an erythroid population and these cells were contained in clusters identified as unrelated cell types. We treated these cells as potential doublets and removed them. In addition, we treated the expression of erythroid-specific markers as an unwanted source of variation.

Identification of subpopulations within the SM and osteoblast clusters

To analyze the relationship between SM and osteoblasts we performed unsupervised clustering (k parameter = 100) and differential gene expression analysis for these cell types as described above, but using the top nine PCs. Differential gene expression analysis parameters were adjusted to account for a higher expected degree of homogeneity within FS1–5, OB, HD, and DM subpopulations in contrast with the analysis at the cell type level; thus, a combined Fisher p value of 0.01 and a minimum lnFC threshold of 0.1 were used to determine significant markers (Haber et al., 2017). In addition, we conducted trajectory analysis with UMAP and approximate graph abstraction with Monocle 3 (Cao et al., 2019). To remove potential confounders due to differences in cell states, and given the presence of one subpopulation enriched for cell cycle genes, this analysis was performed with and without removing from the data the cell cycle-related genes (GO:0007049). Removal of cell cycle genes did not result in changes in the overall trajectory.

Cell type enrichment analysis and intersections between bulk and single-cell RNA-seq

Differentially expressed genes from regional, temporal, and genotypic comparisons determined by bulk RNA-seq analysis were tested for significant (p ≤ 0.05) gene set enrichment against the cell-type-, stage-, and genotype-specific signatures, determined by scRNA-seq analysis using Fisher’s exact tests and using Bonferroni correction for multiple comparisons. This method was also used to test for enrichment of dissociation-associated signatures in the scRNA-seq data.

To compare enriched gene ontology categories between bulk and single-cell RNA-seq datasets, we focused on selecting the top 5 to 10 most significant categories (FDR ≤ 0.05). Only those categories that were present in both data types were selected for plotting the intersections shown in Figure 4.

Alkaline phosphatase staining

Alkaline phosphatase (ALPL) histochemistry was performed on LCM “guide slides” as described previously with modifications (Miao and Scutt, 2002). Frozen sections were first dried at 37°C for 30 minutes, fixed in 4% paraformaldehyde (PFA) in PBS for 15 minutes, washed in three changes of PBS for five minutes each, incubated with 0.1M Tris-maleate buffer pH9.2 (Sigma-Aldrich, T3128) for five minutes, then incubated with 0.4 mg/ml Fast Red TR (Sigma-Aldrich, F8764) and 0.02% naphthol AS-MX phosphate (Sigma-Aldrich, 855) in 0.1M Tris-maleate buffer for approximately 3–5 s, until red staining was just apparent.

Suture width measurement

Frontal suture widths were determined on LCM “guide slide” sections stained for ALPL. The “Segmented Line” tool of ImageJ was used to trace the length of the suture between osteogenic fronts, and the “Measure” tool of ImageJ was used to quantify the linear distance between them. These widths were determined at three locations along the AP axis: in the same coronal plane of the anterior of the eye, within the midregion of the eye, and the posterior of the eye. At least five mice were measured in each genotype group.

MEGENA network analysis

The same filtered, normalized, and covariate-adjusted gene expression data matrix used in the differential gene expression analysis was also used as input for network analysis with the R package MEGENA v1.3.4 (Song and Zhang, 2015). The Pearson method was used to calculate co-expression correlation between gene pairs, with FDR calculation by permutation analysis (100 iterations) and a cutoff of 0.05 to identify significant correlations. Following planar filtered network construction, multi-clustering analysis and multi-scale hub analysis were performed using 100 permutations. Significant modules and hubs were identified at FDR ≤ 0.05 at optimized resolution parameters. Finally, key driver analysis (Zhang et al., 2013) was performed on each module to identify key regulators of the module. Each module was then tested for significant gene set enrichment (FDR ≤ 0.05) against the Molecular Signatures Database (MSigDB) (Subramanian et al., 2005), DEGs among groups and cell-type-specific genes using a Fisher’s exact test and Bonferroni correction to account for multiple testing, and module diagrams showing the co-expression relationships between modules and genes within each module were plotted with Cytoscape (Shannon et al., 2003).

Lentiviral packaging

Third-generation lentiviral expression vectors (pLV[Exp]-EGFP-EF1A) carrying murine genes of interest were obtained from Cyagen Biosciences. C-terminally FLAG-tagged proteins were expressed from the EF1A promoter. A GFP-TA-puromycin resistance protein was expressed from the CMV promoter (vector maps available on request). Third-generation lentiviral packaging vectors (pRSV-Rev, pMDLg/pRRE, pMD2.G) were obtained from Addgene (12253, 12251, 12259). Lentivirus was packaged by co-transfecting individual expression vectors with associated packaging vectors into HEK293T/17 cells (ATCC, CRL-11268) with Lipofectamine 2000 (Invitrogen, 11668019) in Opti-MEM I Medium containing 10% FBS. Cells were maintained at 37°C in a humidified 5% CO2 incubator. On the second day the medium containing the DNA-Lipofectamine 2000 complexes was removed and replaced with 20 mL of complete culture medium without antibiotics. Medium containing lentivirus was harvested on the third and fourth days after transfection and purified by centrifugation at 3,000 rpm for 15 minutes followed by filtration through a 0.45 μm PVDF filter.

Lentiviral cell lines

MC3T3-E1 Subclone 4 cells (ATCC, CRL-2593) were transduced with lentivirus-containing cell culture medium in the presence of 6 μg/ml polybrene, and incubated at 37°C in a humidified 5% CO2 incubator for 24 hours. Transduced cells were selected with 1 μg/ml puromycin for 10–12 days until no further cell death was observed. To obtain a more homogeneous population expressing the target gene at higher levels, each cell line was FACs-sorted by GFP signal selection before passage 5 using an IMI5L cell sorter (BD Biosciences). Cells with a high GFP fluorescence signal intensity were isolated using a threshold of GFP fluorescence based on background fluorescence signals of the non-transduced MC3T3-E1 cell line using FACSDiva 8.0.2 (BD Biosciences). Cells expressing high levels of GFP were collected, expanded, and cryopreserved for future experiments.

Osteoblast differentiation assay

MC3T3-E1 cells expressing genes of interest were plated at 3×103 cells/cm2 in 6-well plates in αMEM complete growth medium (αMEM, 10% FBS, 1% non-essential amino acids, 1% penicillin/streptomycin). At confluence cells were incubated in osteoblast differentiation medium (αMEM, 10% FBS, 1% non-essential amino acids, 1% penicillin/streptomycin supplemented with 10 mM β-glycerol phosphate, 200 μM ascorbic acid, and 0.1 μM dexamethasone). Culture medium was changed every 2–3 days. Mineralization was visualized by staining with alizarin red and von Kossa reagents by standard means.

qRT-PCR

To quantify gene expression in stably-transfected cell lines, MC3T3-E1 cells expressing genes of interest were plated in triplicate at 3×103 cells/cm2 in 6-well plates as described above. Upon confluence RNA was prepared using the QIAGEN RNeasy Kit as described by the manufacturer. RNA concentration and purity were determined by Nanodrop. Single-stranded cDNA was synthesized with oligo(dT) primers using the SuperScript III First-Strand Synthesis System (Invitrogen, 18080051). Primers for genes of interest were obtained from the literature (Demitrack et al., 2015; Ikegawa et al., 2008; Koch et al., 2004), Icahn School of Medicine qPCR Core, or designed using Primer3 (see Table S6). RT-qPCR targets were amplified in triplicate using SYBR Green and Platinum Taq polymerase (Thermo Fisher Scientific, S7567) on a 7900HT Real-Time PCR instrument (Thermo Fisher Scientific, 10966). To quantify differential gene expression and splicing events identified by bulk RNA-seq, RT-qPCR targets were amplified in triplicate from 5 ng of amplified region-specific cDNA prepared in “Bulk RNA-seq library preparation and sequencing.” The relative expression of each gene was calculated by the dCt method, using Rps11 (differential gene expression) or Actb (differential splicing and transfected cell lines) as the reference.

RNA in situ hybridization

Differential gene expression identified by bulk RNA-seq comparisons was validated by RNA in situ hybridization of cryoembedded and sectioned embryos. Embryos were fixed in 4% PFA in PBS at 4°C overnight, washed in PBS, equilibrated with 30% sucrose/PBS, and embedded in OCT compound. All PBS solutions were pre-treated with DEPC. Frozen sections were cut on an Avantik QS11 cryostat at 12 μm and transferred to Superfrost/Plus microscope slides (Thermo Fisher Scientific). Riboprobe templates were generated by PCR with primers from published literature or designed by Primer3 (see Table S6), using cDNA derived from embryonic murine brain RNA. Riboprobes were prepared with DIG RNA Labeling Mix as described by the manufacturer (Sigma-Aldrich, 11277073910). RNA in situ hybridization (ISH) was performed essentially as described (Wilkinson, 1998). Cryosectioned slides were hybridized with riboprobes at 65°C overnight, incubated with anti-digoxigenin-AP antibody (Sigma-Aldrich, 11093274910) at 4°C overnight, and hybridization was visualized by BM Purple staining (Sigma-Aldrich, 11442074001) at room temperature for periods of 1 hour to 5 days depending on the intensity of the color reaction.

Single molecule fluorescent RNA in situ hybridization

Single molecule fluorescent RNA in situ hybridization (smFISH) was performed using the RNAscope Fluorescent Multiplex Reagent Kit (Advanced Cell Diagnostics, 320850) as described by the manufacturer with the following modifications: 4% PFA fixation of sections was performed for one hour, and Proteinase IV digestion was performed for 10 minutes. Probes (Advanced Cell Diagnostics) were for Acta2 (319531-C3), Clec3b (539561-C2), Col8a1 (518071), Cxcl12 (422711-C3), Ibsp (415501), Icam1 (438611) Npnt (316771-C2), and Top2a (491221). E18.5 C57BL/6J sections of 10 μm were from cryoembedded fresh frozen embryos. Widefield images were acquired on an AxioImager Z2M equipped with a 20x/0.8NA Zeiss Plan-Apochromat objective, a monochrome Axiocam 503 camera (Zeiss, 1936 × 1460 pixels, 4.54 μm × 4.54 μm per pixel, sensitivity ~400 nm–1000 nm) and Zen 2 Blue Edition software (version 2.0). For larger fields-of-view tiled images were stitched together using Zeiss Zen Blue software. Z stack images were acquired at optimal sampling rate (meeting Nyquist frequency requirements), as calculated by the software.

QUANTIFICATION AND STATISTICAL ANALYSIS

Specific statistical analyses are described in the STAR Methods section and Figure legends. Image analyses were performed with ImageJ (National Institutes of Health). GraphPad Prism or Microsoft Excel were used to create graphs indicating mean ± standard deviation (SD) and calculate statistical significance using the unpaired Student’s t test, where p values ≤ 0.05 were considered significant. Sample size (n) indicating the number of individual mice or tissue culture wells used is specified in Figure legends.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Anti-Digoxigenin-AP	Sigma-Aldrich	11093274910; RRID:AB_514497
Chemicals, Peptides, and Recombinant Proteins
Cresyl violet	Sigma-Aldrich	C5042
Collagenase Type II	Worthington	LS004176
Dispase II	Sigma-Aldrich	4942078001
Red Blood Cell Lysis Solution	Miltenyi Biotec	130–094-183
SIGMAFAST Fast Red TR/Naphthol AS-MX Tablets	Sigma-Aldrich	F4523–50SET
β-Glycerophosphate disodium salt hydrate	Sigma-Aldrich	G9422
L-Ascorbic acid 2-phosphate trisodium salt	Wako Chemical Ltd.	321–44823
Dexamethasone	Sigma-Aldrich	D4902
Alizarin Red Staining Solution	EMD Millipore	TMS-008-C
Lipofectamine 2000 Transfection Reagent	Invitrogen	11668019
DIG RNA Labeling Mix	Sigma-Aldrich	11277073910
BM Purple	Sigma-Aldrich	11442074001
Critical Commercial Assays
Arcturus PicoPure RNA Isolation Kit	Applied Biosystems	KIT0204
Agilent RNA 6000 Pico Kit	Agilent Technologies	5067–1513
Chromium Single Cell Gene Expression Kit (Version 3)	10X Genomics	1000075, 1000078
Nextera XT DNA Sample Preparation Kit	Illumina	FC-131–1096
Anti-BrdU Antibody	Sigma-Aldrich	RPN202; RRID:AB_2314032
Von Kossa Stain Kit	American MasterTech Scientific	NC9239431
QIAGEN RNeasy Mini Kit	QIAGEN	74104
SuperScript III First-Strand Synthesis System	Invitrogen	18080051
Sybr Green kit	Thermo Fisher Scientific	S7567
Platinum Taq Polymerase	Thermo Fisher Scientific	10966
RNAscope Multiplex Fluorescent Detection Kit v2	Advanced Cell Diagnostics	320850
Deposited Data
Bulk mRNA expression data for the initial set of five replicate WT and Twist1+/− library sets	FaceBase (https://www.facebase.org/chaise/record/#1/isa:dataset/accession=FB0000C)	Accession FB00000805
Bulk mRNA expression data for the initial set of five replicate Fgfr2+/S252W;EIIACre+ library sets	FaceBase (https://www.facebase.org/chaise/record/#1/isa:dataset/accession=FB0000094)	Accession FB00000904
Single cell mRNA expression data for the E16.5 and E18.5 frontal sutures	FaceBase (https://www.facebase.org/chaise/record/#1/isa:dataset/accession=FB00001013)	Accession FB00001013
Experimental Models: Cell Lines
HEK293T/17	ATCC	CRL-11268; RRID:CVCL_1926
MC3T3-E1 Subclone 4	ATCC	CRL-2593; RRID:CVCL_5440
Experimental Models: Organisms/Strains
C57BL/6J mice	The Jackson Laboratory	000664; RRID:IMSR_JAX:000664
Twist1+/− mice	(Chen and Behringer, 1995)
Fgfr2+/NeaS252W mice	(Wang et al., 2005)
EIIACre mice	(Lakso et al., 1996)
Oligonucleotides
See Table S6
Recombinant DNA
Lgr5 ISH plasmid: gift from Linda Samuelson (University of Michigan); AS: NdeI/T7; S: NotI/Sp6	(Demitrack et al., 2015)
Third-generation lentiviral packaging plasmid, pRSV-Rev	Addgene	12253
Third-generation lentiviral packaging plasmid, pMDLg/pRRE	Addgene	12251
Third-generation lentiviral packaging plasmid, pMD2.G	Addgene	12259
pLV[Exp]EGFP:T2A:PuroEF1A > mCpxm2/FLAG	this paper
pLV[Exp]EGFP:T2A:PuroEF1A > mLgr5/FLAG	this paper
pLV[Exp]EGFP:T2A:PuroEF1A > mCherry/FLAG	this paper
Software and Algorithms
10X Genomics Cellranger	10X Genomics	https://www.10xgenomics.com/
mvoutlier R package v2.0.8	(Filzmoser and Gschwandtner, 2015)	https://cran.r-project.org/web/packages/mvoutlier/index.html
Seurat package v3.1.1	(Butler et al., 2018; Stuart et al., 2019)	https://satijalab.org/seurat/, RRID:SCR_016341
Monocle 3	(Cao et al., 2019)	https://cole-trapnell-lab.github.io/monocle3/
gProfileR	(Reimand et al., 2007)	https://cran.r-project.org/web/packages/gProfileR/index.html
STAR	(Dobin et al., 2013)	https://github.com/alexdobin/STAR/, RRID:SCR_015899
featureCounts	(Liao et al., 2014)	http://bioinf.wehi.edu.au/featureCounts/, RRID:SCR_012919
Bioconductor Limma package	(Ritchie et al., 2015)	http://bioinf.wehi.edu.au/limma/, RRID:SCR_010943
RSEM package	(Li and Dewey, 2011)	http://deweylab.biostat.wisc.edu/rsem/, RRID:SCR_013027
LeafCutter v1.0 (Li et al., 2018)	(Li et al., 2018)	https://github.com/davidaknowles/leafcutter/
MEGENA v1.3.4 R package	(Song and Zhang, 2015)	https://cran.r-project.org/web/packages/MEGENA/index.html
ImageJ	National Institutes of Health	https://imagej.nih.gov/ij/
FACSDiva 8.0.2	BD Biosciences
Primer3		http://primer3.ut.ee; RRID:SCR_003139
Cytoscape v3.3.0	(Shannon et al., 2003)	https://cytoscape.org/