Unexpectedly rapid evolution of mandibular shape in hominins

Members of the hominins – namely the so-called ‘australopiths’ and the species of the genus Homo – are known to possess short and deep mandibles and relatively small incisors and canines. It is commonly assumed that this suite of traits evolved in early members of the clade in response to changing environmental conditions and increased consumption of though food items. With the emergence of Homo, the functional meaning of mandible shape variation is thought to have been weakened by technological advancements and (later) by the control over fire. In contrast to this expectation, we found that mandible shape evolution in hominins is exceptionally rapid as compared to any other primate clade, and that the direction and rate of shape change (from the ape ancestor) are no different between the australopiths and Homo. We deem several factors including the loss of honing complex, canine reduction, and the acquisition of different diets may have concurred in producing such surprisingly high evolutionary rates. This study reveals the evolution of mandibular shape in hominins has strong morpho-functional and ecological significance attached.


Evolutionary rate analysis results
With the FULL dataset a significant rate shift applies to the clade including the genus Homo and the Australopiths. Potential shifts are estimated according to variation in the multivariate Brownian rate, yet the same results apply by analyzing RRphylo rates ( fig. S2). We estimated whether the shift located at the clade parental to humans and the australopiths, as applied to RRphylo rates, is significant as well. To this aim, we first calculated the multivariate rate by taking the Euclidean Norm of the vector of individual variables rates (variables are the relative warp scores) per branch.
Then we took the average difference between the absolute values of the multivariate RRphylo rates of the hominin clade and rates attached to the branches of the rest of the tree. Then, we assessed significance by computing 10,000 average differences obtained by randomizing rates across the tree branches. With the FULL dataset, this test indicates that the branches belonging to the hominin clade are significantly higher than background rates at p <0.0001. The average multivariate rate difference is 0.065. With the SMALL dataset the average difference in the multivariate rates between the hominin clade and the rest of the tree is again significant at p <0.0001, yet the computed average multivariate rate difference is slightly smaller than with the FULL dataset at 0.025.

Accounting for phylogenetic uncertainty in node age and topology
The distribution of evolutionary rates depends on the distribution of branch lengths and on the tree topology (Bapst 2014). Every phylogenetic tree represents at best a phylogenetic hypothesis, which should be evaluated against alternative topologies, and branch lengths. To account for phylogenetic uncertainty, we wrote a Rcode that changes the tree topology and branch lengths. For every given species, the function swaps the phylogenetic position up to two nodes distance. For instance, the topology ( (A,(B,C)),D) could be swapped to the forms ((C,D),(A,B)); (((B,D),A),C) and so on. In addition, each node age is randomly set at any age between the age of its parental node, and the age of its oldest daughter node. We applied the tree swapping function 100 times, computed RRphylo rates at each time, and draw the difference in mean absolute rates between the hominin clade and the rest of the tree each time. The distribution of mean absolute rate differences ( fig. S3) still points to higher absolute rates in the hominin clade branches, and it is statistically significant, as assessed by means of a t-test (t = -10.227, df = 11.482, p << 0.0001).

Geometric Morphometrics of Primate mandibles
We used Geometric Mophometrics (Gmm, Rohlf & Marcus, 1993;Klingenberg, 2010) to extract morphological data. This method permits to retrieve shape information of anatomical objects after removing non-shape variation (i.e. as related to size, position and orientation of the objects) by applying Generalized Procrustes Superimposition (GPA, Rohlf and Slice 1990). By using the TpsRelw software ver. 1.53 (Rohlf, 2013b) we performed Relative Warps Analysis on aligned coordinates (RWA, Bookstein 1991;Rohlf 1993;Zelditch et al. 2002) to decompose shape variation into orthogonal axes of maximum variance.
For this study we collected (either by taking pictures directly, from digital sources, or from We used tpsDig2 software (Rohlf, 2013) to digitize 9 landmarks as to adequately describe the lower jaw profile ( fig. S4). Gmm also returns the Centroid Size, an index that permits to get back the information related to size that are removed by GPA. We regressed the natural logarithm of centroid size (lncs) and body mass estimates taken from the literature, to asses whether lncs works good as a proxy for body size. The regression is highly significant and positive (slope = 0.300, R 2 = 0.844, p < 0.001, fig. S5).
Shape variance was decomposed into 14 axes (Relative Warps). We performed the Gmm analyses twice: on the full dataset, and on a dataset deprived from pictures we obtained from literature. The former dataset (FULL) consists of 211 species, the reduced dataset (SMALL) includes pictures for 158 species (145 extant, 13 extinct).

fig. S5 -The regression between the centroid size (y) and body mass (x).
With the FULL dataset, RW1 and RW2 axes describe 38.59% and 19.50% of shape variance, respectively. Positive values of RW1 are associated to jaws with slender mandibular corpus and elongated surface of attachment for the masseter. Negative values of RW2 are associated to mandibles with stronger angle and lower ramus ( Fig. 1 Figure 2).
On the FULL dataset we performed a second GPA eliminating allometry, in order to account for potentially large scaling effects linked to the wide body size range in the data (Zelditch et al. 2012).
To this aim, we regressed the aligned coordinates versus lncs and took the residuals to use to perform GPA. With the 'size free' dataset, the first two RWs explain 32.9% and 24.9% of the shape variance, respectively (table S2, fig. S6, S7). We additionally repeated GPA in the 'size and shape' space, and obtained 15 new axes, in particular RW1 and RW2 explain 94.35% and 2.3% of shape variance respectively ( fig. S8).
Eventually, on the SMALL dataset we obtained 14 RWs. RW1 explains 44.35% and RW2 18.88% of shape variance. Along these axes the deformations are qualitatively the same that we obtained from first GPA with 211 species (fig. S9).