Morgane Dubied
email: morgane.dubied@u-bourgogne.fr

Sophie Montuire
email: sophie.montuire@ephe.psl.eu

Nicolas Navarro
email: nicolas.navarro@ephe.psl.eu

Orcid

Commonalities and evolutionary divergences of mandible shape ontogenies in rodents Running title: Evolution of postnatal ontogenies in rodents

Keywords: Rodents, geometric morphometrics, macroevolution, mandible shape, postnatal growth

In mammals, significant changes take place during postnatal growth, linked to changes in diet (from sucking to gnawing). During this period, mandible development is highly interconnected with muscle growth and the epigenetic interactions between muscle and bone control the spatialization of bone formation and remodeling in response to biomechanical strain. This mechanism contributes to postnatal developmental plasticity, and may have influenced the course of evolutionary divergences between species and clades. We sought to model postnatal changes at a macroevolutionary scale by analyzing ontogenetic trajectories of mandible shape across 16 species belonging mainly to two suborders of Rodents, Myomorpha and Hystricomorpha, which differ in muscle attachments, tooth growth, and life-history traits. Myomorpha species present a much stronger magnitude of changes over a shorter growth period. Among Hystricomorpha, part of the observed adult shape is set up prenatally, and most postnatal trajectories are genus-specific, which agrees with non-linear developmental trajectories over longer gestational periods. Beside divergence at large scale, we find some collinearities between evolutionary and developmental trajectories. A common developmental trend was also observed, leading to enlargement of the masseter fossa during postnatal growth. The tooth growth, especially hypselodonty, seems to be a major driver of divergences of postnatal trajectories. These muscle-and tooth-related effects on postnatal trajectories suggest opportunities for developmental plasticity in the evolution of the mandible shape, opportunities that may have differed across Rodent clades.

Introduction

Understanding developmental mechanisms in evolution is crucial to apprehend the diversification of organismal forms [START_REF] Darwin | On the origin of species[END_REF], Alberch 1980[START_REF] Smith | Developmental constraints and evolution: a perspective from the Mountain Lake conference on development and evolution[END_REF]. It is well recognized nowadays that these mechanisms lead to discontinuities and directionalities in the space of shapes (Alberch 1980[START_REF] Gerber | Not all roads can be taken: development induces anisotropic accessibility in morphospace[END_REF][START_REF] Salazar-Ciudad | Why call it developmental bias when it is just development?[END_REF]. Indeed, they induce biases in the phenotypic outcomes of random mutations [START_REF] Alberch | Developmental constraints in evolutionary processes[END_REF], Hallgrimsson et al. 2006[START_REF] Uller | Developmental plasticity and evolutionary explanations[END_REF] and therefore in the evolutionary trajectories [START_REF] Uller | Developmental bias and evolution: A regulatory network perspective[END_REF][START_REF] Kavanagh | Developmental bias in the evolution of phalanges[END_REF]. The existence of such developmental bias results in some parallelism between developmental and evolutionary trajectories [START_REF] Alberch | Size and shape in ontogeny and phylogeny[END_REF][START_REF] Gould | Ontogeny and Phylogeny[END_REF], Webster and Zeldicth 2005). Despite the abundant theoretical literature about the role of development in evolution [START_REF] Arthur | Developmental drive: an important determinant of the direction of phenotypic evolution[END_REF][START_REF] Uller | Developmental plasticity and evolutionary explanations[END_REF], there are practical problems to gather ontogenetic data at a large scale. Empirical works are still needed to assess whether directions in the morphospace supported by developmental variation are also the most evolvable. Developmental processes influence thus the phenotypic production and the availability of variation to selection, questioning the importance of evolution of the ontogenetic trajectories themselves in shaping the diversity of organismal shapes.

The evolution of ontogenetic trajectories involves a complex interplay of the rate and timing of development of the various parts constituting an organism as well as with the internal and/or external environment. During pre-and postnatal development, organismal growth responds to many genetic, biomechanical, and environmental factors. As organs grow, they compete for space and resources, while maintaining a functional phenotype under a variety of selection regimes [START_REF] Dibner | Metabolic challenges and early bone development[END_REF][START_REF] Nijhout | Competition among body parts in the development and evolution of insect morphology[END_REF][START_REF] Olsen | Bone development[END_REF]. In vertebrates, for example, epigenetic interactions will compensate for, and coordinate, the growth of the organs that form the head, to acquire and/or maintain functions such as occlusion between the lower and upper jaws (Hallgrimsson andHall 2011, Lieberman 2011). Similarly, changes in muscle forces during development regulate the spatialization and intensity of bone remodeling [START_REF] Herring | Muscle-bone interactions and the development of skeletal phenotype: jaw muscles and the skull[END_REF][START_REF] Zelditch | Epigenetic interactions: the developmental route to functional integration[END_REF][START_REF] Zelditch | Modularity of the rodent mandible: Integrating bones, muscles, and teeth[END_REF], which may change not only the structure but also the shape of the bone (Martine-Maza et al. 2016). Organ interactions may drive the occurrence and intensity of developmental plasticity and the indirect response to selection [START_REF] Fusco | Evolving Pathways: Key themes in evolutionary developmental biology[END_REF]. By responding to developmental changes in forces and movements, epigenetic interactions between bones and other tissues could drive morphological variation (Hallgrimsson and Hall 2011). Epigenetics could therefore structure evolutionary trajectories, potentially biasing diversification among taxa (Young andBadyaev 2007, Renvoisé et al. 2017) and thus becoming a major determinant of plasticity-led evolution [START_REF] West-Eberhard | Developmental Plasticity and Evolution[END_REF][START_REF] Uller | Developmental plasticity and evolutionary explanations[END_REF].

Postnatal growth is a key period in mammals, during which skull variation is strongly modified [START_REF] Mitteroecker | The Ontogenetic trajectory of the phenotypic covariance matrix, with examples from craniofacial shape in rats and humans[END_REF]) and then stabilized [START_REF] Zelditch | The ontogenetic dynamics of shape disparity[END_REF]. During this transition from the juvenile to the adult head, major biomechanical changes occur in response to changes in the use of the masticatory apparatus related to weaning, from sucking to gnawing and chewing movements [START_REF] Curley | The Meaning of weaning: influence of the weaning period on behavioral development in mice[END_REF]). An important component of adult anatomy is the eruption and growth of teeth, in relation to the acquisition of mastication. Tooth development partly determines the complexity of the ontogenetic trajectory of the mouse mandible [START_REF] Swiderski | The complex ontogenetic trajectory of mandibular shape in a laboratory mouse[END_REF]. The ever-growing incisor leads to bone remodeling during mandible growth [START_REF] Renvoisé | Developmental mechanisms in the evolution of phenotypic traits in rodent teeth[END_REF]. This simple, highly integrated, skeleton element strongly responds to changes in mechanical constraints related to food hardness [START_REF] Anderson | Adaptive plasticity in the mouse mandible[END_REF][START_REF] Menegaz | Ontogenetic and functional modularity in the rodent mandible[END_REF], or to muscle dystrophy [START_REF] Renaud | Epigenetic effects on the mouse mandible: common features and discrepancies in remodeling due to muscular dystrophy and response to food consistency[END_REF]. The mandible is an ideal model system to observe organ interactions, and thus to test their evolutionary consequences, as its development is tightly integrated with that of teeth and muscles (Atchley 1993, Klingenberg and[START_REF] Klingenberg | Development of the mouse mandible: a model system for complex morphological structures[END_REF].

Rodentia is a very diverse and disparate mammal order, in which postnatal changes can be observed on a large scale. The muscle attachment to the cranium is used to classify rodents into three main suborders: Myomorpha, Sciuromorpha, and Hystricomorpha [START_REF] Hautier | Masticatory muscle architecture in the gundi Ctenodactylus vali (Mammalia, Rodentia)[END_REF][START_REF] Hautier | Evolution of the zygomasseteric construction in Rodentia, as revealed by a geometric morphometric analysis of the mandible of Graphiurus (Rodentia, Gliridae)[END_REF][START_REF] Simpson | The principles of classification and a classification of Mammals[END_REF][START_REF] Wood | Grades and clades among rodents[END_REF]. The mandible presents a wide variety of forms across these suborders (Hautier et al. 2011), which group into two mandibular structures [START_REF] Tullberg | Ueber das System der Nagethiere: eine phylogenetische Studie[END_REF]: hystricognaths (with a strong inwardly projecting angular process and an outwardly projecting condyloid process) or sciurognaths (with an outwardly projecting angular process and an inwardly projecting condyloid process). The insertions of the masseter muscles to the mandible and the cranium were used by early authors [START_REF] Brandt | Beiträge zur nähern Kenntniss der Säugethiere Russland[END_REF][START_REF] Wood | Grades and clades among rodents[END_REF] to describe four different morphologies across rodents: protrogomorphy (small to medium infraorbital foramen, a reduced masseter and a predominant temporalis), sciuromorphy (very small infraorbital foramen and a lateral masseter inserted in the anterior part of the orbital arch), hystricomorphy (infraorbital foramen is very developed and almost circular) and myomorphy (infraorbital foramen is narrow at the base and widened in its upper part). Rodents present also three main types of molar growth [START_REF] Renvoisé | Developmental mechanisms in the evolution of phenotypic traits in rodent teeth[END_REF]: root apex closed after tooth maturation (brachyodonty), ever-growing with open root (hypsodonty), or with no root (hypselodonty); and molars take up a great deal of space in the mandible as they develop, particularly in hypselodont species.

Beyond intricate interconnections with postnatal developing structures (muscles and tooth), shape disparity of the mandible seems to be generated prenatally, at least in Sciuridae, with compensated postnatal ontogenetic variations [START_REF] Zelditch | Disparate postnatal ontogenies donot add to the shape disparity of infants[END_REF], suggesting a relatively reduced role for biomechanical plasticity in the establishment of species differences. Data supporting this result include most species belonging to the Sciuridae, which is relatively conservative in terms of skull shape (Cardini and Thoringhton 2006). Myomorpha and Hystricomorpha are diversified in shape, probably because of their high diversification rate and high diversity [START_REF] Alhajeri | Disparity and evolutionary rate do not explain diversity patterns in Muroid rodents (Rodentia: Muroidea)[END_REF][START_REF] Wilson | Allometric disparity in rodent evolution[END_REF]. Therefore, studying the evolution of their ontogenies will likely provide additional contrasts to better understand the importance of postnatal growth in the set-up of species divergences.

These two suborders and the evolution of their ontogenies will be the focus of this study, mainly because they exhibit diverse life-history traits among species (Wilson et al. 2016;[START_REF] Wolff | Rodent societies: an ecological and evolutionary perspective[END_REF]. Gestation and weaning are shorter in Myomorpha, which generally give birth to large litters many times a year. This behavior involves a short in-utero development and a rapid and therefore mechanically brutal weaning [START_REF] Curley | The Meaning of weaning: influence of the weaning period on behavioral development in mice[END_REF]. Hystricomorpha give birth to small litters a few times a year. Parental care is therefore different, with longer behavioral weaning and less brutal dietary weaning, over a much longer period [START_REF] Wolff | Rodent societies: an ecological and evolutionary perspective[END_REF]. The potential for postnatal plasticity could likely be different between the two suborders. A wide variety of diets is also found within these two groups [START_REF] Nowak | Walker's Mammals of the World[END_REF]. Changes in diet might correlate with evolutionary changes in muscle attachment, direction of mastication, and mandible shape (Alvarez and Pérez 2019).

Geometric morphometrics is used to describe postnatal trajectories of the mandible shape on specifically collected species with identified juveniles (before weaning), instead of using static allometry from a series of specimens varying in size. Comparing these postnatal growths within and between Myomorpha and Hystricomorpha suborders-together with two Sciuromorpha species for comparison-will provide valuable data about the influences of development in evolution, for instance if evolutionary and developmental directions in the morphospace coincide [START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF][START_REF] Gerber | Not all roads can be taken: development induces anisotropic accessibility in morphospace[END_REF][START_REF] Wilson | Allometric disparity in rodent evolution[END_REF][START_REF] Zelditch | Disparate postnatal ontogenies donot add to the shape disparity of infants[END_REF]. It will also provide data about the evolution of ontogenesis in rodents, if these evolutionary changes sustain species divergences, and if postnatal growth is a key period for the establishment of adult disparity. Contrasting the variation in these mandible shape trajectories to life-history traits and to the development of muscles and teeth will underline the importance of the development of these traits in the evolution of the mandible ontogenesis.

Materials and methods

SPECIMENS

Specimens were selected from five different collections: the Royal Belgian Institute of Natural Sciences (IRSNB), the Natural History Museum Geneva (MHNG), the National Museum of Natural History Paris (MNHN), the Natural History Museum Basel (NMB) and the Biogéosciences Laboratory, Dijon (Table S1). Mesocricetus auratus (RjHan:AURA), Meriones unguiculatus (Crl:MON) and Mus musculus (BALB/c) are laboratory individuals raised at the University of Burgundy (Project APAFIS#18405-2019011014262528) for which the age of each individual is known (7 days for juveniles). Juvenile specimens from Museum collections were carefully checked for their putative age with these following criteria about tooth eruption and ossification, which correspond to the observed characteristics made on laboratory 7 day-old juveniles. Molars' occlusal surface is just emerged from the bone and incisors are barely emerged as well. The ossification of mandibular processes is poor (presence of cartilage at the top of the postcondylar process) and the one of the cranium is incomplete (cartilage between the bones constituting the cranium). It results that a large number of registered juveniles were disregarded to keep only undoubtful early juveniles. The final data set for this study includes 105 hemimandibles attributed to 16 species: five for Myomorpha, eight for Hystricomorpha, two for Sciuromorpha, and Oryctolagus cuniculus (Lagomorpha) as the outgroup species (Fig. 1). Only one hemimandible per specimen was used, chosen by state of preservation. For the purpose of checking the effect of the age of juveniles, 27 additional specimens were also sampled and added as supplementary data. They represented 15 days old juveniles of Myomorpha species raised in the lab and a few juveniles of Microtus arvalis older than one week.

Mandibles shorter than 5 cm in length were scanned by µCT-scan (Bruker Skyscan 1174) and reconstructed using Avizo®9.2 (FEI systems). Those longer than 5 cm were scanned with a Shining 3D® Einscan Pro surface scanner. One specimen was processed on the two acquisition pipelines and digitized several times. Preliminary analysis identified no spurious variation in these replicates. Before landmark digitization, 3D models were decimated to 200 000 faces, using the Rvcg R package 0.18 [START_REF] Schlager | Chapter 9 -Morpho and Rvcg -Shape analysis in R: R-Packages for geometric morphometrics, shape analysis and surface manipulations[END_REF].

LIFE HISTORY

Data about dental development [START_REF] Ungar | Mammal Teeth: Origin, Evolution, and Diversity[END_REF], and ecology such as reproductive traits (Wilson et al. 2016;De Magalhaes and Costa 2009), were gathered from the literature. In summary (Fig. 1), gestational duration varies considerably, three to nine times shorter for Myomorpha (19 to 22 days, median = 20), than for Hystricomorpha (66 to 210 days, median = 125). A similar difference is also observed for weaning time, which occurs in Myomorpha between 20 and 26 days after birth (median = 22), while it occurs from 18 to 100 days after birth (median = 60) in Hystricomorpha. The two suborders have a similar number of litters per year (a maximum of five), but on average Myomorpha have four litters per year, while Hystricomorpha have only two. About the skull structure, all species belonging to Myomorpha and Hystricomorpha are either myomorphous or hystricomorphous respectively. Sciuromorpha species are either sciuromorphous (Funisciurus lemniscratus) or myomorphous (Muscardinus avellanarius). All species are sciurognaths except the Hystricomorpha species, which are hystricognaths.

LANDMARKS AND SEMILANDMARKS

Ten landmarks were digitized, together with 33 semilandmarks divided into three curves (with 4 on the coronoid process, 9 along the lunar notch, and 10 between postcondylar and angular processes; Fig. 2) using the Digit3DLand R package 0.1.3 [START_REF] Laffont | Digitalization of 3D landmarks on mesh[END_REF]. Semilandmarks were slid along their tangent by minimizing the bending energy, and back-projected on the 3D curves [START_REF] Gunz | Semilandmarks in three dimensions[END_REF].

DATA ANALYSIS Procrustes alignment, checked for approximation and imbalance effects

Landmarks and sliding semilandmarks were superimposed using a full generalized Procrustes analysis (GPA) with the Morpho R package 2.8 [START_REF] Schlager | Chapter 9 -Morpho and Rvcg -Shape analysis in R: R-Packages for geometric morphometrics, shape analysis and surface manipulations[END_REF]. Tangent space approximation of distances between specimens was checked against Riemannian distance, which has been shown to be generally accurate across studies encompassing a variety of taxa and taxonomic ranges (Klingenberg 2020), including mammalian orders [START_REF] Marcus | Application of landmark morphometrics to skulls representing the orders of living mammals[END_REF]. The reference shape on which landmark conformations are aligned is known to structure variation in the tangent space, thus determining any distortion of distances between individuals (Bookstein 2016). As developmental stages and species were not equally sampled or equally represented between clades, we also checked whether aligning landmark conformations on a weighted average rather than on the grand mean (the usual method) still provides consistent results. The average from each developmental stage of each clade (weighting species in the clade equally) was therefore used as reference shape in a new Procrustes superimposition. The angle between this weighted mean shape and the original mean shape was computed and compared to the angles observed between any two specimens. The Lagomorphs were aligned only as supplementary data.

Developmental, phylogenetic-aligned and postnatal morphospaces

A developmental morphospace, the embedding of both juveniles and adults in the same morphospace [START_REF] Eble | Evolutionary dynamics: exploring the interplay of selection, accident, neutrality, and function[END_REF], was obtained by applying principal component analysis (PCA) to the covariance matrix of the tangent coordinates. This morphospace maximizes shape variation in directions that describe either developmental changes, evolutionary changes or both simultaneously. To compare the developmental structuration of the morphospace and evolutionary divergences, a phylogenetically aligned morphospace was also computed based on the principal component ordination of the phylogenetic covariance matrix ! (Collyer and Adams 2021). The phylogenetic covariance matrix ! was computed according to a Brownian-motion model of evolutionary divergence. This space (PaPCs) maximizes variation in directions that describe evolutionary divergences. This analysis was done with the geomorph R package 3.2.1 [START_REF] Adams | geomorph: Software for geometric morphometric analyses[END_REF]. To observe the structuration of the variation of the postnatal trajectories alone, a PCA was performed on the covariance matrix of the vectors " of differences between juvenile and adult shape averages # of each species. These vectors are the average shape changes occurring during postnatal growth. Scaling these vectors by the difference between juvenile and adult average sizes would have led to a rough estimate of allometric vectors. As Hystricomorpha species are much larger than other species, rescaling would have resulted in trivial size differences.

A tree pruned for the species present in the sample (Wilson et al. 2016) was built with the ape R package 5.3 [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF], and was projected onto the different morphospaces, using the phytools R package 0.7-20 [START_REF] Revell | phytools: an R package for phylogenetic comparative biology (and other things): phytools: R package[END_REF], to observe their phylogenetic structuration.

Ancestral postnatal trajectories, common pattern of growth and evolutionary changes

The phylogenetic mean of postnatal trajectories was computed according to the generalized least-squares estimator " $ !"#$ = ('′! %& ') %& '′! %& " [START_REF] Rohlf | Comparative methods for the analysis of continuous variables: geometric interpretations[END_REF][START_REF] Revell | Size-correction and principal components for interspecific comparative studies[END_REF], with " the vectors of shape changes between juveniles and adults and ! the phylogenetic covariance matrix. This estimate corresponds to the value at the root of the tree [START_REF] Rohlf | Comparative methods for the analysis of continuous variables: geometric interpretations[END_REF]. The " vectors correspond to postnatal shape changes, " $ !"#$ estimates the ancestral pattern of shape changes during postnatal growth according to Brownian motion. Estimate of the postnatal trajectory for each suborder was computed as the least-squares (LS) means " + = ,-, with , the matrix corresponding to the linear contrast for the suborder, and -= (.′! %& .) %& . ' ! %& ", the phylogenetic generalized least-squares estimates of suborder differences.

To check the alignment of the principal directions of the morphospace (developmental or phylogenetically-aligned) with the phylogenetic mean postnatal trajectory " $ !"#$ , or with the suborder estimates of postnatal trajectories " +, angles between these estimates and the shape features described by PCs were estimated as 0 = 123 %& (4 ( • 4 ) ) [START_REF] Klingenberg | Heterochrony and allometry: the analysis of evolutionary change in ontogeny[END_REF][START_REF] Zelditch | Chapter 14 -Morphometrics and systematics[END_REF]), with the normalized vectors 4 being either eigenvectors or normalized trajectories to unit length. The probability that the angle between two random vectors is lower than the observed angle was estimated as the area of the cap of a hypersphere defined by this angle [START_REF] Li | Concise formulas for the area and volume of a hyperspherical cap[END_REF], and computed with the Morpho R package 2.8 [START_REF] Schlager | Chapter 9 -Morpho and Rvcg -Shape analysis in R: R-Packages for geometric morphometrics, shape analysis and surface manipulations[END_REF]. Approximation of standard errors on angles based on " + were computed based on the sampling distribution of - (Meyer and[START_REF] Meyer | Sampling based approximation of confidence intervals for functions of genetic covariance matrices[END_REF][START_REF] Houle | Estimating sampling error of evolutionary statistics based on genetic covariance matrices using maximum likelihood[END_REF], and according to " * ~ 7(" +, 8 + ⊗ ,(.′! %& .) %& ,′), where 8 + is the residual covariance matrix, and ⊗ stands for the Kronecker product. Descriptive statistics on sampled angles were computed with the circular R package 0.4-93 [START_REF] Agostinelli | R package 'circular': Circular Statistics. R Package. Alberch, P[END_REF].

The angle between phylogenetically aligned variation (PaPCA) and the average directions of postnatal changes (" $ !"#$ and " +) was measured as : ( = 123 %& {(4 ( , ,, -4 ( ) ..0 } with , 12 × 5 the matrix of = eigenvectors from the PaPCA, where = was set equal to half the divergence dimensionality (= = 7) and 4 the normalized postnatal vectors [START_REF] Krzanowski | Between-groups comparison of principal components[END_REF]). To assess more formally the degree of independence between developmental trajectories and evolutionary divergences, a common spectral decomposition of the two matrices was performed [START_REF] Krzanowski | Between-groups comparison of principal components[END_REF]). This approach captures the overall similarity from the eigen decomposition of

? = ∑ , 6 , 6 - ! (7&
where , is the matrix of the q first eigenvectors of either the divergences or the postnatal changes. The eigenvalues close to the maximum of p (here p = 2) means that the evolutionary or postnatal changes could be inferred from a linear combination of the eigenvectors of the other. Similarly to above, an angle : ( could be defined to estimate how similar an eigenvector of ? is to the evolutionary or postnatal changes [START_REF] Krzanowski | Between-groups comparison of principal components[END_REF].

Comparison of juvenile and adult disparities

Even if some evolutionary divergences arise in the directions of postnatal growth, adult disparity may or not increase compared to juvenile disparity as changes in the amount or direction could be compensated and finally most of the disparity could be established prenatally as suggested for Sciuridae [START_REF] Zelditch | Disparate postnatal ontogenies donot add to the shape disparity of infants[END_REF]. Disparity levels between developmental stages of Myomorpha and Hystricomorpha were calculated as the Procrustes variance, i.e., the sum of squared Euclidean distances in the tangent space between species means # and their average divided by the number of species in the group (Zelditch et al. 2004, Drake and[START_REF] Drake | Large-scale diversification of skull shape in domestic dogs: Disparity and modularity[END_REF]. Bootstrapping of species means was used to estimate standard errors on disparity [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF][START_REF] Foote | Morphological disparity in Ordovician-Devonian crinoids and the early saturation of morphological space[END_REF][START_REF] Navarro | MDA: a MATLAB-based program for morphospace-disparity analysis[END_REF]). Disparities were also computed on a sequence including successive principal components (PCs) defining the developmental morphospace. This sequence was used to show how taxa spread over the morphospace as defined by the main patterns of shape differentiation among species and developmental stages. Comparison of disparity between the two suborders also assessed whether evolvability varies across anatomies and modes of development.

Evolution of postnatal trajectories

The degree of phylogenetic signal in the postnatal vectors " was assessed relative to their expectation given a Brownian evolution, using the multivariate version of the K statistics (Adams 2014), with 10 000 simulations. The Kmult was computed with the geomorph R package 3.2.1 [START_REF] Adams | geomorph: Software for geometric morphometric analyses[END_REF]. This package also returns effect size as the z-score standardization of the Kmult statistics, given its distribution under the null obtained from simulation. For comparison, the phylogenetic signal was also computed separately for juvenile and for adult shapes.

To assess whether postnatal trajectories have evolved between suborders, the differences between suborder estimates " + with regard to their magnitude or direction of postnatal shape changes were computed. Standard errors on angles were approximated using " * sampling as explained above. Evolutionary divergence in the magnitude of shape changes occurring after birth between the Myomorpha and Hystricomorpha was evaluated from the LS-means contrast of the Euclidean length of the suborder trajectories A" + 89: A -A" + ;9$ A, with a t test using approximate standard error based on " * sampling and the residual degree of freedom of the linear model.

Approximately phylogeny-corrected postnatal changes were computed as " < = C %& " (Arnold 1981, Schabenberger and[START_REF] Schabenberger | Statistical methods for spatial data analysis[END_REF], where " is the vector of shape changes between juveniles and adults of each species and C %& is the inverse square root of the phylogenetic covariance matrix !. To obtain C %& , the eigen decomposition of ! was used as C %& = DE %..0 D′, with D the matrix of eigenvectors, and E the diagonal matrix of eigenvalues of !. Their angles with the ancestral growth " $ !"#$ was measured to assess the degree of divergences from the root of the tree. The pairwise angles between these trajectories within suborders were computed to assess the remaining variation between species once Brownian divergence is factored in, variation that could have resulted from evolutionary divergence from other traits such as life history traits. This question was evaluated more formally using a phylogenetic generalized least-squares analysis (pGLS) with the effects of tooth growth and reproductive traits on postnatal shape changes ". Considering the number of species sampled and differentiation between and within Myomorpha and Hystricomorpha, traits could be strongly aggregated [START_REF] Adams | Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations[END_REF]. Correlations between phylogeny and life history traits were first checked using two-block partial least-squares [START_REF] Rohlf | Use of two-block partial least-squares to study covariation in shape[END_REF], following the approach of [START_REF] Adams | Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations[END_REF]. Because of the strong clustering of gestational and weaning duration and their resulting high correlation with phylogeny, they were pooled within suborders before being analyzed. Significance of the Procrustes sum of squares was evaluated using 10 000 residual permutations with the RRPP R package 0.5.2 [START_REF] Collyer | RRPP: An r package for fitting linear models to high-dimensional data using residual randomization[END_REF]. Residual permutations appear to control the family-wise error rate, even for small trees and for large isotropic dimensionalities, and to reach high power (Collyer andAdams 2018, Adams andCollyer 2018). Two points should however be noted: the sample size used here is smaller than the minimum size in these studies; the simulation scheme used by these authors did not acknowledge the peculiarities of isotropy in the tangent space where correlations between coordinates exist due to the Procrustes superimposition [START_REF] Bookstein | The inappropriate symmetries of multivariate statistical analysis in geometric morphometrics[END_REF][START_REF] Klingenberg | Walking on Kendall's Shape Space: Understanding Shape Spaces and Their Coordinate Systems[END_REF]. The family-wise error rate and power given the characteristics of the data (i.e., tree and mean shape) were computed for a continuous covariate, a three-state covariate, and a binary classification. Covariates and shapes evolve according to a Brownian model and to Procrustes constraints on the tangent space. Simulation results are provided as supplementary materials.

Visualization of shape changes

Shape effects were visualized on two meshes, corresponding to the mean shape for Myomorpha, and then for Hystricomorpha, because of the considerable anatomical differences between the two suborders. Meshes were deformed in relation to shape effect, using thin plate spline with the Morpho R package 2.8 [START_REF] Schlager | Chapter 9 -Morpho and Rvcg -Shape analysis in R: R-Packages for geometric morphometrics, shape analysis and surface manipulations[END_REF]. Deformations were either colorized in relation to the signed distance between the predicted and reference meshes or visualized as animated deformations, which are provided as supplementary movies.

Results

Checking tangent space approximation and imbalance effects

Correlation between Euclidean distances in the tangent space and Riemannian distances was very high (r = 0.999), confirming the accuracy of the approximation. The mean shape of Myomorpha and Hystricomorpha was very similar to the mean shape based on weighted averages (a = 2.55°), and angles between specimens were always larger than this angle, suggesting that the shape space was ideally projected for the sample studied, and was not affected by the imbalance between the suborders. It should be recalled here that the lagomorph was aligned as supplementary data only.

Developmental and evolutionary patterns of shape changes

Figure 3. Developmental morphospace and shape changes associated with the phylogenetic mean postnatal trajectory. A. Developmental morphospace where average juveniles of species are represented as the tip of the phylogenetic branches, postnatal trajectories are represented as arrows from the juvenile to adult shapes. Red is for Hystricomorpha species, blue for Myomorpha species, green for Sciuromorpha species and black for Lagomorpha; B: Shape changes associated with the phylogenetic mean postnatal trajectory mapped on the Hystricomorpha mean shape (AP: angular process, CP: condylar process, DI: diastema, CRP: coronoid process, LN: lunar notch, MF: masseter fossa, PCP: postcondylar process, VRMF: ventral ridge of masseter fossa); C: Shape changes associated with the phylogenetic mean postnatal trajectory mapped on the Myomorpha mean shape. Green areas represent expansion from juvenile to adult shapes, whereas beige areas correspond to a compression effect. Gray shapes correspond to the juvenile mean shape (2) or to the modeling of adult mean shape (3) according to the juvenile mean shape plus the phylogenetic mean trajectory " $ !"#$ . See also Movie S1 for dynamic shape changes.

Variation is concentrated on the first two PCs, which account for 59,9% of the total shape variance (Fig. 3A), with 90% for the first eight PCs. Myomorpha and Hystricomorpha are opposed on PC1, which accounts for 42,7% of shape variance. This pattern is expected, since mandible anatomy differs strongly between the two suborders. Developmental vectors (the arrows in Fig. 3A) appear to be mainly orthogonal to this shape divergence and project consistently on PC2. In this 2D projection, these postnatal trajectories in Myomorpha appear to be oriented in similar directions, while the vectors of Hystricomorpha seem less consistent with this general trend. Sciuromorpha are placed close to Myomorpha and their developmental vectors seem to follow the same main direction. The special case of Muscardinus, which has a myomorphous skull structure but belongs to the Sciuromorpha suborder, does not seem to be more similar to the trajectories of Myomorpha species.

The PC2 axis seems to summarize a growth trajectory common among rodents. The phylogenetic mean of postnatal shape trajectory " $ !"#$ , which corresponds to the ancestral pattern of postnatal shape changes, is in close agreement with the shape changes described by PC2 (a = 41.4°, p < 0.0001). Suborder estimates " + of the postnatal trajectory are at 31.1° ± 7.2 of this common growth pattern described along PC2 for Myomorpha, at 58.9° ± 14.9 for Hystricomorpha, and at 49.6° ± 6.6 for Sciuromorpha. These angles are more similar than expected between two random directions (p < 0.0001). On average, 47.3% ± 0.22 of these vectors map onto PC2, with the remainder spread over all the other PCs. The estimation of the postnatal vectors and derived parameters seem robust given the amount of change and species differences to the number of sampled juveniles (Fig. S1). They may also be sensitive to variation in age of juveniles because of the non-linearity of postnatal trajectories [START_REF] Sheets | Studying ontogenetic trajectories using resampling methods and landmark data[END_REF]. Nonetheless, the vectors are still describing a roughly similar direction in the shape space and for instance differences in angles with " $ !"#$ due to age variation are smaller than the observed angles (Table S2).

Postnatal shape changes described by the ancestral pattern of growth " $ !"#$ (Fig. 3BC; Movie S1) imply elongation of the condylar and postcondylar processes, with strong bending of the lunar notch. The angular process expands and subsides while the anterior part of the ventral ridge of the ramus retreats, increasing the ramus. The diastema elongates and flattens. The alveolar region is proportionally larger in the juvenile mandible, whereas in the adult mandible the alveolar region is of similar proportions to the ascending ramus.

The apparent orthogonality observed between the phylogenetic signal and the postnatal development on the first PCs of the developmental morphospace is actually weaker than supposed from the preceding. Whereas the first phylogeny-aligned component is clearly orthogonal to " $ !"#$ or to the suborder estimates " +, additional phylogeny-aligned components are more similar to some of these trajectories, implying that some divergences may have happen along similar directions than the postnatal development (Fig. 4A). The overall angle with the Hystricomorpha estimate " + is larger (Fig. 4A) and it could be noted the striking difference between suborder estimates on the second component (Fig. 4B). These results suggest that the divergences in this suborder are somehow less related to its main pattern of postnatal changes. The overall congruence between evolutionary divergences and developmental directions is even stronger. Krzanowski's common subspace analysis shows that species divergences share some shape changes with the main patterns of postnatal changes. The first five eigenvalues F of H range between 1.67 and 1.96 (for a maximum of 2) and the individual angles G between the eigenvectors of H and the postnatal or evolutionary matrices ranges from 7.7° and 23.9°.

Levels of adult and juvenile disparities

Given the levels of agreement between evolutionary changes and growth directions, an increase in disparity between juveniles and adults could be expected. However, such an increase is observed only in Myomorpha and is only marginal (Fig 4A, Table S3). On contrary in Hystricomorpha, disparity seems to be established in the early stages, as levels of disparity are similar in juveniles and in adults. Disparity in Myomorpha juveniles accumulates more slowly than for the other three clade

´ age groups (Fig. 4B). In particular, disparity based on the main patterns of differentiation between species and developmental stages (the first PCs) is much lower in juvenile than in adult Myomorpha, a pattern that is not observed in Hystricomorpha. This suggests that Myomorpha juveniles have a relatively undifferentiated mandible among species, given the broad pattern of differentiation of the mandible observed at the scale of the order.

Evolution of postnatal shape changes

The observed variation in postnatal shape trajectories (the arrows in Figure 3) appears consistent with a phylogenetic signal (Fig. 5A-C). Species trajectories cluster within suborders both in direction and magnitude. The main pattern of changes shows an opposition between larger postnatal shape changes on negative values of PC1 and smaller ones on positive values (Fig. 5B). These variations agree with a pattern expected under Brownian divergence (Kmult = 0.42, p = 0.04), as was the case for the spreading of species in the developmental morphospace for both adults (Kmult = 0.91, p < 1 ´ 10 -4 ) and juveniles (Kmult = 0.93, p < 1 ´ 10 -4 ). However, the effect size is much smaller for developmental trajectory than for shape (Z = 1.78 versus 4.60 for adults and 4.96 for juveniles), suggesting a less neutral divergence and species-specific postnatal growth trajectories.

Figure 5. Variation of postnatal trajectories. A. Angles between the postnatal trajectories of species (") and the phylogenetic mean (" $ !"#$ ). The lengths of the vectors are proportional to their norm; B. The principal component ordination of the postnatal trajectories (") together with the projection of the phylogeny; C. Magnitude of postnatal shape changes; D. Angles between the estimates of suborder trajectories (" +) and the common growth pattern " $ !"#$ ; E. Angles between the estimates of suborder trajectories (" +); F. Magnitude of the suborder estimates of postnatal changes; G. Angles between " < within suborders Gray shading in panel A, D, E corresponds to the distribution of angles between random vectors with a probability > 0.05. See Figure 3 for color scheme and species acronyms.

In agreement with this phylogenetic signal, the suborder trajectories " + of agree with the common growth pattern (Fig. 5D). The deviation between these suborder estimates (Fig. 5E) is 54.7° ± 17.4° (Myomorpha vs Hystricomorpha), 37.1° ± 9.4° (Myomorpha vs Sciuromorpha), and 45.8° ± 15.1° (Hystricomorpha vs Sciuromorpha). These angles are more similar than expected between two random directions (p < 0.0001). These elements support a relative conservation of the main pattern of postnatal changes of the mandible in rodents. However, the magnitude of postnatal shape changes ‖" +‖ differs between Myomorpha and Hystricomorpha (t12 = 2.70, p = 0.009), with larger vectors in the latter (Fig. 5F).

Once removed the expected patterns under Brownian divergence, the approximately phylogeny-corrected trajectories " = still present important variation within and among suborders (Fig 5G). These angles between phylogeny-corrected vectors " = are large, indicating species specificities and divergence of individual trajectories. For many species comparisons, this angle is still more similar than expected for random vectors. This evidences a relative conservation of the ancestral main pattern of growth.

To assess the influence of life history traits and tooth growth on this remaining variation in trajectories, a phylogenetic generalized least-squares was performed. Gestation time and duration of pre-weaning period were pooled within suborders prior to analysis as they present a very high correlation with deep node in the phylogeny (Fig. 1). It appears that only type of tooth growth has an effect after controlling for Brownian expectation of evolutionary divergence (Table 1). Simulations show that appropriate type-I error is returned and adequate power could be obtained for large effects (Table S4 and Fig. S2). Pairwise comparisons between types of tooth growth show that hypselodonty has a different effect from brachyodonty (d = 0.132 Z = 3.36, p = 0.002) and hypsodonty (d = 0.117, Z = 1.83, p = 0.05), whereas the effect of brachyodonty is similar to that of hypsodonty (d = 0.087, Z = 0.934, p = 0.17). Figure 6. Shape changes associated with dental growth: brachyodont versus hypselondont. Effects are mapped on the mean shape of myomorph specimens. Green areas represent expansion from brachyodont to hypselodont shapes, whereas beige areas correspond to a compression effect. See also Movie S2 for dynamic shape changes. (AP: angular process, CP: condylar process, DI: diastema, CRP: coronoid process, LN: lunar notch, MF: masseter fossa, PCP: postcondylar process, VRMF: ventral ridge of masseter fossa). 1: labial view; 2: upper view; 3: lingual view. Gray shapes correspond to the least square estimates of the B: brachyodont or C: hyselodont shapes. 1: labial view; 2: lingual view.

Shape changes in relation to hypselodonty (Fig. 6 and Movie S2) show an extension of the condylar process, which rises vertically and becomes thinner, and a thicker diastema, especially in its upper part. The molar area increases in volume, as does the anterior part of the ventral ridge of the masseter fossa. The alveolar region is more robust and extends further on the lingual side.

Discussion

Conservation of the main patterns of growth and evolutionary divergences

Our main results highlight a general pattern of shape changes during postnatal development across Rodentia, with similar local bone growth. The mandible shape in juveniles is generally elongated and flattened, with relatively short processes. These processes later expand, increasing their relative height and thus enlarging the masseter fossa. This ontogenetic modification of the posterior part of the mandible, well described in laboratory mice [START_REF] Swiderski | The complex ontogenetic trajectory of mandibular shape in a laboratory mouse[END_REF], and observed in Sciuridae (Zeldtich et al. 2016), is here shown to be conserved across the main clades of Rodentia. This pattern is undoubtedly related to the development of the muscular complex inserting on the mandible. As these muscles gain volume during growth, their contact surface with the bone must expand. Similarly, the different processes become more robust after weaning to support the attachment tension of the muscles. This musculoskeletal system articulates with the skull (cranium and mandible) to allow mastication, which will impose greater mechanical constraints at this life stage [START_REF] Scott | Teaching an old jaw new tricks: diet-induced plasticity in a model organism from weaning to adulthood[END_REF]. This general pattern of growth appears orthogonal to the main shape divergence between suborders, which opposes the mandible shape of Hystricomorpha with the ones of Myomorpha and Sciuromorpha. This contrast corresponds to gross differences in anatomy, the mandibular structure of [START_REF] Tullberg | Ueber das System der Nagethiere: eine phylogenetische Studie[END_REF], opposing hystricognaths and sciurognaths based on the opposite orientation of the angular and condyloid processes.

Evolutionary divergences between species appear nonetheless much more related to the ontogenesis and its variation than expected from the orthogonality of their main directions, suggesting that developmental processes such as heterochrony could have played a role in species divergences. Based on the sampled species, Hystricomorpha exhibit developmental patterns that are much more genus-specific than those of Myomorpha species, where the trajectories are more similar to the common developmental pattern. This difference seems to be mainly related to the difference in species divergence between suborders. Myomorpha have a common juvenile shape that diverges postnatally among sub-families. In Hystricomorpha, beside the common juvenile pattern, the mandible shape of juveniles presents some species-specific features from the early stages, and the magnitude of postnatal shape changes is smaller. Hystricomorpha being not precisely age controlled and in smaller numbers compared to Myomorpha, variation observed in Hystricomorpha could be more biased by some age or sampling effects (Fig S1 andTable S2). These biases may partly explain some variation observed between the suborders. Nonetheless, the differences between suborders could also be explained by the vast differences between the ecological traits for each suborder (Wolff and Sherman 2009), including reproductive strategy, gestation, and weaning (acting directly on pup development), diet (leading to more or less mechanical stress on the masticatory apparatus), and taxonomic richness (reflecting the diversity of a clade). Myomorpha newborns share common anatomical features, whereas the longer gestational period in Hystricomorpha allows longer in-utero development, producing a skeletonized and differentiated organ [START_REF] Huggett | The relationship between mammalian foetal weight and conception age[END_REF]. In this respect, the peculiarity of gerbils (Meriones) is revealing as this species has a longer gestation than other Myomorpha species and presents at seven days old a more differentiated juvenile mandible as the Hystricomorpha species.

This developmental difference between Myomorpha and Hystricomorpha suggests that juvenile disparity at birth is partly driven by timing differences in the skeletonization and differentiation processes in relation to variation in gestational duration. Therefore, the potential for plasticity after birth is likely to be different between the two suborders since the mandible at birth is not at the same state of ossification and differentiation.

Evolution of postnatal trajectories

Variation in gestational length may explain changes in postnatal trajectories because it modifies the onset of the postnatal period. It could thus explain variation in magnitude of postnatal changes as juveniles could be more or less differentiated and skeletonized at birth. The observed differences in magnitude between suborders could be related to such factor. As developmental trajectories are strongly non-linear [START_REF] Green | Developmental nonlinearity drives phenotypic robustness[END_REF], variation in the onset of postnatal growth will affect the linear approximation of the postnatal trajectory leading to variation in angles in relation to differential duration of the gestation. Again, variation between suborders seems to be related to such an effect. However, once the expected effect of Brownian divergence was factored in, no effect of the gestational duration was observed to explain the residual differences in postnatal shape changes. The separation between suborders is ancient and occurred during the Paleocene [START_REF] Swanson | A phylogenomic rodent tree reveals the repeated evolution of masseter architectures[END_REF]. Deep-rooted phylogenic signal strongly structures traits and aggregates them, as gestational length or weaning timing, within clades. These correlations potentially mask differences related to ecological traits.

The weaning period and dietary diversification occurring within this period could also play a role in the divergences of postnatal trajectories as different behaviors and feeding habits lead rodents to adapt their chewing movement. Rodents can gnaw with their incisors as well as chew with their cheek-teeth [START_REF] Byrd | Mandibular movement and muscle activity during mastication in the guinea pig (Cavia porcellus)[END_REF][START_REF] Cox | Functional Evolution of the Feeding System in Rodents[END_REF][START_REF] Hiiemäe | A cinefluorographic study of mandibular movement during feeding in the rat (Rattus norvegicus)[END_REF]). These movements, acquired more or less rapidly during weaning, could induce new mechanical constraints and thus change the pattern of bone remodeling [START_REF] Jacobs | Rodentia: Extraordinary diversification of a morphologically distinctive and stereotyped order[END_REF]. Changes in the consistency of food are known to induce plasticity of the mandible, with a large effect on the angular and coronoid processes, modifying the correlation between the alveolar and muscle-bearing regions and their biomechanics [START_REF] Anderson | Adaptive plasticity in the mouse mandible[END_REF]. In Myomorpha species, weaning is rapid and therefore mechanically brutal [START_REF] Curley | The Meaning of weaning: influence of the weaning period on behavioral development in mice[END_REF]. This rapid upheaval in diet could explain the major changes observed in morphology. The feeding system must change from the sucking movement to active chewing in only a few days. In mice, bite force increases drastically in the days preceding weaning [START_REF] Ginot | Decoupled ontogeny of in vivo bite force and mandible morphology reveals effects of weaning and sexual maturation in mice[END_REF]). In the Hystricomorpha species, parental care is different, with long behavioral weaning and less brutal dietary weaning, over a much longer period [START_REF] Wolff | Rodent societies: an ecological and evolutionary perspective[END_REF]. As the Hystricomorpha species, at least the ones sampled in this study, are born with a more differentiated shape and present postnatal shape changes of smaller magnitude during a longer period, dietary changes at weaning have likely less impact on the patterns of correlated bone remodeling. Guinea pigs (Cavia) have a short weaning timing of 18 days comparable to the one of the Myomorpha species but present nonetheless a magnitude and direction of postnatal trajectory similar to other Hystricomorpha species. As with gestation, no effect of the weaning period was observed to explain the residual differences in postnatal shape changes once the expected effect of Brownian divergence was factored in, but again weaning period timing is strongly aggregated within each clade.

Only tooth growth, especially molar hypselodonty, appears to produce a special signature on the postnatal shape changes of the mandible. Within Myomorpha, hypselodonty in arvicolines may explain the divergence of their postnatal trajectory, as their unrooted prismatic teeth strongly modify the dental alveolar region. A comparable effect of tooth growth has been shown on individual postnatal ontogenies in mice [START_REF] Swiderski | The complex ontogenetic trajectory of mandibular shape in a laboratory mouse[END_REF]. During development, bone and teeth interact with each other and these interactions may explain important changes either in bone or in tooth. For instance, tooth-bone mechanical interaction explains the lateral offset of molar cusps during development, a repeated evolutionary innovation in the mammalian tooth pattern [START_REF] Renvoisé | Mechanical constraint from growing jaw facilitates mammalian dental diversity[END_REF]. This interaction influences also the mandible shape as it is modified with arrested tooth development in mutant mice [START_REF] Paradis | Jaw growth in the absence of teeth: the developmental morphology of edentulous mandibles using the p63 mouse mutant[END_REF][START_REF] Boughner | P63 expression plays a role in developmental rate, embryo size, and local morphogenesis[END_REF]. Thus, the bone must adapt to the dynamic process of tooth growth, which imposes new biomechanical stresses. This dynamic reconfiguration of bone strain imposes a new spatialization of bone formation and remodeling [START_REF] Martínez-Vargas | Postnatal mandible growth in wild and laboratory mice: Differences revealed from bone remodeling patterns and geometric morphometrics[END_REF], that may have evolutionary outcomes in both tooth and bone shapes.

Our results suggest that evolutionary changes in the postnatal development of the mandible could be related to muscle and tooth development. The mandible responds dynamically to both the growth of the dental alveoli and the development of the masticatory muscles [START_REF] Swiderski | The complex ontogenetic trajectory of mandibular shape in a laboratory mouse[END_REF]. The biomechanics of mastication functionally integrates the alveolar and ramus regions from their direct epigenetic interactions [START_REF] Zelditch | Modularity of the rodent mandible: Integrating bones, muscles, and teeth[END_REF]. Changes in environment (e.g., diet or behavior) could induce changes in the spatialization of biomechanical strain and therefore of bone remodeling, because of the intricate relationship between the mandible and tooth and muscle growth. If tooth or muscle growth drives mandible plasticity, natural selection on these anatomical elements will induce correlated changes in the mandible via epigenetic interactions, which may in turn be subject to genetic accommodation [START_REF] West-Eberhard | Developmental Plasticity and Evolution[END_REF][START_REF] Uller | Developmental plasticity and evolutionary explanations[END_REF]. For example, at a smaller scale, the repeated dietary adaptation in mice observed following their invasion of sub-Antarctic islands seems to agree with a similar scenario, with mandible plasticity as the main process at the early stages of the invasion [START_REF] Renaud | Divergent in shape and convergent in function: Adaptive evolution of the mandible in Sub-Antarctic mice[END_REF].

Comparison with the cranium

Despite a common pattern of postnatal changes in relation to muscle and tooth growth, mandibles are quite different between suborders suggesting that ontogeny did not constrained anatomical diversity. At its highest level (hystricognaths versus sciurognaths), this diversity occurs mainly in an orthogonal direction of the common postnatal growth pattern. This independence of anatomical diversity to developmental constraints was similarly observed for the cranium [START_REF] Wilson | Allometric disparity in rodent evolution[END_REF]Sanchez-Villaga 2010). However, at a lower scale, the diversity of postnatal trajectories (i.e. variation in the magnitude of changes and angles) appears to be specific to the suborders. This observation contrasts with the constrained variation of allometric patterns observed on the cranium [START_REF] Wilson | Allometric disparity in rodent evolution[END_REF]Sanchez-Villaga 2010). Within the head, the cranium is a composite unit composed of multiple bones, and this assemblage is constrained by the development of teeth [START_REF] Renaud | Developmental constraints revealed by co-variation within and among molar rows in two murine rodents[END_REF]), facial muscles [START_REF] Ravosa | Using "Mighty Mouse" to understand masticatory plasticity: myostatin-deficient mice and musculoskeletal function[END_REF]), and internal organs (i.e. the brain, Richtsmeier and Flaherty 2013, or the olfactory and auditory organs, [START_REF] Barone | Anatomie comparée des mammifères domestiques[END_REF]). This complexity may explain the low skull disparity in relation to lineage diversification described in Rodentia [START_REF] Alhajeri | Disparity and evolutionary rate do not explain diversity patterns in Muroid rodents (Rodentia: Muroidea)[END_REF], or in the adaptive radiation of some families [START_REF] Maestri | The ecology of a continental evolutionary radiation: Is the radiation of sigmodontine rodents adaptive?[END_REF], or the relatively overlapping growth pattern of the skull observed among rodent clades [START_REF] Wilson | Allometric disparity in rodent evolution[END_REF]. The higher complexity of skull development is thought to reduce its dependence on muscle and tooth development [START_REF] Swiderski | The complex ontogenetic trajectory of mandibular shape in a laboratory mouse[END_REF], and this difference in complexity probably explains the contrast observed between these two skeleton elements of the head.

Conclusion

In conclusion, the main evolutionary and developmental patterns appear to be along orthogonal directions of the shape space. Beside these main patterns, most of the shape divergences arise along postnatal directions of growth. The filling of the morphospace varies between rodent suborders being mainly during the gestational period in Hystricomorpha and during the postnatal growth in Myomorpha. The postnatal trajectories of Myomorpha subfamilies are closer, but present a much stronger magnitude of changes over a shorter growth period. Among Hystricomorpha, part of the observed adult shape is set up prenatally, and most postnatal trajectories are genus-specific, which agrees with non-linear developmental trajectories over longer gestational periods. Juvenile mandibles present similarities among suborders, with a flattened shape. Mandible shape then diversifies during growth in relation to muscle and tooth development, with epigenetic interactions coordinating the changes in the alveolar and muscle-bearing regions. This functional integration by biomechanical interactions might have favored evolutionary changes driven by developmental plasticity, but the importance of this process is likely to differ between rodent suborders.

Figure 1 .

 1 Figure 1. Summary of data on phylogeny, tooth growth (Ungar 2010) and reproductive traits (gestational duration, weaning time and offspring by litter) for the species sampled for the study (Wilson et al. 2016, De Magalhaes and Costa 2009). Color code for outgroup and suborders: a) Lagomorpha in black; b) Myomorpha in blue; c) Sciuromorpha in green; d) Hystricomorpha in red.

Figure 2 .

 2 Figure 2. Landmarks and semilandmarks on 3D mandibular surface. Black dots are manual 3D landmarks, and small yellow dots represent the semilandmark template. (AP: angular process, CP: condylar process, DI: diastema, CRP: coronoid process, LN: lunar notch, MF: masseter fossa, PCP: postcondylar process, VRMF: ventral ridge of masseter fossa). A: labial view; B: occlusal view; C: lingual view.

Figure 4 .

 4 Figure 4. Developmental changes in disparity. A. Comparison of levels of disparity between juveniles and adults within suborders. B. Accumulation curves of juvenile and adult disparities along the principal components of the developmental morphospace. Blue is for Myomorpha and red is for Hystricomorpha.

Table 1 .

 1 Phylogenetic generalized least-squares of postnatal trajectories.

	Source	Df SS	MS	r 2	F	Z	Pr(>F)
	Gestation length 1 0.02228 0.022283 0.05460 1.0948 0.37464 0.3675
	Weaning age	1 0.03433 0.034328 0.08411 1.6866 1.10556 0.1411
	Litter per year	1 0.02846 0.028462 0.06974 1.3984 0.79606 0.2203
	Tooth growth	2 0.08831 0.044156 0.21639 2.1695 1.94286 0.0252
	Residuals	10 0.20353 0.20353 0.49871

Acknowledgements

We are grateful to all the curators who allowed us to study their collections: Violaine Nicolas (National Natural History Museum Paris), Loïc Costeur (Natural History Museum Basel), Olivier Pauwels (Royal Belgian Institute of Natural Sciences), Laurent Vallotton and Manuel Ruedi (Natural History Museum of Geneva). Mathilde Tissier (Hubert Curien Pluridisciplinary Institute) is also thanked for providing access to Cricetus cricetus specimens. We thank Lauriane Poloni (Biogéosciences) for help with CT scanning, and Remi Laffont (Biogéosciences) for help in data curation. We are grateful to Philip Cox, Vincent Debat, Lionel Hautier and Miriam Leah Zelditch for helpful comments on early versions of the manuscript. We are also grateful to Carmela Chateau-Smith for English editing. The Gismo platform and morphOptics are acknowledged for providing access to 3D scanning. This research project was supported by the AP EPHE-2019.

Data availability

The landmark data that supports the findings of this study are available on Dryad.

Author contributions

All authors contributed to the study design. MD collected species in museum and conducted the phenotyping. MD and NN analyzed the data. All authors contributed with edits and comments to analyses and wrote the manuscript.

Conflict of interest

Authors declare no conflicts of interest.

Ethical approval

All procedures on laboratory animals were approved by the local ethical committee of the University of Burgundy in Dijon (Project APAFIS#18405-2019011014262528). Error rate and power of the phylogenetic generalized least-squares using Procrustes sum of squares and residual permutations Shapes were simulated in the tangent space following its linear approximation [START_REF] Bookstein | The inappropriate symmetries of multivariate statistical analysis in geometric morphometrics[END_REF][START_REF] Klingenberg | Walking on Kendall's Shape Space: Understanding Shape Spaces and Their Coordinate Systems[END_REF] given the mean shape of our sample and an isotropic variation. This simulation scheme respects the constraints that the Procrustes superimposition imposes on the variation. This is the main difference to the simulations of [START_REF] Adams | Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations[END_REF] where multiple traits were simulated with a similar effect on each dimension. Also, our focal tree is half smaller than the tree used in this previous study. The norm of the simulated shape change was scaled to a certain proportion of the total Procrustes variance (at most 50%). 100 simulations of shape effect and residuals were done and for each of them 100 Brownian evolution of the covariate were simulated. Two kind of covariate were simulated: a continuous one and a three-state character. Analysis with a random binary classification was also realized. Probability of the effect was estimated with phylogenetic GLS within the package geomorph 3.2.1 [START_REF] Collyer | Phylogenetically aligned component analysis[END_REF]) and evaluated using 1000 residual permutations with the package RRPP 0.5.2 [START_REF] Collyer | RRPP: An r package for fitting linear models to high-dimensional data using residual randomization[END_REF]. The family-wise error rate was computed on as the proportion of significant pGLS when no covariate effect was simulated and again this error rate was averaged across the 100 shape configurations (Table S4). For each of the three models, the expected FWER is lower than 5%. S4. Error rate for the pGLS given the 16 species tree Effect fwer 2.5 th quantile 97.5 th quantile Quantitative covariate 0.042 0.005 0.08 3-state covariate 0.049 0.01 0.11 Binary class 0.039 0.01 0.08 Power was computed as the proportion of the simulations where the probability was lower than 0.05 and the average across the 100 shape configurations was reported (Figure S2).

Supplementary Table

Figure S2.

Power of pGLS given the 16 species tree and the three kind of covariate or grouping. X-axis on the left panel reported the ratio of the norm of the effect b (effect to be simulated) to the total variance (the sum of eigenvalues), and the observed r 2 on the right. Shade area represent 95 th interval of the 100 shape simulations and lines their averages. Power reaches 80% for an effect size equal to 0.35 (quantitative covariate), 0.37 (3-state covariate), and 0.39 (binary classification). Movies S1. Shape changes according to the ancestral estimate of postnatal growth trajectory apply to the estimate of the juvenile mean shape of either the Myomorpha or the Hystricomorpha.

Movies S2. Shape changes according to the effect of the hypselodonty estimated from the pGLS and applied to the Myomorpha mean shape (based on a brachyodont shape).