Effects of Multi-Generational Soft Diet Consumption on Mouse Craniofacial Morphology

Variations in craniofacial morphology may arise as a result of adaptation to different environmental factors such as soft diet (SD), which lessens functional masticatory load. Prior studies have shown that changes in the masticatory muscle function associated with a switch to short-term SD led to changes in craniofacial morphology and alveolar bone architecture. However, the long-term effects of SD and the associated adaptive changes in craniofacial shape are unclear. Our novel study set out to profile prospective skull changes in mice fed with SDs over multiple generations using three-dimensional (3D) geometric morphometric analysis (GMA). Our results revealed that short-term SD consumption led to a significant decrease in craniofacial size, along with numerous shape changes. Long-term SD consumption over 15 continuous generations was not associated with changes in craniofacial size; however, shape analysis revealed mice with shortened crania and mandibles in the anteroposterior dimension, as well as relative widening in the transverse dimension compared to the average shape of all mice analyzed in our study. Moreover, changes in shape and size associated with different functional loads appeared to be independent – shape changes persisted after diets were switched for one generation, whereas size decreased after one generation and then returned to baseline size. Our study is the first to study the role of prolonged, multi-generational SD consumption in the determination of craniofacial size and shape.


INTRODUCTION
Can exclusive and prolonged consumption of soft diets (SDs) alter our and our offsprings' craniofacial morphology? Human skulls have undergone morphological transformations over centuries with changes being most pronounced in areas associated with masticatory muscle attachments (Jantz, 2001;Godde, 2015;Manthey et al., 2017). Changes in human skull form (shape and size) have been associated with transition to softer diets, from ancestral hunter-gatherer culture to that of agriculture and farming (Larsen, 2006;Pinhasi et al., 2015). Furthermore, transition to softer diets following the Industrial Revolution has been hypothesized to be one of the primary reasons for adaptive changes to craniofacial structures (Larsen, 2006;von Cramon-Taubadel, 2011;Rando et al., 2014). Neo-Lamarckian evolution suggests that the environment can directly alter phenotype and this acquired trait is passed down in a heritable manner through environmental epigenetics and epigenetic transgenerational inheritance (Skinner, 2015). Furthermore, parental diets have been shown previously to potentially influence the presence of a wide range of craniofacial deformities in offsprings (Hassan et al., 2019;Wang et al., 2019).
Morphological changes due to prolonged, multi-generational SD, if diet is indeed an important factor, are clearly not isolated to a single generation. Furthermore, 2D morphometric analysis is insufficient to capture variation across complex craniofacial structures as the essential units of morphology are shape and size, which are difficult to accurately assess with linear and angular measurements. The limitations of 2D analysis may be overcome by utilizing micro-computed tomography (microCT) with three-dimensional (3D) geometric morphometric analysis (GMA) (Bouvier and Hylander, 1982;Sylvester, 2013). 3D GMA utilizes the placement of landmarks to analyze 3D shape variation between specimens in an unbiased manner and is central to evaluating functional and evolutionary hypotheses (Slice, 2005;Freudenthaler et al., 2017). However, landmark-based GMA still suffers from limitations in the characterization of organismal form. More complex structures such as the condyle and glenoid fossa are more challenging to accurately quantify using traditional landmarks, as it can be challenging to identify sufficient landmarks to represent 3D shape. Thus, semi-/pseudo-landmarks were introduced as additional tools to characterize complex 3D shapes (Mitteroecker et al., 2004;Gunz et al., 2005). These approaches enhanced the representation of complex shapes, by densely sampling the regions that may not have many discrete points of geometric homology within or between them but represent homologous structures across specimens.
Thus, the main advantage of semi-landmarks is the ability to readily represent homologous curves and surfaces by sets of points to establish geometric homology between different subjects or specimens Goswami et al., 2019).
The aim of our novel study was to identify the effects of prolonged, multi-generational SD consumption on mouse craniofacial morphology using 3D GMA. Toward this end, we fed mice HD (normal mouse chow) vs. SD over 15 continuous generations and performed 3D GMA and 2D linear morphometrics on mouse crania and mandibles. Semilandmarking was utilized on the condyle. Short-term SD resulted in a decrease in mouse craniofacial size. Prolonged SD did not affect the size of the cranium and mandible but shape changes, such as widening in the transverse dimension and decrease in the anteroposterior dimension of both the cranium and mandible were noted.

Animals and Experimental Design
FVB mice (Charles River), which are inbred and generate a relatively high number of pups (6-8) per litter, were utilized for these studies. At the initiation of our study, pregnant FVB mice at embryonic (E) day 0.5 (i.e., vaginal plug was observed) were maintained on HD (normal mice chow) or switched to SD (PicoLab 5058, LabDiet, Deans Animal Feeds, Redwood City, CA, USA). This diet regimen was maintained continuously for 15 (F15SD) generations. With each new generation, non-littermate males and females (two mating pairs) were randomly selected (runts were excluded) and paired for mating to continue the HD or SD mouse lines. HD Con and F15SD mice were sacrificed at the same time. We then switched HD Con and F15SD mice to SD (F1SD) and HD (F15SD-F1HD), respectively, for one generation. F1SD and F15SD-F1HD mice were sacrificed at the same time. HD and SD are equivalent in nutrient content, differing only in hardness (Deans Animal Feeds). Sample sizes in these groups were HD Con (N = 12; 6 males and 6 females), F15SD (N = 12; 6 males and 6 females), F1SD group (N = 8; 4 males and 4 females), and F15SD-F1HD (N = 10; 5 males and 5 females) ( Figure 1A). All mice were weaned at 21 days. At 6 weeks of age, mice were weighed, sacrificed, and fixed in 4% paraformaldehyde at 4°C for 48 h. All aspects of animal care and experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at the University of California San Francisco (UCSF) and performed under animal research protocol number AN164201.

Image Processing and Landmark Data Collection Protocol
Fixed mouse heads were imaged by microCT using a SkyScan 1076 MicroCT at the Small Animal Tomographic Analysis Facility (SANTA) located at Seattle Children's Research Institute, United States. Specimens were scanned at 17.2 micron resolution (55 kV, 150 mA, 0.5 mm Al filter). Reconstructions were generated using NRecon (Version 1.6.9.4) with consistent Frontiers in Physiology | www.frontiersin.org 3 July 2020 | Volume 11 | Article 783 A B C FIGURE 1 | Continued thresholding parameters, and converted to 3D volumes. Skull segmentation and the generation of bony surfaces from microCT data were performed with Avizo Lite (version 9.1.1).

Landmarking
We established three sets of landmarks: Set 1 included 13 paired bilateral landmarks that characterized the morphology of the mandible, Set 2 included 43 landmarks that characterized the morphology of the cranium (Figure 1B), and Set 3 included 50 semi-landmarks distributed across the right condyle of all specimens -nine of these points were placed along the condylar head and neck along an imaginary line passing from the tip of the coronoid process to the gonial angle ( Figure 1C). We applied semi-landmarking techniques to quantify 3D condylar shapes since semi-landmarks can significantly improve the characterization of complex structures .

Principal Component and Canonical Variate Analyses
Landmark coordinates were exported as text files and imported to MorphoJ software to perform GMA (Klingenberg, 2011). For each set, landmarks were superimposed using generalized Procrustes analysis (GPA) and the resulting Procrustes coordinates were assessed using principal component analysis (PCA) and canonical variate analysis (CVA) (Slice, 2005). These procedures minimize the influence of size and adopt a single orientation for all specimens by shifting the landmark configurations to a common position, scaling them to a standard size, and rotating them until a best fit of corresponding landmarks is achieved. GPA, therefore, locks the effects of scale, translation, and rotation but does not eliminate the allometric shape variation that is related to size (Rohlf, 1996;Slice, 2005). For this study, variations in craniofacial shape for these samples were assessed using PCA. We carried out PCA using the residuals of multivariate regression of Procrustes coordinates on centroid size to investigate the shape variation independent of size (Darroch and Moismann, 1985;Falsetti et al., 1993). PCA of Procrustes coordinates is based on an eigenvalue decomposition of a covariance matrix, which transforms the Procrustes coordinates into scores along with principal components (PCs). In most cases, the first few PCs explain most of the variance in the dataset. Each observation is scored for each principal axis and the scores of an observation along the principal axes map that observation into the morphospace using MorphoJ (Klingenberg, 2011).
In addition, we carried out CVA to identify shape features that maximize the separation between groups using variances within each group. Changes in the shape of crania and mandibles were displayed visually using the wireframe outlines in MorphoJ compared to the average shape along each canonical variate with the outlines at the extremes. Lastly, we calculated the centroid size (the square root of the sum of squared Euclidean distances from each landmark to their own centroid) for the cranium and mandible. Generally, this measurement is the standard for assessing size in GMA (Le and Kume, 2000;Schwarze et al., 2019).

Statistical Analysis
Descriptive statistics were calculated for linear measurements and centroid sizes. For significance tests, we used one-way ANOVA. Graphpad Prism (Version 8.0.2) was used to perform tests and graphics. To determine whether shape differences were statistically significant, values of p were also calculated using Procrustes ANOVA and Mahalanobis distances (Supplementary Table S1). To determine inter-observer landmark reproducibility, two observers (MGH and HK) independently located the landmarks on 10 randomly selected samples (Supplementary Figure S2). A 10,000 round permutation test was performed on the Procrustes distance between the observers' landmarked samples, testing for mean overall shape differences between them. For linear measurements, intra-observer reliability was analyzed using the paired-t test and Bland-Altman method (Supplementary Figure S3). The main observer (MGH) took a second set of measurements (seven randomly selected samples) 1 week after the first set to confirm the absence of intra-observer error.

Effects of SD on Cranial Morphology
A one-way ANOVA was conducted to illustrate the effect of SD on craniofacial morphology. There was a significant difference  an, anterior; di, distal; la, lateral; me, medial; po, posterior; pr, proximal. in mean centroid size [F(3,38) = 14.25, p < 0.0001] between the four groups. Post hoc comparisons using the Tukey's test were carried out. Switching diets for a single generation significantly decreased cranial size in F1SD and F15SD-F1HD mice (Figure 2A; Table 2). Notably, there were no significant differences in the cranial measurements between HD Con and F15SD mice (Figure 2A; Table 2).
The shape of the mouse crania differed significantly (Procrustes ANOVA, F = 14.72, p < 0.0001; Supplementary Table S1). When the shape of the cranium was considered (i.e., without the effects of centroid size), the majority of the variation was accounted for by PC1 and PC2 (17.9 and 15.4%, respectively, Figure 2B, Supplementary Figure S4). F15SD mice possessed a shorter midface were wider transversely due to outwardly displaced zygomatic arches and taller vertically due to raised posterior cranial vaults (Figures 2B,C). Conversely, HD Con mice were associated with forward displacement of the midface and were narrower transversely due to inwardly displaced zygomatic arches and shorter vertically due to depressed posterior cranial vaults (Figures 2B,C).
CVA revealed aspects of cranial shape variation capable of discriminating among the four groups, following adjustment for centroid size (Figure 2D). The first and second canonical variates (CV1 and CV2) accounted for 96% of the variation (data not shown). CV1 accounted for more than 75% of the observed shape variation. CV1 positive scores were associated with smaller crania, depressed posterior cranial vaults, anterior displaced midface, and narrowed transverse widths due to inwardly displaced zygomatic arches (Figure 2E). Accounting for 20% of the shape variation, CV2 separated the four groups as well, with HD Con mice at one end, F15SD mice at the other, and F1SD and F15SD-F1HD mice in between ( Figure 2D). Clusters with positive CV2 values such as F15SD mice possessed wider zygomatic arches and shortened midface in the anteroposterior dimension ( Figure 2E). Conversely, HD Con mice were clustered at negative CV2 values and possessed narrowed zygomatic arches, raised posterior cranial vaults, and longer alveolar ridge relative to the average shape ( Figure 2E). F1SD and F15SD-F1HD mice were clustered close together in between HD Con and F15SD mice ( Figure 2D).

Effects of SD on Mandibular Morphology
Similar to the cranium, there were statistically significant differences in the mean centroid size [F(3,38) = 16.76, p < 0.0001] between groups as determined by one-way ANOVA. The biggest difference in size was again observed with switching of diets in a single generation, F1SD and F15SD-F1HD mice. A Tukey's post hoc test revealed that the mandible centroid size, as well as linear measurements of F1SD and F15SD-F1HD mandibles were significantly lower than HD Con and F15SD mandibles ( Figure 3A; Table 2). There were no significant differences in BGW between the groups but the IMW was significantly lower in F1SD, F15SD, and F15SD-F1HD compared to HD Con mandibles.
The shape of the mouse mandibles differed significantly (Procrustes ANOVA, F = 16.76, p < 0.0001; Supplementary Table S1). PCA of allometry-free shape data showed clear separation of HD Con and F15SD mandibles. The first two PCs of the PCA accounted for 45% of the shape variance ( Figure 3B, Supplementary Figure S4). There were differences in mandible shape that discriminated between the groups: F1SD and F15SD-F1HD mice possessed longer mandibular bodies, wider and shorter mandibular rami, posteriorly displaced coronoid processes, and shallow mandibular notches ( Figure 3C). Notably, F15SD mice had shorter mandibular bodies, narrower and elongated mandibular rami, anteriorly displaced coronoid processes, and a deep mandibular notch ( Figure 3C).
CVA revealed distinct mandibular shape variation among the four different groups (Figure 3D). CV1 and CV2 were associated with 91% of the variation (data not shown). Clusters with lower values of CV1, such as F1SD and F15SD-F1HD mice, presented with longer mandibles and smaller transverse dimensions due to inwardly displaced mandibular condyles and gonial angles ( Figure 3E). Furthermore, clusters with positive values of CV2, such as F15SD mice, possessed shorter mandibles and elongated and narrower rami ( Figure 3E).

Effects of SD on Condyle Morphology
Using semi-landmarks and GMA to analyze condyle morphology, multivariate regression of Procrustes coordinates on centroid size showed that 12.3% of shape variation was due to static allometry in condyle morphometry between the four groups. The shape of the mouse condyles differed significantly (Procrustes ANOVA, F = 2.98, p < 0.0001; Supplementary Table S1). PC1 and PC2 were associated with 42% of shape variation in the mandibular condyle of all specimens (Figure 4A; Supplementary Figure S4). Clusters with higher values of PC1 represented F15SD mice with wider condylar necks in the anteroposterior dimension and outwardly displaced condylar heads, in comparison to HD Con mice that were located at the negative end values of PC1 ( Figure 4B). F1SD and F15SD-F1HD mice were clustered in between HD Con and F15SD mice. CVA was consistent with PCA outcomes and clearly distinguished the four groups ( Figure 4C). CV1 and CV2 were associated with 86% of the variation (data not shown).
Positive CV1 values such as that observed for control mice possessed inwardly displaced condyles ( Figure 4D). Clusters with negative CV2 values such as in F15SD mice presented with wider condylar necks and outwardly displaced condylar heads ( Figure 4D).

DISCUSSION
Our study, to our knowledge, is the first to investigate changes in mouse craniofacial morphology due to prolonged, multigenerational exposure to SD using 3D GMA. Prior studies using animal models have analyzed craniofacial adaptation to SD and load reduction in a single generation using 2D linear measurements (Corruccini and Beecher, 1982;Kiliaridis et al., 1985;Dias et al., 2011;Denes et al., 2013Denes et al., , 2018Anderson et al., 2014;Goto et al., 2015;Kono et al., 2017). Mice are preferred models because of the use of isogenic strains, tight environmental control of multiple specimens, relatively short gestation periods, and extensive information from prior genetic and environmental studies (Fish, 2016;Liu, 2016;Van Otterloo et al., 2016).

Short-Term Effects of SD on Craniofacial Morphology
Our results showed that short-term SD (F1SD mice) significantly decreased craniofacial size (Figures 2, 3). These results are consistent with studies analyzing the effects of low masticatory function on craniofacial structures in one generation using 2D linear measurements (Kiliaridis, 2006;Odman et al., 2008;Guerreiro et al., 2013;Hichijo et al., 2015;Fujita and Maki, 2018;Ödman et al., 2019). Interestingly, one study demonstrated that SD has no effect on the rat cranial length and width despite using the same linear measurements (Katsaros et al., 2002). With respect to size-free shape changes, our results were partially consistent with prior studies (Kiliaridis et al., 1985;Kono et al., 2017). Our F1SD mice possessed transverse widening of zygomatic/temporal regions, flattening of the neurocrania, less convex angular processes, shallow antegonial notches, inward tipping of the condylar processes, and lengthening of the mandibular bodies. Prior studies largely showed flattening of the neurocrania and less convex angular processes but not our other findings (Kiliaridis et al., 1985;Kono et al., 2017). The differences in our results and prior studies may be due to differences in experimental design. In prior studies such as by Katsaros et al. (2002), rats were given SD at weaning (28 days) for a 4-week duration. In our study, F1SD mice were exposed to SD beginning at E0.5 (parents were given on SD) until 6 weeks old. Although F1SD mice did not directly eat SD from E0.5 to ~21 days of age, they developed in cages with SD and were indirectly exposed to SD through the parents, and it is unclear the exact age that F1SD would begin to eat SD. Thus, it is interesting to posit whether indirect exposure to SD via parents or early exposure to SD may influence offspring morphology. The condyles and condylar processes of F1SD mice showed decreases in the anteroposterior dimension, findings that are in agreement with prior studies in which masticatory hypofunction decreased condylar dimensions (Kiliaridis et al., 1999;Scheidegger et al., 2018). It should be noted that condylar heights were not analyzed in our study because the lower boundary of the semi-landmarks was based on an imaginary line passing through the tip of the coronoid process and the angle of the mandible, which are anatomical sites prone to functional load adaptations.

Long-Term Effects of SD on Craniofacial Morphology
Our experimental design exposed mice to 15 continuous generations of SD to observe long-term morphological effects. Our design would eliminate the confounding influence of transient weight loss presumably due to adaptation to a new diet, which has been noted upon switching to SD (Kiliaridis et al., 1985). The fact that no differences were recorded in mean cranium and mandible centroid sizes nor in the linear measurements of the craniofacial structures between HD Con and F15SD mice demonstrated that long-term SD had no impact on skull size. Moreover, we did not observe any differences  F1SD and F15SD-F1HD mandible centroid sizes and linear measurements were smaller than HD Con and F15SD. Post hoc comparisons using the Tukey test with one-way ANOVA (A). Data represents the means ± SD. (B) PCA of mandibular shape: HD Con (red), F1SD (brown), F15SD (blue), and F15SD-F1HD (green). The ellipses represent the confidence range of the means of each group. (C) Wireframe illustrations demonstrate shape changes associated with the extreme positive and negative PC1 and PC2 values in shape (red line) and average shapes (dotted black line). (D) CVA shows clear separation between HD Con and F15SD mice confirming shape differences between the two groups. (E) Wireframe illustrations demonstrate shape changes associated with the extreme positive and negative CV1 and CV2 values in shape (red line) and average shapes (dotted black line). *p < 0.05; **p < 0.005; ***p < 0.0005; ****p < 0.0001.
in weights between HD Con and F15SD male and female mice (Supplementary Figure S5). Our experimental design is the likely reason that changes in craniofacial morphology in F15SD mice differ from all prior single generation studies. F15SD mice possessed oval neurocrania with depressed posterior cranial vaults, outward displaced zygomatic arches, and posteriorly displaced maxillary alveolar bones. F15SD mandibles were shortened in the anteroposterior dimension and were associated with increases in the transverse dimension, decreases in mandibular alveolar bone heights, and accentuated anti-gonial notches. F15SD condyles presented with increases in the anteroposterior dimension and outward displacement, which were coordinated with transverse widening of zygomatic/temporal bones and mandibles. The lateral displacement of the zygomaticotemporal region was coordinated with lateral displacement of the mandibular posterior region and condyle. Interestingly, the changes in F15SD skull morphology were similar in one prior study in pigs fed SD (Larsson et al., 2005). The pigs were fed SD at weaning (5 weeks old) until 22 months of age, the longest duration of time that an animal had been fed SD in a single generation study. The pigs showed significantly wider arches and fewer posterior crossbites compared to HD-fed animals, similar to that observed in our F15SD mice.

CONCLUSION
To our knowledge, our study is the first to analyze mice on prolonged, multi-generational SD to determine the effects of diet consistency on craniofacial form. Our results suggest that changes in craniofacial shape and size in response to decreased functional loads may be independent. Short-term SD resulted in smaller cranial and mandibular sizes likely due to transitional adaptation of mice to new diets. Prolonged SD did not affect the size of the cranium and mandible but altered the shape. Shape changes associated with prolonged SD consumption included widening in the transverse dimension, as well as decreases in the anteroposterior dimension of both the cranium and mandible. The lateral displacement of the zygomatico-temporal region of the cranium was coordinated with lateral displacement of the posterior mandible including the condyle.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ Supplementary Material and upon request.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee (IACUC) at UCSF and performed under animal research protocol number AN164201.

AUTHOR CONTRIBUTIONS
MH and AJ designed the study. BZ maintained the mice. MH and HK performed the landmarking. TC performed the microcomputed tomography. MH conducted the geometric morphometric analysis. NY reviewed the analysis. MH and AJ wrote the manuscript. All authors contributed to the article and approved the submitted version.