Bite force and its relationship to jaw shape in domestic dogs

The constituents of the masticatory system have been described in some detail in the domestic dog (Barone, 2010; Budras, 2007; Curth et al., 2017; Evans and DeLahunta, 2010; Hoppe and Svalastoga, 1980; Johnson, 1979; Miller et al., 1965; Penrose et al., 2016; Robins and Grandage, 1977; Thomas, 1979; Tomo et al., 1993). During mastication, the lower jaws (i.e. the mandibles) move up or down relative to the upper jaw (here we use this term to refer to the cranium and face, following the Nomina Anatomica Veterinaria nomenclature; International Committee on Veterinary Gross Anatomical Nomenclature, 2017) by rotation about the temporomandibular joint that receives the condylar process of the mandible. These movements are driven by contractions of the jaw adductors. Acting like a lever, the forces are transmitted to the teeth, generating the bite force (Kim et al., 2018). The macroscopic arrangement of muscle fibres (i.e. muscle architecture) directly determines muscle force production. A good overall measure of this architecture is the physiological cross-sectional area (PCSA), which takes into account muscle volume, fibre length, fibre type and pennation angle (Haxton, 1944).

The extraordinary variability in the size and shape of the head (Brassard et al., 2020; Coppinger and Coppinger, 2001; Drake and Klingenberg, 2010; Miller et al., 1965; Selba et al., 2019; Wayne, 1986, 2001), and jaw muscle architecture (Brassard et al., 2020) between dog breeds raises questions about the impact of this variability on the function of the masticatory system and bite performance. Differences in skull shape between breeds have been suggested to be associated with differences in jaw strength (Case, 2013) and bite force (Ellis et al., 2008, 2009) as the shape of the neurocranium drives the size of the jaw muscles and the length and shape of the jaws determine the out- and in-lever arms of the system.

A few studies have investigated bite force in domestic dogs using the dry-skull method or in vivo measurements (Ellis et al., 2008, 2009; Kim et al., 2018; Lindner et al., 1995). However, quantitative data on muscle architecture that could be used to improve these models are scarce. Ellis et al. (2008, 2009) used two-lever models or multivariate regression modelling to estimate bite forces (Ellis et al., 2008, 2009; Kim et al., 2018; Lindner et al., 1995). The jaw is modelled as a two-lever system: jaw muscle cross-sectional area of the major jaw-adducting muscles, and the moment arms (the perpendicular distance between the point of application of the force and the temporomandibular joint) of the muscles (in-levers) and of bite points about the temporomandibular joint (out-levers) are approximated from skull dimensions taken from photographs. Skull length and skull width are then considered as a proxy of shape and size. Estimations obtained using the equations provided by Kiltie (1984) and Thomason (1991) were used and adjusted by values recorded in vivo on 20 dogs of various breeds during stimulation of the m. temporalis and m. masseter under general anaesthesia (Ellis et al., 2008). Ellis et al. (2008, 2009) also established an equation from multivariate regression analysis to estimate bite force independently of any lever model, using cranial measurements and the body mass of the same dogs for which bite force was recorded in vivo. However, the authors did not consider the muscle cross-sectional area and the effective moment arms of the forces in their equations. Moreover, in these two dimensional (2D) methods, the PCSA of the temporal muscle is often underestimated, while that of the m. masseter and m. pterygoideus is overestimated (Davis et al., 2010). The regression model using body mass was based on only 20 dogs of different breeds. Using these equations, negative bite forces were obtained for small brachycephalic dogs (with a short and wide skull), which demonstrates that the equation is not appropriate when applied outside of the range of values for which it was developed. However, no studies to date have explored the covariation between bite force and bone shape.

In the present paper, we aimed to explore the diversity of bite force in dogs as well as the relationships between the three dimensional (3D) shape of the upper and lower jaws and bite force. To do so, we used a biomechanical model based on 3D lever mechanics using muscle data (fibre length, pennation angle, muscle mass and PCSA) and the 3D coordinates of origin and insertion of jaw adductor muscles obtained from dissection of 47 dogs of various breeds. We also report in vivo measurements recorded from three trained Belgian shepherd dogs (Malinois) to validate the output of our bite model. We used a combination of geometric morphometric techniques and comparative methods to: (1) assess the variability in bite force in dogs and test for differences between morphotypes; (2) test which components of the jaw adductor system are the best predictors of bite force; and (3) describe the pattern of covariation with the overall shape of the upper and lower jaws.

Dogs can be classified based on the shape of their head by using the cephalic index, a ratio between skull width and length (Helton, 2011; Koch et al., 2003). There are three morphotypes: dolichocephalic (relatively long skulls) are opposed to brachycephalic (broad skulls) dogs, and mesocephalic (moderate skulls) dogs are intermediate. As a result of selection determined by standards (Fédération Cynologique Internationale, FCI; http://www.fci.be/en/Presentation-of-our-organisation-4.html), different breeds contain animals with specific traits/metrics that can therefore be assigned to one of these three morphotypes. Excessive artificial selection has resulted in some ‘hypertypes’ (where some characters within a dog breed are developed to excess; Triquet, 1999; Guintard and Class, 2017) among brachycephalic or dolichocephalic groups, showing exaggerated morphotypes. Given that our sample is small and gathers non-pure-breed dogs, we here focused on the impact of the morphotype (Roberts et al., 2010) only and compared brachycephalic with mesocephalic and dolichocephalic dogs. As previously stated by numerous authors (Ellis et al., 2008, 2009; Kim et al., 2018; Lindner et al., 1995), the combined variability in size and morphology – pertaining to both skull shape and jaw muscle architecture – among morphotypes probably significantly explains the variability in estimated bite force. For example, as suggested by Ellis et al. (2009), we expected bite forces to be higher in large brachycephalic dogs than in other morphotypes. Here, we aimed to describe the relationships between the overall morphology and bite force.

Moreover, we expected intensive breeding for aesthetic reasons or functional ability to potentially have perturbed the functional relationships between the different components of the feeding system as diet no longer imposes constraints on the jaw system. Indeed, the domestication of dogs has led to a release from ecological constraints, which may have increased the diversity in a large array of genes as aberrant phenotypes were no longer selected against (Björnerfeldt et al., 2006). Recent dog breeds largely feed on processed food requiring little or no chewing. Given that bite force is usually a good indicator of dietary diversity as it directly determines prey size and feeding ecology (Aguirre et al., 2002; Cornette et al., 2013, 2015; Dollion et al., 2017; Felice et al., 2019; Firmat et al., 2018; Forbes-Harper et al., 2017; Herrel and Holanova, 2008; Herrel et al., 2001, 2005; Huber et al., 2009; Kerr et al., 2017; Maestri et al., 2016; Marcé-Nogué et al., 2017; Nogueira et al., 2009; Sagonas et al., 2014; Santana et al., 2010; Van Daele et al., 2009; Verwaijen et al., 2002; Young and Badyaev, 2010), one would expect a disruption between bite force and bone shape, resulting in low coefficients of covariation (or possibly non-significant coefficients) across dogs as a whole.

The dataset is composed of 47 dog heads (Table 1). Breeds were estimated based on morphological similarities with dogs from existing standards. Beagles are the most represented, with 10 specimens. Given that most of the breeds are represented by only one specimen and that most of the dogs are likely to be cross-breeds, our sample does not allow any conclusion at the breed level.

To test for the effect of the morphotype (brachycephalic, mesocephalic or dolichocephalic), the cephalic index (CI) was calculated following Roberts et al. (2010): skull width/skull length×100 (Fig. 1). Skull width was measured between the two zygomatic arches (landmark 37 and the symmetric landmark to the sagittal plane), and skull length was measured from the anterior tip at the end of the suture of the nasal bones (landmark 2) to the most posterior point on the occipital protuberance (landmark 14). The brachycephalic dogs have the most elevated values of CI, the dolichocephalic dogs have the lowest values, and the mesocephalic dogs have intermediate values. Given that there is no clear consensus on the boundary between groups (Roberts et al., 2010), we chose limits to ensure that the three groups were similar in size. Specimens within the same breed can thus be classified into two adjacent morphotypes because they fall on either side of a morphotype boundary. This may be linked to age-related changes that can push skulls from one type to another. Dogs with a CI <0.70 were considered brachycephalic and the dogs with a CI <0.60 were considered dolichocephalic; dogs with an intermediate CI were mesocephalic.

Age was estimated based on the aspect of the cranial sutures (degree of closure), tooth wear and bone texture. Group A corresponds to dogs with molar teeth still erupting, a very porous mandible (minute interstices were observable with the naked eye) and unclosed cranial sutures (4–6 months according to Barone, 2010). Group B corresponds to dogs with the basispheno-basioccipital suture still open (<8–10 months according to Barone, 2010) and a still-porous mandible. Group D corresponds to old dogs with a closed interfrontal suture and worn denture (>3–4 years). Group C corresponds to intermediate adults (from 10 months to 3 years). We chose to keep the youngest individuals in our analyses to increase the morphological variability in the sample but they are far under-represented (only one dog in group A and one dog in group B) More detailed information about the sample is given in Table S1.

3D geometric morphometric analysis was used to describe the patterns of morphological variation. R version 3.6.0 (2019-04-26; http://www.R-project.org/) was used for all statistical analyses.

For the lower jaw, we considered 25 landmarks, 190 sliding semi-landmarks on curves and 185 sliding semi-landmarks on surfaces that were placed on 3D reconstructions of the lower jaw obtained using photogrammetry (Brassard et al., 2020; the 3D models of the mandibles are available on request from the corresponding author). A 3D sliding semi-landmark procedure was performed to transform all the landmarks into spatially homologous landmarks (Bookstein, 1997; Gunz et al., 2005). For the upper jaw, 54 landmarks were recorded on one side (left or right) of the upper jaw using a MicroScribe MX (Revware). A mirror function was then applied to obtain the symmetrical landmarks relative to the sagittal plane using the function ‘mirrorfill’ from the package ‘paleomorph’. This resulted in a total of 108 landmarks. Fewer landmarks were used for the upper jaw because the shapes are more easily quantified with a smaller number of landmarks. The landmarks are represented in Fig. 1 and described in Table S2.

For further visualization of shape changes, we used the 3D models of the upper and lower jaws of a beagle, obtained using photogrammetry. One-hundred and forty photographs were taken using a Nikon D5500 Camera (24.2 effective megapixels) with a 60 mm lens, by turning around the dorsal and ventral views of the upper jaw (Fau et al., 2016). One-hundred photographs per side were taken for the lower jaw. 3D models of the jaws were obtained after merging the two sides of each jaw, using the Agisoft PhotoScan software (©2014 Agisoft LLC, St Petersburg, Russia).

A generalized Procrustes analysis (GPA; Rohlf and Slice, 1990) was performed to obtain the Procrustes shape of each jaw using the function ‘procSym’ (Dryden and Mardia, 2016; Gunz et al., 2005; Klingenberg et al., 2002). The centroid sizes of each jaw were used as an estimation of size.

In vivo bite force data were recorded for three pure-breed Belgian shepherd dogs (Malinois) from a training club for defence dogs (Beauvais, France). The dogs are trained to bite for competitions and shows. These dogs are thus expected to bite relatively hard as a result of artificial selection by breeders for this purpose. We used a piezoelectric isometric Kistler force transducer (9311B, range ±5000 N; Kistler Inc., Winterthur, Switzerland) linked to a charge amplifier (type 5058A5, Kistler Inc.) similar to the one used for large turtles in Herrel et al. (2002). The transducer was mounted in a custom set-up and fixed on a wooden stick and covered by hessian fabric to protect the teeth of the dog and to provide a known bite substrate (see Movie 1). We performed several consecutive trials (at least 4) for each animal. The dogs did not all bite at exactly the same position on the transducer, so we had to correct the recorded bite force for each trial by taking into account the distance between the real location of the bite that we recorded in vivo and the sensor (corrected BF=recorded BF×distance between the sensor and the bite plates of the transducer/distance between the location of the bite and the sensor). For this purpose, we filmed each trial so we could extract the location of the bite relative to the sensor. The dogs bit at a gape angle of about 40–45 deg while grabbing the hessian cover with the last premolar teeth. We retained the maximal corrected bite force recorded across all trials for analyses. Head length and width were measured from photographs. The obtained values were used to validate the model output for shepherd dogs of a similar size biting on the first lower molar tooth at 40 deg.

To estimate bite force, we used a 3D lever model similar to that described by Herrel et al. (1998a,b). The movement of the lower jaw near occlusion is mainly rotational so we did not consider any translational movement. 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 bite force.

At static force equilibrium, for each side (left or right), the sum of the moments of the external forces is zero (positive moments of all the muscles on one side plus negative moments of the bite force on one side). The moment is a vector that corresponds to the vectorial product of the moment arm and the force. The magnitude of the moment thus corresponds to the product of force magnitude and the shortest distance between the centre of the system and the line of action of the force. The muscular forces were established by multiplying the reduced PCSA by a conservative muscle stress estimate of 30 N cm⁻² (Herzog, 1994). The maximal bite forces were then deduced from the sum of the muscle moments for each side, and doubled, considering that the adductor muscles on both sides are contracting maximally during maximal effort biting.

To calculate the effective length of the lever arm for all the muscle moments, we used the 3D coordinates of origin and insertion of each muscle (Fig. S1) that we recorded with a microscribe. We first chose a reference frame with a centre located at the right temporomandibular joint (Fig. 2). The x-axis runs through the rostral border of the right mandible, just medially to the first incisor, and the y-axis is directed towards the top of the skull and perpendicular to x. The approximate centroid of the origin and insertion areas of the muscles were used, based on observations from our dissections (Brassard et al., 2020). We also recorded the 3D coordinates of three possible points of application of the food reaction force: at the incisor teeth, at the lower canine tooth and at the lower carnassial tooth (point between the fourth upper premolar tooth P⁴ and the first lower molar tooth M₁). We chose these locations as they are important during feeding in canids. The first point is at the first incisor tooth (BPi in Fig. 2), the second one is just behind the lower canine tooth (BPc in Fig. 2) 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. 2). The last two points are compatible with the bite points chosen by Ellis et al. (2008, 2009).

The bite force vector is opposed to the food reaction force. Given that we do not know the direction of 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; Herrel et al., 1998b; Cleuren et al., 1995), we calculated the moment of the food reaction force for a large range of angles thereof (set to vary between −40 and −140 deg with respect to the lower jaw). We calculated the bite forces for several mouth opening angles (0, 20 and 40 deg). The magnitude and orientation of the joint forces were estimated as well because, at static equilibrium, the sum of the external forces (muscle and 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 bite force. Model output consists of the magnitude of the 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 bite force is available on request from the corresponding author.

As previously suggested for muscle data, bite force is expected to be highly dependent on size (Brassard et al., 2020). To test whether skull length or the centroid size of the jaws is a driver of bite force, we performed linear regressions on bite force and log₁₀ of skull length or log₁₀ of centroid size of the upper or lower jaw using the function ‘lm’. The importance and significance of the correlation between bite force and the centroid size of the bone were explored using the function ‘cor.test’. The residuals of these two regressions are further considered as ‘residual bite forces’.

To test for differences in bite force between morphotypes (brachycephalic, dolichocephalic and mesocephalic), an ANOVA and a linear model were calculated on residual bite forces using the functions ‘anova’ and ‘lm’. Post hoc tests were performed using the function ‘TukeyHSD’.

To investigate the muscular drivers of bite force (to determine which muscle contributed the most to the variation in bite force among the muscle bundles we dissected), we performed a linear regression of the bite force on the mandibular centroid size and main muscle mass, fibre length, pennation angle and PCSA using the function ‘lm’. For this analysis only, we considered the three main adductor complexes to increase statistical power and avoid noise in the data: the masseter complex, the temporal complex and the pterygoid complex. Among each complex, fibre length and pennation angle were averaged while mass and PCSA were summed. Data were log₁₀-transformed before analyses. We considered the calculated bite force for a food reaction force orientation of 90 deg and a gape angle of 20 deg. The best-fitted models were obtained from stepwise model selection by AIC using the function ‘stepAIC’ from the package ‘MASS’. To compare the contribution of each bundle to the bite force, we calculated the ratio of the moment exerted by each muscle and the moment exerted by the bite force at the lower carnassial tooth. We performed Friedman tests and post hoc tests using the functions ‘friedman.test’ and ‘posthoc.friedman.nemenyi.test’ from the package ‘PMCMRplus’ to compare the contribution of each bundle for the four gape angles and to compare the contribution of the three main muscular complexes for a gape angle of 0 deg.

To test whether the shape of the upper or lower jaws is a driver of bite force, we performed Procrustes ANOVA on jaw shape and bite force or residual bite force using the function ‘procD.lm’ from the package ‘geomorph’ (Adams and Collyer, 2016, 2017; Anderson, 2001; Anderson and Braak, 2003; Collyer et al., 2015; Goodall, 1991).

We explored the covariations between bite force (block 1) and the Procrustes coordinates of the upper or lower jaw (block 2). The patterns of covariation were explored using two-block partial least square (2B-PLS) analysis with the function ‘pls2B’ from the package ‘Morpho’ (Rohlf and Corti, 2000). The 2B-PLS method constructs pairs of variables that are linear combinations of the variables within each of the two blocks and that maximize the covariance between blocks (Rohlf and Corti, 2000). With this method, PLS coefficients and P-values are generated. PLS coefficients are the coefficient of correlation between PLS scores (between blocks) and thus reflect the intensity of the covariation (we refer to these coefficients as coefficients of covariation, r-PLS, below). P-values reflect the significance of the covariation for each new axis. They are calculated by comparing the singular value with those obtained from 10,000 permuted blocks. We did not consider phylogeny (Parker et al., 2004) in our analyses because we had no indication of a pure breed membership.

For these analyses, bites force and muscle data were log₁₀-transformed. 2B-PLS analysis was thus conducted between the shape of the upper or lower jaw and bite force or residual bite force. We also performed analysis for brachycephalic dogs only, and for mesocephalic/dolichocephalic dogs only. As differences in the number of variables and the number of individuals influence the PLS coefficient, a Z-score was calculated to compare the levels of functional integration between different types of dogs with the function ‘compare.pls’ from the package ‘geomorph’ (Adams and Collyer, 2016). The deformation of the mandible of a beagle to the consensus of the GPA was used as a reference for all visualizations. The beagle was chosen because it was the dog that was closest to the centre of the PCA describing variation in mandibular shape in our previous study (Brassard et al., 2020).

The model outputs for all specimens are detailed in Table S1. Although the small sample size did not allow us to describe the intra-breed variability, except for beagles, we have indicated the breeds in the Results so that future studies can expand upon our results.

The magnitude of the bite force ranged widely depending on the gape angle, bite point and orientation of the food reaction force (Table 2, Fig. 3). Mean bite force decreased when the gape angle increased. Moreover, the mean bite force estimated on the carnassial tooth was more elevated than that estimated on the canine tooth or incisor teeth. For example, for an angle of the food reaction force of 90 deg, mean bite force ranged from 412 to 531 N at the incisor teeth, from 478 to 617 N at the canine tooth and from 846 to 1091 N at the carnassial tooth, for a gape angle ranging from 40 to 0 deg. A shift of the food reaction force away from the perpendicular axis caused an increase in bite force. Contrary to mean bite force, mean joint force increased when the gape angle increased (Table 2), for all orientations of the food reaction force. For example, for an angle of the food reaction force of 90 deg, it ranged from 1018 to 1065 N at the incisors, from 990 to 1040 N at the canine tooth and from 863 to 912 N at the carnassial tooth for a gape angle ranging from 40 to 0 deg. The angle of the joint force decreased when the gape angles increased (Table 2), ranging from 131 to 139 deg for incisors, from 132 to 141 deg for the canine tooth and from 140 to 154 deg for the carnassial tooth.

For a given gape angle and bite point, a great variation in bite force exists among dogs (Fig. 4). We observed similar patterns of variation for all gape angles. To merge in vivo measurements recorded on the Malinois shepherd dogs with the bite forces estimated from dissection in the following descriptions, we focused on the model outputs obtained for a gape angle of 40 deg and we considered an angle of the food reaction force of 90 deg. Estimated bite force ranged from 100 to 1092 N on the incisor teeth, from 116 to 1268 N on the canine teeth and from 214 to 2172 N on the carnassial tooth. The three in vivo measurements were included in the overall variability of the predicted bite force. In vivo values were relatively close to our estimations for other dogs that are similar in shape (dog 1: 1094 N, dog 2: 688 N, dog 3: 903 N). The in vivo measurements thus validate our model. The dogs with the highest bite forces in our sample were the largest brachycephalic dogs, such as the rottweiler (2172 N on the carnassial tooth for one individual) and the pitbull (2051 N). The dogs with the lowest bite forces were the smallest dogs, belonging to the toy group. If we consider the 10 beagles we dissected, calculated bite forces ranged from 262 to 466 N on the incisor teeth (mean: 359 N), 302 to 481 N on the canine tooth (mean: 375 N) and 501 to 902 N on the carnassial tooth (mean: 709 N). Bite force was correlated to the length of the skull (R²=0.33, P<0.001), as well as to the mandibular centroid size (R²=0.54, P<0.001) and that of the upper jaw (R²=0.41, P<0.001). Brachycephalic dogs produced higher bite forces than dolichocephalic and mesocephalic dogs when scaled to the same skull length (P<0.05 when testing for differences between brachycephalic and mesocephalic or dolichocephalic dogs), which suggests that the shape of the upper jaw is an important driver of bite force and that brachycephalic dogs produce higher bite forces. However, the Leonberger dog from our sample seems to break this trend, as it produced a bite force that was as high as that of large brachycephalic dogs. There was no significant difference between mesocephalic and dolichocephalic dogs (P>0.05).

The contributions of each adductor muscle to the total moment of the bite force are indicated in Table 3 and Fig. 5. The muscles that contributed the most to bite force are the m. temporalis pars superficialis (TS; 22%) and m. temporalis pars profunda (TP; 25%). The moment exerted by the temporal complex (SZ, TS and TP) was, on average, responsible for 50% of the moment of the bite force, while the m. masseter (MS, MP, ZMA and ZMP) was responsible for around 37% and the pterygoid (P) 13% (P<0.001 between all muscle groups according to the post hoc tests). We noticed that the more important the angle of mouth opening, the higher the contribution of the m. temporalis pars profunda (TP) and m. zygomaticomandibularis anterior (ZMA), and the lower the contribution of the m. masseter pars superficialis (MS) and m. temporalis pars suprazygomatica (SZ), which was supported by significant results of the post hoc tests after a Friedman test (P<0.05). The contribution of the m. pterygoideus (P) remained almost unchanged. The contributions of the m. temporalis pars superficialis (TS), m. masseter pars profunda (MP) and m. zygomaticomandibularis posterior (ZMP) did not significantly differ.

The centroid size of both the upper and lower jaw is a driver of bite force (upper jaw: R²=0.41, P<0.001; lower jaw: R²=0.54, P<0.001).

Stepwise multiple regressions with calculated bite forces at the carnassial tooth as a dependent variable and muscle mass, fibre length and PCSA as independent variables retained a significant model with mandible size (β=−0.40, P=0.02), m. masseter mass (β=0.62, P<0.001), m. temporalis (β=−0.18, P=0.08), m. pterygoideus (β=0.18, P=0.09) and m. masseter fibre length (β=−0.35, P<0.001), m. pterygoideus fibre length (β=−0.19, P=0.002), m. temporalis PCSA (β=0.44, P<0.001), m. masseter pennation angle (β=−0.10, P<0.2) and m. pterygoideus pennation angle (β=−0.10, P=0.04) as best predictors (R²=0.97; P<0.001).

The Procrustes ANOVA showed that the shape of the upper jaw was not correlated to the absolute bite force, while the absolute bite force explained 16% of the variation in the shape of the lower jaw (P<0.001). However, the residual bite force was significantly correlated to the shape of both the upper jaw (R²=0.17, P<0.001) and the lower jaw (R²=0.11, P<0.001).

A summary of the results of the 2B-PLS analysis is given in Table 4. Only the main results will be described further. Further visualizations are presented in Fig. S2.

We observed significant covariation between the shape of the upper and lower jaws and absolute bite force (Table 4; Fig. S2). Significant covariation with high coefficient of covariation was also observed between residual bite force and bone shape (r-PLS=0.65 for the lower jaw versus 0.62 for the upper jaw; Fig. 6). A comparison of the Z-scores indicated that there was no significant difference in the level of covariation between the upper and lower jaws and between absolute and residual bite force (P>0.05). The covariation remained significant if we distinguished brachycephalic from mesocephalic/dolichocephalic dogs in the analyses with absolute bite force, but was no longer significant in the analyses with residual bite force. This suggests that size and the diversity in shape are important drivers of the covariation.

For all analyses, the first axis of the 2B-PLS explained more than 99% of the total covariance. The large brachycephalic dogs occupied the left part of the scatterplot and were related to high (or relatively high) bite forces. The distinction between brachycephalic and mesocephalic/dolichocephalic dogs along the first PLS axis was even clearer for the PLS with residual bite force. This reinforces the idea that the ability of large brachycephalic dogs to produce high bite forces is related to the specific shape of both the upper and lower jaws. The beagles occupied a small part of the scatterplot, even though a small variability in the covariation was observed. Age did not seem to drive the observed covariation.

Below we describe only the covariation between bone shape and the residual bite force (Fig. 6), and we refer to Fig. S2 for covariation with the absolute bite force. Dogs that produce a low bite force for their size have an elongated, flat and straight mandibular body in the sagittal plane, a small and narrow coronoid process with a shallow masseteric fossa, a medially short and small condylar process of the mandible and a weak angular process (Fig. 6A). The upper jaw is fox-like (Fig. 6B): it is elongated, with a proportionally long, flat and narrow snout, and a reduced braincase. The zygomatic process of the frontal bone and the post-orbital constriction are not very pronounced and the zygomatic arches are narrower. The perpendicular plate of the palatal bone is reduced and the retro-articular process of the temporal bone is reduced. In contrast, the dogs that can produce a high bite force for their size have a very robust mandible with a relatively large, coronoid process with a deep masseteric fossa, a shortened, ventrally and laterally curved mandibular body, a big, medially extended and caudally curved condylar process of the mandible, and a more pronounced angular process. The upper jaw is more massive, with a proportionally shorter, wider and laterally very marked snout, and a bigger and more rounded braincase. The zygomatic arches are much larger and more distant from the cranium, and the area that bears the frontal process of the zygomatic bone is more craniodorsally elevated. The perpendicular plate of the palatal bone and the retro-articular process of the temporal bone are well developed. For both the upper and lower jaws, the cheek teeth (premolar and molar teeth) are more cranially located for the dogs that produce the highest bite forces.

The objectives of this study were: (1) to assess the variability in bite force in the domestic dog considering a wide range of breeds/morphotypes; (2) to test which components of the jaw adductor system are the best predictors of bite force; and (3) to describe the pattern of covariation between bone shape and bite force.

In this study, we used a biomechanical model to explore the effect of the great variation in size and shape of the jaw on bite force in dogs. We provide the first estimations of bite force using individual PCSA and 3D coordinates of attachment of the jaw adductors obtained from dissection in a sample of domestic dogs of various breeds. These data are complementary and consistent with those already available in the scientific literature (values recorded in vivo or under anaesthesia, estimation using the dry-skull method or regression methods calibrated by measurements obtained under anaesthesia; Ellis et al., 2008, 2009; Lindner et al., 1995). For comparative purposes, we here focus on estimates for a gape angle of 30 deg. We estimated bite forces from 124 to 1380 N (mean: 520 N) on the canine tooth and from 229 to 2364 N (mean: 919 N) on the lower carnassial tooth. The strongest biters in our sample were the rottweiler and pitbull, with values exceeding 2000 N at the carnassial tooth. Lindner et al. (1995) provided in vivo measurements on 22 dogs of various breeds ranging from 13 to 1394 N (there is no mention of the bite point). Ellis et al. (2008) measured values under anaesthesia on 20 dogs of various breeds ranging from 147±6.9 to 926±8.1 N on the canine tooth and from 574±83.2 to 3417±43.1 N on the lower carnassial tooth. Ellis et al. (2009) estimated bite forces up to 4468 N on the lower carnassial tooth with their regression models and up to 3338 N with the dry-skull method as described by Thomason (1991). It is difficult to ascertain whether our estimates are more accurate than those of Ellis et al. (2008, 2009). Their studies were based on a simple 2D analysis of bite force supported by experimental determination of bite force on 20 dogs of a range of shapes and breeds. The present work included the muscle architecture and craniomandibular shape, combined with in vivo measurements of bite force from individuals of one breed (Malinois). Our results demonstrated excellent correspondence for dogs of similar size and shape, thus validating our model. Similar validation for other breeds is needed, however, to be able to confirm that our model provides reliable results irrespective of breed or shape.

The m. temporalis was found to contribute to half of the estimated bite force. Its contribution even tended to increase when the gape angle increased, which is related to the significantly increasing contribution of its deep bundle, and shows that it provides a performance advantage at large gapes. This is consistent with the need for carnivores to produce high bite forces at large gapes (Herrel et al., 2008; Turnbull, 1970). As demonstrated for other species (Herrel et al., 1998b; Bourke et al., 2008; Cleuren et al., 1995; Dumont and Herrel, 2003; Herrel et al., 2008; Kerr et al., 2017), the bite force and angle of reaction force in the joint (relative to the upper jaw) decrease as the mouth opening angle increases, and as food reaction forces move away from the orthogonal to the lower jaw. The reaction force, in contrast, increases.

Our results show that bite force is extremely variable in dogs and that it increases as size increases, as expected. However, we found significant differences in the residual bite force between brachycephalic dogs and the two other morphotypes. These results are consistent with the results of Ellis et al. (2009) and are coherent with lever mechanics. A short out-lever transmits a high force for a small movement: breeds with a shorter lower jaw will produce a relatively higher bite force (Slater et al., 2009). Conversely, a long lever arm (related to a longer lower jaw, as in dolichocephalic dogs) transmits a lower force but the amplitude of its movement is larger. We found no significant difference between the mesocephalic and dolichocephalic dogs, but this is probably due to our sample which contained very few extreme dolichocephalic dogs. The limits between the groups are further somewhat arbitrary and depend on the definitions used by different authors (Ellis et al., 2009; Miller et al., 1965; Roberts et al., 2010). To make homogeneous groups, we chose a relatively high value of the cephalic index of 0.6 to distinguish mesocephalic and dolichocephalic dogs. However, Miller et al. (1965) stated that the mean value was 0.39 for dolichocephalic dogs versus 0.52 for mesocephalic dogs and 0.81 for brachycephalic dogs. This could explain why we did not find a difference between these two groups. It would be interesting to add some graioid dogs (greyhounds) and other large dolichocephalic dogs (Leonberger) in future analyses.

The huge diversity in bite force observed can be expected to be related to the use of the dogs and the task they have been bred for. The distinction between the large brachycephalic dogs and the other dogs (small brachycephalic, mesocephalic and dolichocephalic dogs) in our sample might indeed be the result of different selection practices. Most of the large brachycephalic dogs are historically dedicated to the protection of humans (such as the rottweiler), whereas small brachycephalic dogs are dedicated to companionship, and mesocephalic or dolichocephalic dogs are dedicated to herding or hunting. For skills such as protection or attacking, breeders try to improve biting or gripping ability. Thus, it is not surprising to observe relatively higher bite forces in large brachycephalic dogs that were bred for defence/attack rather than dogs bred for herding or hunting, which are more commonly dolichocephalic. Previously, trade-offs between these two functions have been demonstrated in the musculoskeletal system (Cameron et al., 2013; Helton, 2011), suggesting that running/scent hound dogs are likely to be poor biters. Interestingly, the Leonberger, despite being dolichocephalic, is classified in group 2 of the FCI, with molossoid breeds and mountain-type dogs such as the rottweiler and cane corso. The breed is derived from a mixture of a Newfoundland, a grand St Bernhard and a Pyrenean mountain dog. It is thus not surprising that this dog produces bite forces as elevated as those of the other large brachycephalic dogs in our sample. To investigate the influence of inbreeding and genetic heritage on the functional abilities of the masticatory apparatus, a much bigger sample, including at the intra-breed level is, however, required.

A high intra-breed variability was also observed. If we consider the 10 beagles we dissected, calculated bite forces ranged from 262 to 466 N on the incisor teeth (mean: 359 N), 301 to 543 N on the canine tooth (mean: 418 N) and 559 to 1018 N on the carnassial tooth (mean: 790 N). Age, size and sex are probably important drivers of this variability. Moreover, changes throughout life may also influence bite force, as mammals are very plastic even at late life-history stages (Scott et al., 2014). Differences related to pathologies may further influence cortical bone modelling or muscle architecture. Diet and training probably also influence muscle development and bone shape and need to be taken into account to understand the intra-breed variability in bite force. The influence of training is a fundamental issue that would be worth exploring further. Indeed, exercise can improve masticatory function (Bourke et al., 2008; He et al., 2013; Kiliaridis et al., 1995; Kim et al., 2018; Lindner et al., 1995; Shirai et al., 2018; Thompson et al., 2001) but few studies have investigated this in detail.

The shape of the upper and lower jaws is significantly related to the absolute or residual bite force. The 2B-PLS analysis indicated that the curvature of the mandibular body, the relative size of the coronoid process and its processes, as well as the shape of the zygomatic arches relative to that of the neurocranium (which determines the space that is available for the adductor muscles to pass through) influence bite force variation. Overall, the presence of hypertypes with malocclusion between the lower and upper jaws (frequent in the small brachycephalic dogs), did not seem to alter the patterns of covariation much. However, there was significant variability in the data and caution is needed when interpreting the functionality of specific shapes, in particular with regard to the small brachycephalic dogs. In the 2B-PLS analysis with absolute bite force (Fig. S2), these dogs slightly diverged from the overall pattern. The shape of the upper or lower jaw alone can therefore lead to an overestimation of the absolute bite force in these small dogs, hence the need to take into account size in the estimates. The visualizations provided by the 2B-PLS analysis with residual bite force and the shape of the upper jaw show that similar shapes along axis 1 can produce very different relative bite forces, especially in brachycephalic dogs (e.g. the papillon dog versus the cane corso dog; Fig. 6): the small brachycephalic dogs produce much lower relative bite forces compared with other dogs of similar shape along axis 1. This supports previous observations by Ellis et al. (2009). The authors suggest that there is an interaction between shape and size for bite force and that shape may not be a significant factor in determining bite force in small brachycephalic dogs. Further studies including more hypertypes (small brachycephalic dogs) would be necessary to confirm that domestication did not completely disrupt the patterns of integration. Finally, it would be interesting to compare the coefficients with those of wild or commensal canids (wolves, dingoes) to test whether domestication has led to a decrease in the functional integration in the masticatory apparatus.

The shape of the lower jaw appears to be a better predictor of absolute bite force than the shape of the upper jaw (i.e. the results for the upper jaw are not significant for the Procrustes ANOVA; the variability of the point cloud along axis 1 of the 2B-PLS of the upper jaw is greater). This is consistent with the specialization of the lower jaw towards a single function (mastication), while the upper jaw has to cope with many functional demands related to the sensory organs, protection of the brain, etc. The strong relationship observed between the lower jaw and bite force should enable us to make predictions of bite force based on bone remains from the fossil record. Our study focused on the relationships between the overall shape of the bones and bite force, but did not allow us to explore the relationships between bone structure and the loads imposed during mastication. Finite element analysis may be an interesting complementary approach for this purpose (Bourke et al., 2008; Kim et al., 2018; Penrose et al., 2020; Wroe et al., 2007). Moreover, exploring the link between bone cortical thickness and bite force may be of interest to track functional variation according to load resistance (Cox et al., 2015; Kupczik et al., 2007; Rayfield, 2007; Ross et al., 2005).

The use of data on the muscle architecture obtained from dissections enabled us to describe the functional links between the muscular and bony components of the jaw system. The extreme variability in bite force in dogs is related to the extreme variability in size and shape, with brachycephalism conferring a mechanical advantage, as well as a great variation in muscle architecture (PCSA). Overall, it seems that the masticatory system is strongly integrated in dogs and that the strong relationships between the lower jaw and bite force observed are promising in terms of predictions using the 3D shape of the mandible only (which may be interesting, for example, in archaeology).

We thank Manuel Comte, Mickaël Godet and Frederic Lebatard for their help in managing specimens and their helpful discussions about the preparation of the skulls. We also thank Marie-France Varlet, member of the Société Centrale Canine and President of the Belgian Shepherd Dog Club (France). 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.