The Fgfr2c ^C342Y/+ mice were a generous gift from Dr. Eswarakumar, who previously described the generation of this model (Eswarakumar et al., 2004). The model was maintained on a CD1 background for viability and breeding. Embryos were collected at embryonic days E13.5, 15.5, and 17.5. PCR for sex and genotype was performed using a tail snipped at collection. We do not anticipate any differences due to sexual dimorphism based on our previous work on multiple mouse models (Motch Perrine et al., 2017; Motch Perrine et al., 2019; Lesciotto et al., 2022), still we have indicated the balance of males and females used in this study in Table 1. Mouse litters were produced in compliance with animal welfare guidelines approved by Icahn School of Medicine at Mount Sinai and Pennsylvania State University Animal Care and Use Committees (ISMMS #07-0757 and PSU #46558). To more precisely define developmental age, E13.5 and E15.5 embryos were staged using the Embryonic Mouse Ontogenetic Staging System, EMOSS (Musy et al., 2018). E17.5 embryos have developed beyond the capabilities of this system so are considered to be an average of 17.5 days post conception. Average developmental age at E13.5 was 341.14 h for Fgfr2c ^C342Y/+ embryos and 338.4 h for Fgfr2c ^+/+ embryos. At E15.5, average developmental age was 371.93 h for Fgfr2c ^C342Y/+ embryos and 372.11 h for Fgfr2c ^+/+ embryos.

Phosphotungstic acid (PTA)-enhanced Microcomputed Tomography (microCT) images for MC and microCT images for bone analyses were acquired by the Center for Quantitative Imaging at the Pennsylvania State University (https://iee.psu.edu/labs/center-quantitative-imaging) using the General Electric v|tom|x L300 nano/microCT system. Five Fgfr2c ^C342Y/+ and 5 Fgfr2c ^+/+ embryos were scanned at E13.5, 5 Fgfr2c ^C342Y/+ and 6 Fgfr2c ^+/+ embryos were scanned at E15.5, and 4 Fgfr2c ^C342Y/+ and 5 Fgfr2c ^+/+ embryos were scanned at E17.5. Image data were reconstructed on a 2,024 × 2,024 pixel grid as a 32-bit volume, and were reduced to 16-bit volume for image analysis using Avizo 2019.3 and 2020.2 (Thermo Fisher Scientific). Scanning parameters varied from 60–100 kV and 75–170 μA, to accommodate age group and type of scan performed. Voxel sizes ranged from 6.9 to 15 microns (µm) for bone scans and 4.5–8 µm for PTA-enhanced scans.

3D models of MC were created from PTA-enhanced microCT (Lesciotto et al., 2020) images of specimens at E13.5, E15.5 and E17.5 (Table 1) using a previously reported automatic segmentation approach (Zheng et al., 2020). Due to the complex 3D cartilage structures and variations in our images as well as large volumes of image data, an effective and efficient automatic method is needed for our cartilage segmentation task. Recently, deep learning methods have achieved great successes in biomedical image segmentation (Ronneberger et al., 2015). Compared with common CT images, each of our 3D microCT volumes is very large (about 1,500 × 2,000 × 2,000 voxels per image), and consequently, it will be extremely difficult to rely on manual annotation to produce a sufficiently large amount of labeled data for training a deep learning segmentation model. Hence, our segmentation approach must be able to accommodate a large image size and make the most out of sparsely labeled training data. First, we chose a very sparse subset of 2D image slices (e.g., 2%–10%) in our 3D microCT training volumes for manual annotation that represents and covers the rest (unannotated) 2D slices of the training volumes well. Then, sparse manual annotations were conducted by experts (MKP and SMMP) using Avizo 2020.2 (Thermo Fisher Scientific). Next, the annotated 2D slices were used to train a judiciously designed fully convolutional network (FCN) model to identify target objects (called pseudo-labels or PLs) in all the unannotated 2D slices of the training volumes. In this way, we sought to bridge the large gap between our sparse annotation and full annotation of the entire 3D training volumes. Since our FCN model is not trained with a standard full annotation protocol (due to the sparse annotation used), its identified PLs on the unannotated 2D slices may be unreliable and noisy. Thus, we estimated the reliability of the PLs by computing the associated uncertainty maps of the PLs, which quantify the pixel-wise prediction confidence. With the expanded training set (i.e., 3D images with both manual labels and pseudo labels) and guided by the uncertainty maps, the FCN model was trained iteratively to distill more stable knowledge about the training data, thus becoming more robust and generalizing better to unseen data. Finally, the segmentation results were integrated along the three orthogonal planes (orthogonal to each of the x, y, and z directions) to further boost the overall segmentation accuracy (by average ensemble of segmentation results along the three directions). Results of automated segmentation were reviewed by an expert anatomist and segmentation was corrected manually if necessary. Volume and area were estimated using the Material Statistics of Avizo 2020.2 (ThermoFisher Scientific). Volume and area were compared between the genotypes within age categories with independent-samples Mann Whitney U test performed in IBM SPSS Statistics v 27 (IBM Corp). Once created, the 3D models of MC were quantified by collecting the 3D coordinates of landmarks designed for age-specific analyses (Figure 2; Table 2). Landmark collection as well as volume and area calculations were performed in Avizo 2019.3 and 2020.2 (Thermo Fisher Scientific). Each sample was digitized twice by the same observer, checked and corrected for gross error, and remaining measurement error was minimized by averaging the coordinates of the two trials. A maximum of 5% intraobserver error in landmark placement was accepted. Morphological differences were assessed as described in Section 2.5. PTA-enhanced microCT images were used to hand segment tooth germs visualized in Figure 1.

Isosurfaces of mandibular bone were created by automated segmentation of bone from 16-bit microCT images using a minimum threshold of 70–100 mg/cm³ partial density hydroxyapatite based on phantoms scanned with the samples using Avizo 2020.2 (ThermoFisher Scientific). Data were processed using a median filter in Avizo 2020.2 (ThermoFisher Scientific) to remove noise and reviewed by an expert anatomist. No manual corrections to form were necessary for MB. Volume and area were estimated for each specimen using the Material Statistics module of Avizo 2020.2 (Thermo Fisher Scientific). Volume and area were compared between the genotypes within age categories with independent-samples Mann Whitney U test performed in IBM SPSS Statistics v 27 (IBM Corp). These isosurfaces were used to assess mandibular form by recording the 3D coordinates of 14 (7 bilateral) landmarks at E15.5 and 20 (10 bilateral) landmarks at E17.5 (Figure 2; Table 3) in Avizo 2019.3 (Thermo Fisher Scientific). Each sample was digitized twice by the same observer, checked and corrected for gross error, and measurement error was minimized by averaging the coordinates of the two trials. A maximum of 5% intraobserver error in landmark placement was accepted. Morphological differences were assessed as described in Section 2.5. Sample sizes for each genotype and embryonic age are listed in Table 1.

To statistically determine differences in MC and mandibular form between genotypes, we used Euclidian Distance Matrix Analysis (EDMA) (Lele and Richtsmeier, 1995; Lele and Richtsmeier, 2001). EDMA converts 3D landmark data into a matrix of all possible linear distances between unique landmark pairs and tests for statistical significance of differences between shapes using a boot-strapped hypothesis testing procedure. Differences of specific linear distances are statistically evaluated by a 90% confidence interval produced through a non-parametric bootstrapping procedure (Lele and Richtsmeier, 1995). Rejection of the null hypothesis of similarity for linear distances enables localization of differences to specific dimensions. EDMA analyses were performed using WinEDMA (Cole, 2002).

Scans targeting soft tissues such as MC (PTA-enhanced microCTs) and scans targeting hard tissues such as MB were performed on each specimen as described in Section 2.2. MB is very difficult to fully segment from PTA-enhanced microCTs, so therefore both scan types were used. For figures with superimposition of MB and MC, isosurfaces of MC and MB were created as described in Sections 2.3,2.4, respectively. Landmarks 1, 2, and 15 (Table 3) were taken on MB in both scans in Avizo 2020.2 (ThermoFisher Scientific) for each individual for which a superimposition was required. A rigid transformation, which moves the points of the first set as close as possible onto the points of the second set (the sum of the squared distances between corresponding points is minimized) was computed from these landmarks for aligning the model (bone microCT scan) to the reference (PTA-enhanced microCT scan) in Avizo 2020.2 (ThermoFisher Scientific). The superimposition was adjusted slightly by hand if necessary following landmark transformation of images. Superimpositions were used only for figure creation. Analysis of individual tissues was performed on individual scans targeting the tissue of interest.

Specimen were collected (n = 3 per age per genotype) and fixed in 4% paraformaldehyde for 24 h. Care was taken to balance between sexes, and E17.5 specimens were decalcified in 0.25 M EDTA at pH 7.4 for 12 days. Isolated heads were processed for paraffin-based histology per established protocol (Behringer et al., 2013; Motch Perrine et al., 2021). Using a rotary microtome, 7 µm serial, coronal sections were cut and mounted on Superfrost-Plus slides (Thermo Fisher Scientific), and the range of slides containing Meckel’s cartilage was determined for each specimen. Slide ranges were then divided into four equal regions (Anterior, Mid 1, Mid 2, Posterior, Figure 1), and one slide was selected per region per specimen for each stain. For immunohistochemistry, slides were subjected to epitope retrieval using an autoclave (121°C, 10 min) and sodium citrate buffer pH 6.0 before blocking with 3% hydrogen peroxide and 1% goat serum/bovine serum albumin. Sections were incubated with primary antibodies at room temperature for 2 h, Proliferating Cell Nuclear Antigen (PCNA, AbCam, ab18197, 1:2,000) or for 1 h, alkaline phosphatase (ALP, ab108377, 1:500). After washing in PBS, slides were incubated with HRP conjugated secondary antibody for one hour at room temperature (ab6721, 1:500) and Diaminobenzidine (DAB, Vector Laboratories) chromagen was used to identify immunoreactive structures (Durham et al., 2019). Safranin-O and tartrate-resistant acid phosphatase (TRAP) staining were performed per established protocol (Howie et al., 2018). Slides were imaged using a Leica BX50 microscope, DFC450 camera, and LAS-X imaging software (Leica Biosystems). Regions of interest (MC, mandible) were identified and analyzed using ImageJ color deconvolution and masks to count stained areas by color (Supplementary Figure S1) (Ruifrok and Johnston, 2001; Haub and Meckel, 2015; Durham et al., 2018; Landini et al., 2021). Careful gating for cell size (3–30 µm) was used to reduce noise in the form of stained areas smaller or larger than a cell. The same gating parameters were applied to every image analyzed. At least two images were captured per slide for each stain for each specimen investigated. Osteoclasts were identified by TRAP positivity and multiple nuclei and were counted by a single investigator (ELD) two times with significant correlation between counts. Counts were averaged for investigation. Non-parametric Mann-Whitney U tests were used where necessary to compare genotypes at each age, and then by each region using SPSS 25 software (IBM).

Specimens (n = 3 per genotype) collected at E15.5 were fixed in 4% paraformaldehyde overnight, embedded in OCT compound, and frozen. Samples were sectioned at 10 µm of thickness on the cryostat. For immunofluorescence, sections were permeabilized in PBST (0.2% Triton X 100 in PBS) for 5 min. Antigen retrieval was performed at 100°C for 10 min in citrate buffer (Sigma-Aldrich, C9999). Sections were then blocked with 3% bovine serum albumin in PBS. Anti-p-ERK1/2 (Cell Signaling Technology, 4370) was diluted at 1:100 with PBS and incubated on the sections for 1 h. After washing with PBS, secondary antibody (1:500, Thermofisher, A32732) was incubated for 1 h. The sections were rinsed briefly in PBS and stained with 0.5 μg/ml 4′,6-Diamidino-2-Phenylindole (DAPI; Thermo Fisher Scientific) in PBS before mounting in anti-fade mounting medium (Vector Labs, H-1700). Slides were imaged using a Nikon Eclipse E 600 fluorescence microscope. Images were analyzed using ImageJ.

Visual inspection of superimpositions of bone and cartilage segmentations indicated only subtle morphological differences between Fgfr2c ^C342Y/+ and Fgfr2c ^+/+ embryos from E13.5 through E17.5 (Figure 3). Very thin areas of bone with little mineralization may not be detected at the indicated bone threshold in microCT images, as compared to traditional histological staining techniques.

At E13.5, Fgfr2c ^C342Y/+ and Fgfr2c ^+/+ embryos do not exhibit significant differences in MC form. However, at E15.5 MC volume is significantly larger in these embryos as compared to Fgfr2c ^+/+ littermates (Table 4) and MC is significantly longer in the anteroposterior dimension (Figure 4A). MC of Fgfr2c ^C342Y/+ embryos is not larger overall at E17.5 but linear distances characterizing the most anterior and most posterior aspects of MC are larger relative to unaffected littermates (Figure 4B).

Effects of the Fgfr2c C342Y mutation on chondrocytes within MC are not qualitatively apparent, especially at the earliest time points investigated here (Figure 5A). Since asymmetry was not noted in our assessment of MC form, we did not separate left and right in our histological assessment. The area of MC measured in coronal sections of E13.5 embryos revealed increased area of MC in Fgfr2c ^C342Y/+ Crouzon syndrome embryos (p = 0.001) as compared to Fgfr2c ^+/+ littermates. This was driven by significantly increased area in Mid 1 (p = 0.003) and Posterior (p ≤ 0.001) regions of MC in Fgfr2c ^C342Y/+ embryos at E13.5. No statistically significant differences in cross-sectional area were identified at the later timepoints (E15.5, E17.5) (Figure 5B). Since cartilage growth can result from increased matrix deposition (Figure 5C), net increase in cells (Figure 5D), or an increase in chondrocyte hypertrophy (cell diameter) (Figure 5E), we investigated these characteristics in MC and identified more matrix area in Fgfr2c ^C342Y/+ embryos at E13.5 (p ≤ 0.001) relative to Fgfr2c ^+/+ littermates. This was the result of greater matrix in all regions in Fgfr2c ^C342Y/+ embryos, but Mid 1 (p = 0.003), Mid 2 (p = 0.001), and Posterior (p = 0.001) regions had statistically significantly more cartilage matrix as compared to Fgfr2c ^+/+ littermates (Figure 5C). Interestingly, there were no differences in chondrocyte number between genotypes at E13.5 or E15.5. The number of chondrocytes within MC was statistically significantly decreased in Posterior of Fgfr2c ^C342Y/+ embryos at E17.5 (Figure 5D). Our assessment of the diameter of chondrocytes within MC indicated no differences by genotype at E13.5 and E15.5, however there were significantly smaller cells in Anterior of MC in Fgfr2c ^C342Y/+ embryos at E17.5 (Figure 5E). This result correlates with the trend toward more cells in Fgfr2c ^C342Y/+ embryos in Anterior at E17.5 (Figure 5D).

Because we identified some significant differences in cell number between genotypes, we investigated proliferation to see if the mutation influenced chondrocyte and/or bone lineage cells during initial ossification. E15.5 and E17.5 embryos were used because the bone was robust enough at these ages for analysis. Differences in proliferation in MC are subtle (Figure 6A; Supplementary Figure S1), however an overall decrease in proliferation in chondrocytes from Fgfr2c ^C342Y/+ embryos was identified at E15.5 (Figure 6B, top) relative to Fgfr2c ^+/+ littermates. This relationship was no longer significant by E17.5 (Figure 6B, row 2). Within the forming mandible adjacent to MC, proliferation was increased in Anterior, and decreased in Posterior of E15.5 Fgfr2c ^C342Y/+ embryos relative to Fgfr2c ^+/+ littermates (Figure 6B, row 3). At E17.5 an overall decrease in proliferation in the mandibles of Fgfr2c ^C342Y/+ embryos was driven by significantly less proliferation in Posterior (Figure 6B, bottom).

As demonstrated by EDMA of 14 anatomical landmarks (Figure 2B; Table 3), the mandible is generally smaller in the antero-posterior dimension of Fgfr2c ^C342Y/+ embryos relative to Fgfr2c ^+/+ littermates at E15.5 (Figure 7A). EDMA of 20 anatomical landmarks (Figure 2B; Table 3) indicated differences between genotypes are concentrated in the posterior of the developing mandible corresponding to an area that stretches from the posterior portion of Mid 1 through the rostral portion of Posterior of MC that includes the coronoid, condylar, and angular processes (Figure 7B). There was no overall difference in bone volume at E15.5 or E17.5 (Table 4).

To assess differentiation of osteoblasts, which are responsible for morphological changes in the mandible of Fgfr2c ^C342Y/+ embryos, we investigated the levels of ALP-positive osteoblasts, which are responsible for bone mineralization (Hall, 2014), but identified no significant effects of this mutation on the ratio of alkaline phosphatase positive osteoblasts to total cells in the mandible (Figures 8A,B). In contrast to the lack of effect on the number of bone building osteoblasts, bone resorbing osteoclasts (measured by TRAP staining and multi-nucleation) were reduced in Fgfr2c ^C342Y/+ embryos at E17.5 relative to Fgfr2c ^+/+ littermates. This overall reduction in osteoclasts was driven by a significant reduction in osteoclasts in Mid 1 of Fgfr2c ^C342Y/+ embryos (Figures 8 C,D).

To assess how Fgfr2c C342Y alters the ERK pathway in the lower jaw, we analyzed the expression of p-ERK1/2 in MC and mandible at E15.5, the earliest timepoint when the phenotypic effects of Fgfr2c C342Y mutation on MC and mandible were observed, prior to significant degradation in MC (Sakakura et al., 2007). p-ERK1/2 was detected in the mandible, tongue, tooth, peripheral chondrocytes in MC and perichondrium of MC by immunostaining (Figure 9). Little phosphorylated ERK was identified within MC, however there was significant activation in the perichondrium particularly in Anterior, Mid 2, and Posterior of MC, though this activation was reduced in E15.5 Fgfr2c ^C342Y/+ embryos as compared to Fgfr2c ^+/+ littermates. For mandibular bone, the expression level and pattern of p-ERK1/2 in Anterior and Mid 1 of Fgfr2c ^C342Y/+ embryos were similar to Fgfr2c ^+/+ littermates, however, ERK activation is increased in Mid 2 and Posterior. Our assessment of the downstream effects of the Fgfr2c mutation indicates that the ERK pathway is activated particularly in the posterior portion of the mandible of Fgfr2c ^C342Y/+ embryos at E15.5 (Figure 9), matching the area of greatest morphological differences at that age (Figure 7).