Masticatory system integration in a commensal canid: interrelationships between bones, muscles and bite force in the red fox

Skull morphology has been demonstrated to be the complex product of phylogeny, development, mechanical processes, compromises produced by competing demands, and epigenetic constraints (Bels and Herrel, 2019; Smith, 1993). The head is involved in many fundamental functions in vertebrates (protection of the sensory organs and brain, and functioning of the digestive and respiratory tracts; Santagati and Rijli, 2003), and has thus been the subject of many studies. In particular, the relationships between the mechanical components of the jaw have been a subject of major interest as they contribute to the understanding of the evolutionary processes that have driven variation in the jaw system (Cornette et al., 2015a). The mandibles rotate up or down relative to the cranium, the movements being driven by contraction of the jaw adductors, thus generating the bite force, which is an excellent indicator of the performance of the jaw system (Anderson et al., 2008; Binder and Valkenburgh, 2000; Dessem and Druzinsky, 1992; Nogueira et al., 2009). Numerous studies have documented the biomechanics of the jaws in a variety of organisms and attempted to link this to prey capture mode, dietary specialisation, competition or non-feeding and environmental variables (e.g. Bels, 2006; Bels and Herrel, 2019; Bels et al., 2012; Cornette et al., 2015b; Fabre et al., 2018; Gueldre and Vree, 1990; Hannam and Wood, 1989; Hartstone-Rose et al., 2012; Herrel and Aerts, 2004; Herrel et al., 1998b, 2008; Herring et al., 2001; Nogueira et al., 2009; Perry et al., 2011; Slater and Van Valkenburgh, 2009; Tseng and Flynn, 2015a,b,c, 2018).

Previous studies focusing on canids have documented the functional relationships of the cranium with the adductor muscles (Penrose et al., 2016, 2020) and the mandible (Curth et al., 2017; Curth, 2018). The relationships between these structures and bite force have often been explored using bite force estimations based on skull measurements (i.e. dry skull method: Ellis et al., 2008, 2009; Forbes-Harper et al., 2017; Thomason, 1991). Yet, this approach does not take into account the macroscopic arrangement of muscle fibres (i.e. muscle architecture: muscle volume, fibre length, fibre type and pennation angle). A good overall measure of this architecture and a representative estimate of muscle force is the physiological cross-sectional area (PCSA; Haxton, 1944; Martin et al., 2020). Moreover, differences in bone shape and lever arms may change the relative arrangement of muscles on the skull with respect to the teeth or the temporomandibular joint, which will also influence bite force (Taylor and Vinyard, 2013). As such, the inclusion of data on muscle architecture and the position of the muscles on the skull may influence estimates of the magnitude of the bite force and jaw-closing speed (Penrose et al., 2020). Muscle dissections hence enable better estimation of bite force. Previous studies have further shown good correspondence between in vivo data and results from biomechanical models based on dissection data (Herrel et al., 1999, 2008; Meyers et al., 2018). Other studies have used finite element analyses (Bourke et al., 2008; Penrose et al., 2020; Wroe et al., 2007) to explore possible adaptations to predation (Radinsky, 1981; Slater et al., 2009; Van Valkenburgh and Koepfli, 1993). Unfortunately, direct in vivo measurements of bite force in canids are scarce. Ellis et al. (2008) recorded bite force data under anaesthesia for 20 dogs of various breeds, Lindner et al. (1995) recorded in vivo data on 22 dogs of various breeds in vivo and Brassard et al. (2020a) recorded in vivo data on four Malanese dogs trained for attack. To date there are no in vivo data for bite force in other canids such as red foxes.

Foxes (Vulpes vulpes) are the most widespread wild canids in the world (Schipper et al., 2008). The success of the red fox in a wide variety of ecological contexts suggests that they present functional adaptations to diverse environments and resources. Foxes are moderately carnivorous opportunists that feed mostly on prey of small or average size (lagomorphs, birds, murid rodents and voles, but also larger prey such as sheep and wild pig) and fruits (Cavallini and Volpi, 1996; Dell'Arte and Leonardi, 2005, 2009; Doncaster et al., 1990; Rosalino and Santos-Reis, 2009). They have long and narrow jaws that can close quickly by a lengthening of the jaw out-lever arm. They thus have an advantage for capturing small and fast moving prey, albeit at the expense of a decrease in bite force at the tip of the jaw (Christiansen and Wroe, 2007; Herring and Herring, 1974; Perry et al., 2011; Radinsky, 1981; Santana, 2016; Slater et al., 2009; Van Valkenburgh and Koepfli, 1993), unless the size and orientation of the jaw musculature compensate (Jaslow, 1987). Yet, the morphological and functional variability of this species has only been briefly described. Forbes-Harper et al. (2017) estimated bite force using the dry-skull method for a population of over 300 Australian red foxes. However, no previous studies have explored the variation in muscle architecture in the red fox, and how the 3D shape variation in the cranium or mandible is related to jaw muscle morphology and consequently bite force. Further, no in vivo bite force measurements are available in the literature, yet these are essential to validate any biomechanical model used to calculate bite force.

Unlike in dogs (Canis lupus familiaris) and despite their commensal nature across a large part of the range, phenotypic variation in red foxes is mostly driven by natural selection. Consequently, they are an excellent model to compare with domestic dogs where phenotypic differences are almost exclusively the result of artificial (intentional) selection. As the jaw system plays a major role in feeding and predation, it is likely to be a highly integrated system that is subject to natural selection pressures. Wild species are only rarely compared with domestic dogs although many studies have attested to the consequences of domestication on the morphological variability of the cranium and mandible (Curth, 2018; Curth et al., 2017; Drake and Klingenberg, 2008, 2010; Machado et al., 2018; Sánchez-Villagra et al., 2016; Selba et al., 2019). Consequently, our understanding of how domestication may impact the functional properties of the jaw system remains limited.

Here, our aim was to explore the variability in and the interrelationships between the shape of the cranium, the shape of the mandible, muscle architecture and bite force in red foxes, to assess the functional impact of the variability in shape and muscle architecture. To do so, we took the following steps. (1) We explored the variability in the shape of the cranium (skull without the mandible) and mandible of 68 foxes Vulpes vulpes from France using 3D geometric morphometrics. (2) We quantified the variability in the jaw muscle architecture (pennation angle, mass, PCSA) of 65 of these animals by means of dissection. (3) We explored the variability in estimated bite force, using a 3D static biomechanical model based on the origin and insertion of the adductor muscles and the PCSAs obtained from dissection, and for which we compared model output with in vivo measurements. Then, we characterised the interrelationships between the components of the masticatory system by testing the (4) covariation between cranial and mandible shape, (5) correlation and covariation between muscle and shape data, (6) and also correlation and covariation between estimated bite force and muscle and shape data. Finally (7), we compared the results obtained here for foxes with results obtained previously for dogs (Brassard et al., 2020a,b,c) to compare the variability in cranial and mandible shape and jaw muscle architecture as well as the patterns of integration between the domestic dog and the commensal red fox. We predicted that the variability in cranial and mandibular shape would be lower in the red fox than in domestic dogs based on previous studies (Drake and Klingenberg, 2010). We also expected stronger correlation and covariation between bone shape and muscles or estimated bite force in the red fox as its jaw system is principally under the influence of natural selection. This study should thus contribute to a better understanding of the evolutionary processes that drive jaw biomechanics in wild canids and the modifications induced by artificial selection in the domestic dog.

The dataset was composed of the heads of 68 fresh-frozen Vulpes vulpes (Linnaeus 1758), including 64 red foxes from the wild and four silver foxes from a wildlife virology testing centre. All ethical guidelines and recommendations on animal welfare were followed. The silver fox is a melanistic form of the red fox that has been selected over decades by humans for its fur, and belongs to the same genus and species Vulpes vulpes (Dugatkin, 2018; Trut, 1999). Detailed information on the sample is available in Table S1. Sixty-five of these heads were dissected. Three heads were not dissected but directly prepared for shape analyses because they were preserved in formaldehyde, which may impact muscle architecture data. Fifty-eight crania and sixty-eight mandibles were well preserved enough to be used for shape analyses, after cleaning and drying. Foxes were classified into several age groups depending on the aspect of the cranial sutures. Four foxes were juveniles with deciduous teeth, a very porous mandible and unclosed cranial sutures. Young foxes represent foxes for which the basispheno-basioccipital suture is still open (<8–10 months for dogs according to Barone, 2010) and the mandible is still porous. Old foxes represent foxes with a closed interfrontal suture and worn dentures (>3–4 years). The other foxes were intermediate adults.

The cranium, mandible and jaw adductor muscles act jointly to produce jaw motion and estimated bite force. We used a similar biomechanical model to the one described by Herrel et al. (1998a,b) that takes into account the 3D coordinates of origin and insertion and the PCSA of the jaw muscles to calculate the moments exerted by each muscle and to deduce the estimated bite force, joint force and angle of the joint force (Fig. 1C).

In this model, the cranium is positioned in a reference frame whose centre is located at the right temporomandibular joint, whose x-axis runs through the long axis of the mandible to the first incisor tooth, and whose y-axis is directed towards the top of the cranium and perpendicular to x. The z-axis runs from the midline outwards perpendicular to the other two axes. In this reference frame, the 3D coordinates of origin and insertion of the adductor muscles were recorded using a microscribe. We approximated the centroid of the origin and insertion areas of the muscles based on observations from our dissections. The 3D coordinates of three bite points (point of application of estimated bite force, which is the opposite of the food reaction force) were also recorded. The first point is at the first incisor tooth (BPi in Fig. 1), the second one is just behind the lower canine tooth (BPc in Fig. 1) and the last one is located on the caudal part of the lower carnassial tooth, which corresponds to the contact area between P⁴ and M₁ (BPm in Fig. 1). These locations were chosen because they are essential during feeding in canids.

In this 3D lever model, the lower jaw rotates around the condylar process of the mandible (the centre of the system) following a gape angle of 0–40 deg. We do not take into account translational movements as they are negligible. All bite points then rotate in an arc for which the radius corresponds to the shortest distance from the condylar process of the mandible to the point of application of the estimated bite force. At static force equilibrium, the sum of the moments of the external forces (force in the joint, force at the bite point and force exerted by each muscle) is zero. In other words, the sum of the vectorial products of the in-lever moment arms and the adductor muscle forces (for both sides, which are considered symmetric) is equal to the vectorial product of the out-lever moment arm and the estimated bite force. The magnitude of each moment corresponds to the numeric product of force magnitude and the shortest distance between the centre of the system and the line of action of the force (i.e. the effective lever arm or moment arm). The magnitude of the muscular forces was established by multiplying the reduced PCSA by a conservative muscle stress estimate of 30 N cm⁻² (Herzog, 1994). The effective lever arms were calculated from the recorded coordinates.

Given that we do not know the direction of the estimated bite force (opposite to the food reaction force, which may depend upon the shape, texture and position of the food item as well as the shape and position of the teeth; Cleuren et al., 1995; Aerts et al., 1997), we calculated the effective out-lever arm for a large range of angles thereof (set to vary between −40 and −140 deg with respect to the lower jaw; indicated as ‘AFRF’ in Fig. 1C). We calculated the estimated bite forces for several mouth opening angles (0, 20 and 40 deg; indicated as ‘gape angle’ in Fig. 1C). The magnitude and orientation of the forces in the joint (indicated as ‘JF’ in Fig. 1C) were estimated as well because, at static equilibrium, the sum of the external forces (muscle and estimated bite forces) is zero.

The input for the model, therefore, consists of the PCSA of the jaw muscles, muscle origins and insertions, mouth opening angle, and the point of application of the estimated bite force. Model output consists of the magnitude of the estimated bite forces, the magnitude of the joint forces, and the orientation of the joint forces at any given orientation of the food reaction forces. An R script for the calculation of the estimated bite force is available on request from the corresponding author. Only the estimated bite forces of the foxes with a well enough preserved cranium were estimated (60 individuals).

In vivo bite force data were recorded on 10 awake restrained silver foxes (five males and five females) at the Laboratory for Rabies and Wildlife in Nancy (France). We used a piezoelectric isometric Kistler force transducer (range ±500 N; 9203, Kistler Inc., Winterthur, Switzerland) linked to a charge amplifier (type 5995A, Kistler Inc.), similar to the set-up presented in Herrel et al. (1999) and Aguirre et al. (2002). The distance between the two steel bite plates was adjusted so that the foxes bit at a gape angle of about 20 deg when biting at the front of the jaw and 30 deg when biting at the molars. The tips of the bite plates were covered with a thin medical cloth tape (which was changed between each animal) to avoid direct contact of the teeth with the metal. The foxes were placed on a table and manually restrained and the transducer placed either at the level of the incisor teeth or behind the major cusps of the carnassial teeth, therefore making contact with upper premolar tooth P⁴ and molar tooth M₁, and lower molar teeth M₁ and M₂. We performed five consecutive trials for each animal and retained the maximal bite force recorded across the trials for analyses. One-sided Welch's tests were performed to compare the mean of the in vivo bite forces with the mean of the bite forces estimated using the biomechanical model for a gape angle of 20 and 30 deg, for the two bite points (on the incisor and molar teeth).

All statistical analyses were run in ‘R’ version 3.6.0 (2019-04-26). The patterns of morphological variation and covariation with muscle data or estimated bite force were explored using geometric morphometric analyses. 3D models of all the mandibles and one cranium were obtained from photogrammetry using the ‘Agisoft PhotoScan’ software (^©2014 Agisoft LLC, St Petersburg, Russia). Twenty-five landmarks, 190 sliding semi-landmarks on curves and 185 sliding semi-landmarks on surfaces were placed on the mandible of each specimen (Fig. 2; Table S1) using the software ‘Landmark’ version 3.0.0.6 (^©IDAV 2002-2005; Wiley et al., 2005). The landmarks were slid and transformed into spatially homologous landmarks using a sliding semi-landmark procedure implemented in the ‘Morpho’ package (version 2.7) in R (Bookstein, 1997; Gunz et al., 2005; Schlager, 2013). Fifty-four landmarks were placed on one side of the cranium, using a microscribe (Fig. 2; Table S1). A mirror was then applied to obtain the symmetrical landmarks compared with the sagittal plane, using the function ‘mirrorfill’ from the package ‘paleomorph’, and leading to a total of 108 landmarks.

Generalised Procrustes analyses (GPA; Rohlf and Slice, 1990) were performed using the function ‘procSym’ (Dryden and Mardia, 2016; Gunz et al., 2005; Klingenberg et al., 2002) from the package ‘Morpho’. Allometries in bone shape and muscle morphology (PCSA and mass) were explored using the function ‘procD.lm’ (Adams and Collyer, 2016; Adams and Collyer, 2017; Anderson, 2001; Anderson and Braak, 2003; Collyer et al., 2015; Goodall, 1991) from the package ‘geomorph’. Allometry-free coordinates were calculated using the function ‘CAC’ (Mitteroecker et al., 2004) and ‘showPC’. Residual muscle data and residual estimated bite forces were obtained from the regression of the log₁₀-transformed muscle or estimated bite force data on the log₁₀-transformed centroid size of the mandible [or cranium whenever appropriate for the subsequent two-block partial least-squares analyses (2B-PLS) analyses], using the function ‘lm’.

Principal component analyses (PCA) were performed using the function ‘procSym’ from the package ‘Morpho’ based on the mandibular or cranial coordinates of all aligned specimens, on allometry-free coordinates, on the PCSA of all muscles, on the residual PCSA of all muscles, on muscle mass and on residual muscle mass. The deformation of the mandible or cranium to the consensus of the GPA was used as a reference for all further visualisations.

To investigate the drivers of cranial and mandibular shape variation, we performed Procrustes ANOVA with permutation procedures on the coordinates from the GPA using the function ‘procD.lm’ from the package ‘geomorph’. We considered the mass of the three main adductor muscle groups, their fibre lengths, pennation angles, PCSA and the centroid size (of the cranium or mandible) as explanatory variables, to increase statistical power. For each muscle complex, we considered the sum of the masses or PCSAs and the mean of the fibre lengths or pennation angles of the constituent bellies. Data were log₁₀-transformed. We performed Procrustes ANOVA using the function ‘procD.lm’ from the package ‘geomorph’. We performed several multiple or simple regressions with log₁₀-transformed centroid size, age, sex and muscle data as explanatory variables. For these analyses, we considered the three main muscle complexes (masseter, temporalis and pterygoid). We used the ‘shape.predictor’ function and Avizo 8.1.1. software to visualise the effect of the variation in the PCSA of the temporal, masseter and pterygoid muscles on the shape of the cranium and mandible.

To explore the patterns of covariation, we used 2B-PLS analyses with the function ‘pls2B’ (Rohlf and Corti, 2000). P-values were calculated based on 1000 permutations. We tested the covariation between shapes or allometry-free shapes and raw or residual muscle masses or PCSAs. Z-scores were finally calculated to compare the PLS coefficients with the function ‘compare.pls’ from the package ‘geomorph’.

The proportions of the masseter, temporal and pterygoid muscles to the total mass of the adductor muscles were calculated as the sum of the masses of all the bundles belonging to a functional group.

To identify the relative contribution of the moment exerted by each of the adductor muscles on the moment of the estimated bite force, we calculated the mean and standard deviation of this ratio for each bundle and four gape angles (0, 20, 30 and 40 deg).

To identify the main drivers of estimated bite force variation, we performed multiple linear regressions using the function ‘lm’ with the masses, fibre lengths, pennation angles and PCSAs of the main muscles, and the centroid size of the mandible as explanatory variables. We considered the estimated bite force for a food reaction force orientation of 90 deg and a gape angle of 20 deg. The data were log₁₀-transformed. For this analysis, we considered the three main muscular complexes (masseter, temporal, pterygoid) to increase statistical power. For each complex, we considered the sum of the masses or PCSAs of the constituent bellies and the mean of the fibre lengths or pennation angles. The best-fitted model was obtained from stepwise model selection by AIC using the function ‘stepAIC’ from the package ‘MASS’. In another analysis, we looked for the best model to explain estimated bite force by mandibular centroid size, age and sex.

To test whether mandible or cranial shape are correlated to estimated bite force, we performed Procrustes ANOVA. The patterns of covariation between mandibular or cranial shape (block 1) and estimated bite force at the three bite points (block 2) were explored using 2B-PLS analyses.

To compare the variability in cranial and mandible shape between dogs and foxes, we performed a PCA on the Procrustes coordinates of the merged coordinates from this study and those of previous studies using the same landmarking protocol (Brassard et al., 2020b,c). To quantify and compare the level of morphological variation between the two species (called disparity), we used the function ‘morphol.disparity’ from the package ‘geomorph’ (Foote, 1993; Zelditch et al., 2012). Morphological disparity is estimated as the Procrustes variance in each species, using residuals of a linear model fit (the sum of the diagonal elements of the group covariance matrix is divided by the number of observations in the group). The differences between species were statistically evaluated through 10,000 permutations, where the vectors of residuals were randomised among groups.

We compared the proportions of the masseter, temporal and pterygoid muscles to the total mass of the adductor muscles as well as the relative contribution of the moment exerted by each of the adductor muscles on the moment of the estimated bite force between dogs and foxes using bilateral Welch's two-sample t-tests.

To test whether the integration was greater in foxes than in dogs, we compared the PLS coefficients obtained here from the 2B-PLS between cranial or mandibular shape (allometry-free or not) and absolute or residual muscle data or estimated bite force for the red fox with those obtained previously for dogs (Brassard et al., 2020a,b,c). For this purpose, we calculated Z-scores with the function ‘compare.pls’ from the package ‘geomorph’.

Model outputs are detailed in Table S1.

Visualisations of the variability in cranial or mandibular shape (or allometry-free shapes) in red foxes are provided in Fig. S1. The multiple Procrustes ANOVA (Table 1) show that size, age and sex explain 17.4% of the variation in cranial shape and 16.3% of the variation in mandible shape. Shape strongly correlates with size (cranium: N=58, R²=0.061, P<0.001; mandible: N=68, R²=0.055, P=0.002). Bigger individuals have a longer snout, a lower and smaller braincase, a more marked postorbital constriction, narrower but thicker and more anteriorly oriented zygomatic arches (Fig. S2A), more developed coronoid, condylar and angular processes on the mandibular ramus, and a straighter ventral border of the mandibular body (Fig. S2B). The effect of age was significant on the allometry-free shape of the mandible but not on that of the cranium (Table 1). The effect of sex on the residual shape (once the effects of size and age are removed) was significant for the cranium only (Table 1).

The pennation angles are around 0 deg in the digastric, 30–40 deg in the temporalis and masseter, and 40 deg in the pterygoids. Muscles from the temporalis complex have very long muscle fibres (from 21±5 mm for the SZ to 24±5 mm and up to 50 mm for the TS) compared with muscles from the masseteric and pterygoid groups that have shorter fibres (from 9.6±1.5 mm for the ZMP to 15.5±8 mm for the MS, and 11.5±3 mm for the PM; see Table S1 for details). The mass of the lateral pterygoid muscle represents only around 9% of the mass of the pterygoid complex (from 0 to 25%) and 0.77% of the total mass of the adductor muscles (from 0 to 2.5%). The masseter represents 27±4%, the temporalis 64±8% and pterygoid muscles 9±3% of the total mass of the adductor muscle. Muscle masses were strongly correlated (r>0.8 for all groups).

Visualisations of the variability in jaw muscle architecture are reported in Fig. S3A. The first axes were strongly correlated with mandible centroid size, age and sex (P<0.001 in all cases, N=58). Silver foxes were included within the variability of the red foxes. However, it can be noticed that they plot with the youngest red foxes, on the left part of the scatterplots with the scaled masses (Fig. S3B) and scaled PCSAs (Fig. S3D). This suggests that these four adult/old silver foxes have smaller MS, TS, TP and P than the average of all the foxes of our sample.

A comparison of bite force estimated on the incisor and carnassial teeth at a gape angle of 20–30 deg (and for an AFRF of 90 deg) with the in vivo data showed good correspondence (Fig. 3; Table S1). On the incisor teeth, the mean bite force estimated for a gape angle of 20 deg (206±55 N) was not significantly inferior to the mean of the in vivo bite force (243±67 N; Welch's unilateral test P=0.06), whereas model outputs for a gape angle of 30 deg (191±53 N) slightly underestimate the in vivo bite force (Welch's unilateral test P=0.02). The mean in vivo bite force at the carnassial teeth (337±86 N with a maximum of 484 N in vivo) was significantly lower than the mean of the model outputs for a gape angle of 20 deg (434±111 N; Welch's unilateral test P=0.003) and for a gape angle of 30 deg (403±107 N; Welch's unilateral test P=0.02). However, we could not compare the two methods for the same individuals, except for the specimen Ny-R5, which was included in both the model and the in vivo measurements and further showed that estimated bite forces are very close to the maximal forces recorded in vivo (incisor teeth: 183 N in vivo and 192 N in our model estimation; molar teeth: 435 N in vivo and 408 N in our model estimation for a gape angle of 20 deg and an AFRF of 90 deg). This suggests that the model output gave a reliable estimate of in vivo data.

As expected, the mean bite force estimated at the carnassial tooth was higher than the one estimated on the canine or incisor teeth (Table 1; Fig. S3). For example, at a gape angle of 0 deg for an angle of the food reaction force of 90 deg, mean estimated bite force ranged from 174 N (gape angle 40 deg) to 230 N (gape angle 0 deg) at the incisor teeth, from 198 to 261 N at the canine tooth and from 368 to 486 N at the carnassial tooth. For a given AFRF, mean estimated bite force decreased when gape angle increased (Table 2). A shift of the food reaction forces away from the perpendicular axis caused an increase in estimated bite force (Fig. S3), which is consistent with previous observations (Dumont and Herrel, 2003). In contrast, the mean joint force increased when the gape angle increased and when the point of application of the food reaction force was closer to the incisor teeth, ranging from 534±123 N on the incisor teeth and 463±109 N on the carnassial tooth for a gape angle of 0 deg, to 553±128 N on the incisor teeth and 484±114 N on the carnassial tooth for a gape angle of 40 deg (Table 2). The more elevated the angle of the food reaction force, the more elevated the force in the joint (Fig. S5). The angle of the joint force decreased when the gape angle increased and when the point of application of the food reaction force was closer to the incisor teeth, ranging from 140±5 deg on the incisor teeth and 152±7 deg on the carnassial tooth for a gape angle of 0 deg, to 131±3 deg on the incisor teeth and 138±4 deg on the carnassial tooth for a gape angle of 40 deg (Table 2). This aligns the joint force more with the orientation of the joint capsule. The more elevated the angle of the food reaction force, the more elevated the force in the joint.

The shape of the cranium strongly covaried with that of the mandible (PLS1: 30% of total covariance, coefficient of covariation r_PLS=0.78, P=0.02, Fig. 4A; PLS2: 19% of total covariance, r_PLS=0.69, P<0.001; PLS3: 14% of total covariance, r_PLS=0.73, P<0.001; PLS4: 9% of total covariance, r_PLS=0.74, P=0.002; PLS5: 6% of the total covariance, r_PLS=0.61, P=0.003). A mandible with a large body that remains straight towards the anterior end with a more caudally inclined coronoid process and a smaller angular process is associated with a shorter and higher cranium, lower and slightly larger zygomatic arches, a larger braincase, more anteriorly positioned orbital processes and a more oblique snout (Fig. 4A). These patterns of covariation partly match our observation of the deformations along the allometric slopes (Fig. S2). Indeed, linear regressions performed on the PLS1 scores of each block and the log₁₀ of the centroid size indicate that covariation was driven by the centroid size of the mandible (R²=0.09, P=0.02) and cranium (R²=0.13, P=0.003). These results support the results of the multiple regression analyses: most of the shape change is driven by size (and possibly age). The covariation between the allometry-free shapes of the cranium and mandible was significant for secondary axes (PLS1: 33% of total covariance, r_PLS=0.79, P=0.054; PLS2: 23% of total covariance, r_PLS=0.66, P<0.001) and reveal that straight mandibles with large muscle insertions are related to relatively shorter but wider crania, with greater muscle insertions (more space available for the temporal muscle to pass through; Fig. 4B).

The multiple Procrustes ANOVA showed that muscle PCSA and mass explained 10.3% of the variation in cranial shape and 11.5% of the variation in mandible shape once size, age and sex were taken into account (Table 1). The simple regressions indicate that cranial shape is more closely associated with the relative volume occupied by the temporalis and masseter muscles than to the PCSA (which takes into account the arrangement of the fibres in the muscle) of these muscles (Table S1). Muscle PCSA did not predict variation in cranial shape, however. In contrast, mandible shape was associated with both the volume and PCSA of the temporalis and masseter muscles.

Because muscle data were strongly correlated, an increase in the PCSA of the masseter, temporal or pterygoid muscles was associated with a similar variation in shape for the upper jaw as well as for the lower jaw (Fig. 5). The PCSA of the masseter, temporalis and pterygoid muscles was related to the area of insertion of the three muscles: the dorsal tip of the coronoid process, the masseteric fossa and the angular process. The shape of the braincase seems to be more related to variation in the PCSA of the temporalis muscle. An increase in the PCSA was related to a change in the convexity of the temporal bones and the shape of the sagittal crest. The PCSA of the masseter drives the shape of the zygomatic arch more specifically, although the shape of the braincase is also impacted. The pterygoid bone does not seem to be impacted much by the PCSA of the three main adductor muscle groups.

Additionally, the 2B-PLS showed significant covariation between muscle data (scaled or not) and the shape of the mandible, and irrespective of whether allometries were considered or not (Table 3). The same observations can be made for the covariation between muscle volume and the shape of the cranium. However, there was no significant covariation between raw muscle PCSA and cranial shape, but the covariation was significant for scaled PCSA and/or allometry-free shape. Ontogeny/age seems to play a major role in the intensity of the covariations. However, even after removing the youngest individuals (four foxes classified as juveniles), covariation remained significant and strong. Covariation was significantly less important between muscle data and the ramus of the mandible only than for the complete mandible (mass: Z=2.07, P=0.02; PCSA: Z=1.17, P=0.01). There was no significant difference between the covariations obtained for the mandible and those for the cranium (P>0.05 in all cases). The covariations drastically decreased when shape and/or muscle data were scaled, which suggests the strong importance of size in the patterns of integration.

The outputs of the biomechanical model show that the temporal represents around 50% of the total moment of the estimated bite force, the masseter 40% and the pterygoid 10% (Table S1). When the gape angle increases, the contribution of the masseter decreases, while that of the pterygoid (and the temporalis) increases. Detailed contributions of all muscle bundles and for several gape angles are reported in Table S1. For example, for a gape angle of 0 deg, M. masseter superficialis contributes to 16% of the moment of the estimated bite force, while M. temporalis superficialis contributes to 24%.

The best model explaining variation in estimated bite force with mandibular centroid size and muscle data (PCSA, mass, fibre length, pennation angle of the three main muscle groups) is the model that considers the PCSA of the masseter (0.40), the PCSA of the temporalis (0.40) and the PCSA of the pterygoids (0.14) and the pennation angle of the masseter muscle (0.12, adjusted R²=0.82, P<0.001).

The centroid size of the mandible and age are also important drivers of estimated bite force as suggested by the results of the linear regression of estimated bite force with size and age (R²=0.60, P<0.001). Moreover, males have significantly higher estimated bite forces (Welch’ two-sample t-test P=0.025) than females, probably because of their larger size. There was no significant difference between sexes in estimated bite force relative to size (of the mandible or cranium, Welch’s two-sample t-test P=0.9). The Procrustes ANOVA (Table 1) showed that 6–7% of the variation in mandible shape is related to estimated bite force or residual estimated bite force (P<0.001).

We also observed significant covariation between mandibular shape and estimated bite force (r_PLS=0.64, P<0.001) or residual estimated bite force (r_PLS=0.56, P=0.002) and between the allometry-free mandible shape and the residual estimated bite force (r_PLS=0.56, P<0.002; Table 3). Once the juvenile foxes were removed, the covariation with residual estimated bite force was still significant but lower (r_PLS=0.50, P<0.01; Table 3). The first PLS axis of the 2B-PLS between mandible shape and estimated bite force (Fig. 6A) was strongly related to the size of the individuals (P<0.001), and consequently also the age of the individuals (P<0.001). Foxes on the left part of the scatterplot (smaller and/or younger) had a proportionally shorter ramus and a longer body, which is more ventrally rounded, a smaller and straighter coronoid process, a small angular process, and a thick body under the carnassial tooth. These foxes had low estimated bite forces. Foxes with high bites forces, on the right part of the scatterplot, had a proportionally large coronoid process and a reduced body length, a more caudally oriented coronoid process with a deeper masseteric fossa, bigger angular and condylar processes, and a more angular and ventral border of the body. Shape deformations were overall similar in the 2B-PLS between mandibular shape and the residual estimated bite force (Fig. 6B). In contrast, the Procrustes ANOVA and 2B-PLS analyses performed with cranial shape and estimated or residual estimated bite force showed that there was no significant correlation or covariation between them (Tables 1 and 3).

The PCAs describing the variation in cranial and mandibular shape for both dogs and the foxes from this study show that dogs and foxes are clearly separated along the first PC axis (accounting for 32% of the total variance for the mandible and 54.6% for the cranium; Fig. 7). Almost all the foxes are located on the right part of the scatterplot. They have straight and flat mandibles, with a small and triangular coronoid process, and a low, long and straight cranium, in contrast to dogs, which have a more curved body and a more rounded and larger cranium, with a reduced snout. The segregation is stronger for the cranium than for the mandible as a few were very close or even overlapped the morphological space of dogs in terms of mandible shape. In particular, the Dobermann, a relatively dolichocephalic dog, has a ‘fox-like’ mandible. The results of the disparity tests indicate that the shape disparity of the cranium is greater in dogs than in red foxes (Procrustes variance: 0.0094 in dogs versus 0.0015 in foxes, P<0.001). The same applies to the shape disparity in the mandible (Procrustes variance: 0.0064 in dogs versus 0.0040 in foxes, P<0.001), but the difference between species is smaller.

The proportions of the masseter, temporalis and pterygoid muscles to the total mass of the adductor muscle were similar in foxes and dogs (in foxes: respectively 27±4%, 64±8%, 9±3%; dogs: 27±5%, 63±10%, 10±3%; Welch's t-tests P>0.10 for each muscle group). MANOVA showed that there was a difference in the PCSA of the adductor muscles of the foxes compared with dogs, which was mostly explained by size (P<0.001 on raw data, P>0.05 on residuals). The first axis of the PCA combining data for dogs and foxes shows that the variability in muscle architecture of foxes clearly overlaps with that of dogs. Values for the red fox were similar to those for small dogs and beagles (Fig. S4A). The disparity tests further indicate that the raw muscle masses and PCSAs were more variable in dogs than in foxes (variance for mass: 0.93 in dogs versus 0.25 in foxes; for PCSA: 0.54 in dogs versus 0.22 in foxes; P<0.001). Residual volume and residual PCSA (residuals of the regression with the log₁₀ of the mandibular centroid size) cover the same area in foxes as described by dogs which suggests a great variability in the relative importance and strength of the jaw muscles in red foxes despite a lower shape variability (Figs S3B and S4B). Tests further indicate that the difference in variability was less significant for residual mass (variance for mass: 0.26 in dogs versus 0.16 in foxes, P=0.02) and no longer significant for residual PCSA (variance for mass: 0.25 in dogs versus 0.19 in foxes, P=0.051).

The patterns described by estimated bite force, joint force or the angle of the joint force according to the angle of the food reaction force are similar in foxes to the ones described previously for dogs (Fig. S5). As for dogs, the muscles that contribute most to the estimated bite force are the temporalis (50%) and then the masseter (40%). Yet, the MS and TS contribute proportionally more to the estimated bite force in red foxes than in dogs (for a gape angle of 0 deg, for the MS: 16% in the fox, 13% in dogs, Welch's bilateral t-test P<0.001; for the TS: 24% in the red fox, 22% in dogs, Welch's bilateral t-test P<0.01). The Z-scores show that the covariation between mandible shape and estimated bite force was significantly lower in foxes than in dogs (foxes: r_PLS=0.63; dogs: r_PLS=0.75; P associated with the Z-score, P_Z=0.002). The same was observed with residual estimated bite forces (foxes: r_PLS=0.55; dogs: r_PLS=0.65; P_Z=0.035) and allometry-free shapes (foxes: r_PLS=0.55; dogs: r_PLS=0.69; P_Z=0.02).

In this study, using three-dimensional geometric morphometrics and jaw muscle dissection, we evidenced substantial morphological variability in cranial shape (which was associated with size and sex), mandible shape (which was associated with size and age), muscle architecture and estimated bite force (using a 3D static biomechanical model based on dissection data) in red and silver foxes.

Our study showed that the masticatory system is highly integrated in foxes: the shape of the mandible and cranium were strongly correlated and were both strongly correlated with muscle mass. Differences in shape and muscle architecture resulted in a wide range of estimated bite forces that probably offer different possibilities of adaptation according to the ecological context (e.g. more or less commensal). However, only mandible shape appeared to correlate with muscle cross-sectional area and estimated bite force, indicating that mandibular shape is a better predictor of bite force than cranial shape. We also evidenced that there was less variability in the cranium or mandible of foxes compared with that in dogs, yet the difference between species was much smaller for the mandible. Despite much greater variation in skull shape and size in domestic dogs, variation in muscle PCSA relative to size was equally great in foxes and dogs and we observed similar patterns of covariation between bone shape and muscles or estimated bite force. Future research is needed to investigate in greater detail the effect of the environmental variation on shape, muscles and estimated bite force in domestic, commensal and wild canids.

Our biomechanical model showed excellent correspondence with the in vivo measurements, especially for the specimen where we had both in vivo and model data. Several parameters may account for the slight differences we observed between the mean of model outputs and the mean of the in vivo bite forces. The in vivo bite forces were recorded for silver foxes and not red foxes. Further, we could not precisely control the gape angle (between 20 and 30 deg). This highlights the inter-individual differences that exist and the importance of deriving individual-specific models based on dissection data to be able to accurately estimate bite forces (Gröning et al., 2013). The estimated bite force ranged from 200 N on the incisor teeth to 450 N on the carnassial tooth for a gape angle of 20 deg and an AFRF of 90 deg. Higher estimated bite forces were recorded/calculated on the carnassial tooth (compared with the estimated bite forces at the incisor or canine teeth), as a result of the shorter out-lever arm (because the bite point is positioned closer to the attachment sites of the adductors; Dumont and Herrel, 2003; Ellis et al., 2008; Ellis et al., 2009; Greaves, 2000; Greaves, 2002; Herrel et al., 2008; Spencer, 1998). The higher the gape angle, the lower the estimated bite force and the angle of the joint force, and the higher the joint force. The same has been observed for dogs and other species of mammals (Brassard et al., 2020a; Dumont and Herrel, 2003; Herrel et al., 2008; Kerr et al., 2017; Santana, 2016). There are no in vivo data available in the literature for wild red foxes, unfortunately. Forbes-Harper et al. (2017) estimated bite force in Australian red foxes using dry skulls only and found forces ranging from 170 to 342 N (mean: 239 N). Whereas our estimations are in the same range, differences between populations are likely, especially when comparing invasive with native populations. Future studies are needed to compare estimations obtained with the dry skull method and those obtained using muscle data obtained from dissections.

As expected, most of the variation in estimated bite force (81%) was explained by muscle data, which were used for the construction of the biomechanical model, rather than variation in lever arms driven by the shape of the bony elements. Size and age were also important drivers of estimated bite force (alone, they explained 60% of the total variation). Additionally, males produced significantly higher estimated bite forces (but not relative to their size), probably as a result of their larger size. This was demonstrated previously with a sample of over 300 Australian red foxes (Forbes-Harper et al., 2017).

The finding that foxes tend to have the largest moment about the temporomandibular joint axis produced by the temporalis is consistent with the fact that, like other carnivores, they need to produce high bite forces at high gape angles (Greaves, 1985; Slater et al., 2009). Interestingly, the contribution of the muscles to estimated bite force does not reflect the contribution in muscle mass, as the masseter contributes relatively more to estimated bite force for its volume compared with the temporalis. This is because of the muscle architecture: the M. temporalis underperforms because of its longer muscle fibres, while the M. masseter overperforms thanks to shorter fibres. This is consistent with the observations of Penrose et al. (2020) across different species of canids. The significant volume of the M. digastricus (9% of the total volume) is probably associated with the need for a fast jaw opening during prey capture (Curtis and Santana, 2018).

Estimated bite force is explained by mandible and cranial shape, but these factors explain relatively little of total shape variation (7% and 2%, respectively), unlike what has previously been shown in studies on shrews (Cornette et al., 2015a) and other vertebrates (Dollion et al., 2017; Fabre et al., 2014). As all muscles were strongly correlated, the deformations associated with the PCSA of all the muscular groups or estimated bite force were similar and involved the area of origin or insertion of the muscles on the bone (Fig. 5). In particular, deformations involved the tip of the coronoid process (insertion of the TS), the rostral border of the angular process (insertion of the SZ), the masseteric fossa (MP, ZMA and ZMP), the angular process (PM) and, on the cranium, the sagittal crest, the temporal fossa and the post-orbital constriction (M. temporalis). Areas of mechanical constraints were also highlighted (between the orbital processes and on the snout at the midline, and the ventral curvature of the mandible), which is consistent with observations of the distribution of stress related to intrinsic loads, as described for other canids (Slater et al., 2009).

Correlations between cranial shape and residual PCSAs were not significant, contrary to the correlations between cranial shape and residual mass (Table 1). Thus, skull shape seems to be more closely related to space constraints rather than reflecting the modelling of the bones to mechanical constraints imposed by the external muscle loadings, contrary to the mandible. The coefficients of covariation between muscle mass and shape were more elevated for the cranium than for the mandible. The same observation has previously been made for dogs and strepsirrhines (Brassard et al., 2020c; Fabre et al., 2018). This supports the hypothesis that, in foxes, the cranium presents a lesser degree of functional plasticity compared with the mandible, probably because it has to cope with additional functional demands, such as the protection of the sensory systems and brain (Fabre et al., 2014; Figueirido et al., 2011). The mandible is thus a better predictor of functional demands and estimated bite force than the cranium. It is also possible that there are stronger correlations between internal bone structure and cranium and muscle PCSA or estimated bite force because cortical thickness may be a better proxy of external loads than the overall shape (Bouvier and Hylander, 1981; Daegling and Hotzman, 2003; Slizewski et al., 2013), but this remains to be tested.

Size explained most of the variation in both cranial and mandible shape. Sexual dimorphism was evident in the allometry-free shape for the cranium, but not for the mandible. This suggests that the observed sexual dimorphism may not be related to feeding but rather to other functions as the cranium also protects the sensory organs and the brain, for example (Fabre et al., 2014; Figueirido et al., 2011; Radinsky, 1981; Santagati and Rijli, 2003). Age correlated with the allometry-free mandibular shape but not the allometry-free cranial shape, probably because the strong effect of size on cranial shape hides further effects of age. However, analyses including additional juvenile foxes would be necessary to investigate the influence of age more exhaustively. Size, age, sex and muscle architecture, however, together explained relatively little of the total variation in both cranial and mandible shape (only 28% of the total variation in both cranial and mandible shape). These results are not surprising considering that low coefficients of correlation are typically found in mammals, including humans (Toro-Ibacache et al., 2016). This suggests that other factors – possibly developmental factors – constrain the shape of the cranium and mandible (Drake and Klingenberg, 2010; Wayne, 1986).

As food mechanical properties are known to influence cortical bone modelling and remodelling (Bouvier and Hylander, 1981; Ionova-Martin et al., 2011; Lieberman et al., 2004; Ravosa et al., 2015, 2016; Scott et al., 2014a,b), diet is another parameter that is worth taking into account. Unfortunately, we had no information about the diet of the specimens we dissected, so we could not test whether it impacts bone shape (or estimated bite force). The relationships between diet and cranial shape and estimated bite force were investigated for Australian foxes using the dry-skull method (Forbes-Harper et al., 2017) and showed that diet does impact cranial morphology but only to a small degree. As the foxes used in this study were mostly from the countryside in south-western France, we could not test for differences between urban or more rural foxes. These parameters would be of interest to explore in future studies.

Dogs and foxes clearly differ in both cranial and mandibular shape and size. Moreover, the disparity in shape is much lower in the fox than in dogs, which is particularly remarkable for the cranium. This is consistent with previous results showing that domestication has resulted in an increase in shape variability in dogs (Darwin, 1868; Drake and Klingenberg, 2010; Heck et al., 2018). The increased variability in mandibular shape in dogs is significant but less obvious, however (Fig. 7).

Interestingly, foxes exhibited a similar variability in muscle architecture once size was removed from the analyses (residual mass and PCSA) compared with dogs and it nearly spanned the same morphological range in foxes as in dogs. In foxes, the lateral pterygoid represents, on average, 1% of the total volume, which is consistent with previous observations in foxes (Penrose et al., 2020: 0.27% for one individual only) or even dogs (Brassard et al., 2020b: 0.67±0.43%, N=48). The proportions of the masseter, temporalis and pterygoid muscles in foxes are similar to what was observed in dogs (Brassard et al., 2020b). Additionally, muscle contributions to estimated bite force in foxes are overall equivalent to what was observed in dogs. The fact that the superficial layers of M. masseter and M. temporalis are more developed and contribute more to the estimated bite force in the red fox than in dogs is consistent with the necessity of a small prey hunter (foxes) to close the jaws quickly. This is also in line with the long muscle fibres of the temporalis, allowing fast jaw closure. Foxes bite less hard than dogs on average, but this is mostly the result of their smaller size and the long jaw out-lever due to their relatively longer jaws.

As expected, cranial and mandibular shape strongly covaried in foxes (r_PLS=0.78) but the covariations were not significantly stronger than those observed for dogs (r_PLS=0.81), in contrast to our predictions. Overall, we observed patterns of covariation that were similar to those for dogs and that were possibly driven by muscle constraints. This covariation was probably driven by the muscles that link the upper and lower jaws, which is supported by the results of the Procrustes ANOVA (Fig. 5). Additionally, we even evidenced a surprisingly lower covariation between mandible shape and estimated bite force in foxes compared with dogs (whether size is taken into account or not), but this was probably due to the greater level of variability in mandible shape in dogs that favours stronger covariation. All these results support the hypothesis that the masticatory apparatus is strongly integrated in canids, and that domestication did not lead to a disruption of the functional links of the jaw system despite the increased variability in shape and the commonly observed malformations in dogs.

We thank the Veterinary school ONIRIS (Nantes, France) and the Laboratory for Rabies and Wildlife (Nancy, France) for providing fox heads for dissection. We are grateful to Manuel Comte, Mickaël Godet and Frédéric Lebatard for their help in managing specimens and their helpful discussions about the preparation of the skulls. We also thank Arnaud Delapré for his help with photogrammetry. We are very grateful to the two anonymous reviewers for their comments and their valuable contribution to the manuscript.