Materials, landmarks and geometric morphometrics

The elbow-joint shape in the 139 specimens (92 living specimens distributed among the families Felidae, Hyaenidae, Viverridae, Mustelidae and Canidae of the order Carnivora (Table 1) and 47 fossil specimens of the family Canidae including representatives of the three subfamilies, Hesperocyoninae†, Borophaginae† and Caninae17,18,19 (Table 3) was recovered by digitizing six landmarks in 2D on the anterior surface of the humerus distal epiphysis20,21,22 from high-resolution digital images (Fig. 1a). The digital images were taken with a tripod and following a standardized protocol for avoiding lens distortion and parallax.

We sampled those families of living carnivores with large representatives such as canids, felids and hyenids, because the inclusion of viverrids and mustelids do not add crucial information in terms of body size or predatory behaviour. However, we also sampled the wolverine because it is the largest mustelid (Gulo gulo), an ambush predator, and the large African (Civettictis civetta) and Indian (Viverra zibetha) civets.

The six landmarks were digitized into 2D Cartesian coordinates (x,y) using TPSdig V. 211 (ref. 38; Supplementary Data 1) but we modelled the inter-landmark distances by means of an outline to obtain clearer shape transformation models in subsequent multivariate analysis. In addition, we superimposed the 2D outlines into a three-dimensional scanned surface of a jaguar, Panthera onca (139959, AMNH) specimen to obtain clear interpretations of the shape transformation recovered in the axes derived from different multivariate approaches. All the specimens digitized were aligned using Procrustes and projected onto the tangent space39. The size of the specimens was represented by their centroid size (Cs), which allows post-hoc assessment for allometric effects40. Centroid size is the square root of the sum of squared distances of each landmark from the centroid of the configuration40. The calculation of Procrustes Coordinates (Proc) and Cs were performed with MorphoJ41.

We used the same six landmarks as Andersson20, as they are established morphological indicators of forearm pronation and supination. However, whereas Andersson20 analysed solely the Euclidean distances among landmarks, we used an approach based on Geometric Morphometrics. Our approach not only captures the relative distances among the landmarks, but also the topological information contained within them.

The influence of size and phylogeny in elbow-joint shape

We assembled a phylogenetic consensus tree (Supplementary Fig. 1) following various published sources using Mesquite42. We incorporated branch lengths in our composite phylogeny in million years before present43,44. In the case of living species, branches were scaled according to node dates estimates based on Nyakatura and Bininda–Emonds45. In the case of extinct taxa, fossil occurrence dates were compiled from various sources based on species locality and age information (Table 2), and the branch lengths were estimated from the first and last appearance of taxa (Supplementary Fig. 1). We used this phylogenetic tree to assess phylogenetic patterning in our data.

As we are only interested in quantifying the presence of phylogenetic signal, instead of the strength of the phylogenetic signal, we used the permutation approach developed by Laurin46, extended for multivariate analysis by Klingenberg and Gidaszewski47, and applied to shape data by other authors (for example, refs 48, 49, 50, 51, 52, 53, 54), to simulate the null hypothesis of complete absence of phylogenetic signal in elbow shape. The mean species shapes are randomly distributed as the tips of the phylogeny in 10,000 permutations, and for each permutation, the tree length was computed. If the resulting tree length computed for each permutation was greater than the one obtained with the original data, the null hypothesis of absence of phylogenetic structure was rejected. A P-value was used for assessing the presence of phylogenetic signal in shape47.

A multivariate regression analysis55 of shape (that is, using Proc) on size (that is, using Cs) was performed to test the influence of allometry. In addition, we applied independent contrasts analysis56 to take into account the phylogenetic relationships of the species under study. Following this, we regressed the contrast for shape (Proc) on the contrast for size (Cs) using MorphoJ41. The statistical significance of both multivariate regression analyses was tested with a permutation test against the null hypothesis of complete independence of shape on size57.

We performed a multivariate regression analysis55 of the independent contrast of the Proc against the s.d. of the standardized contrast (that is, the square root of the corrected branch lengths) with MorphoJ41, following Díaz-Uriarte and Garland25. The values of the s.d. were obtained from the PDAP module for Mesquite42,58. The significance of this regression was evaluated with a permutation test against the null hypothesis of complete independence between the two variables. This test was specifically performed to explore the adequacy of: (i) the model used for tree topology; (ii) the branch lengths used; and (iii) the model of Brownian motion for our tip data25.

Inferring predatory behaviour in extinct canids

To determine those elbow-joint shape features that best distinguish among the three present-day predatory modes (ambush; pounce-pursuit; and pursuit), we performed a CVA from the Proc of elbow shape in modern predators. We used CVA instead of Principal Components Analysis (PCA) as performed by Andersson20 because to find those elbow-shape features that best distinguish among these predatory groups of specimens, the use of CVA is more appropriate than PCA. In fact, while CVA is a classification method, PCA is an ordination one. In any case, it is worth mentioning that although Andersson21 also used CVA, the analysis was performed to differentiate among the three subfamilies of canids instead to infer their predatory behaviour.

We classified all of the living species into one of the three present-day predation modes following previously published sources23,59 as follows: (i) ambush predators stalk their prey and may pursue them over short distances, and the forelimbs may be used to grapple with large prey; (ii) pounce/pursuit predators usually hunt small prey using either a pounce or short chase, and rarely grapple with their prey; And (iii) pursuit predators usually chase their prey for a long distance (>30 m), and may hunt cooperatively to bring down large prey, but do not grapple with their prey (Table 1). Therefore, we consider here the pounce-pursuit category of Van Valkenburgh23 instead of only ‘cursorial’ and ‘non-cursorial’ predators as Andersson20 differenced. However, these three hunting types are usually correlated with prey size, while ambushers can take prey of all sizes, pounce predators usually take small prey and pursuit predators usually take large prey. Both ambush and pursuit predation can be directly linked to forearm mobility and elbow function in the case of large prey, but a predator will not ‘pounce’ on large prey nor ‘pursue’ small prey. We consider the cheetah (Acinonyx jubatus), the only highly cursorial felid, to be a pursuit predator, even though this species still retains some ability to supinate as they often swat their prey (but they do not grasp it). The reason is because the cheetah is a solitary hunter but does not have the endurance of the canid and hyenid pursuit predators and it will bring down its prey with the swipe of a forelimb. However, the cheetah is a felid that is clearly more specialized towards running than others in its family, pursuing its prey for a distance of up to 200 m and it has been classified as a pursuit predator by other researchers15,59.

As the predatory categories used here are dynamic and not static, our predatory grouping scheme established by Van Valkenburgh23 represents an arbitrary division of a continuous range. Furthermore, many carnivorous taxa range over some predatory types (for example, the cheetah), and so the assignment of a given species to a particular category is a ‘best fit’ designation rather than an exclusive one. In summary, our predatory categories are thus a significant (albeit necessary) simplification of the complex range of predatory behaviour.

The reliability of group separation was assessed by the pairwise Procrustes distances and Mahalanobis distances among all possible pairs of groups using the pooled within-group covariance matrix for all the groups jointly. A permutation test was computed for the Mahalanobis and Procrustes distances of all pairwise comparisons. The statistical significance of these pairwise differences was assessed with 10,000 permutations using MorphoJ41. Both canonical functions obtained from the living sample in CVA were later applied to fossil taxa using their Proc of elbow shape. However, although we will never be sure if the established predatory categories based on living carnivores can be extrapolated to extinct species, this uncertainty is alleviated if fossil forms have extant relatives. This is particularly the case of the family Canidae, as extant forms phylogenetically bracket them, where the morphological correlates of behaviour are known.

We used the direct method of leave-one-out cross-validation for assessing the percentage of probability of living canids to belong for one of the three present-day predation modes. The classification of fossil species into one of the three predatory groups was inferred according to their proximity to group centroids of the three predatory groups of extant carnivorans with SPSS v.19.

We did not use the stepwise discriminant analysis60 because we have more cases per group than predictor variables; 12 variables or six bi-dimensional landmarks digitized in 35 pounce-pursuit predators, 30 ambush predators and 22 pursuit predators (Table 1). Thus, there is no statistical reason to perform a stepwise approach instead of the direct method because CVA only tends to over fit differences in those cases where there are more variables than the number of cases within groups61.

The morphospace depicted from the scores of the living specimens into both canonical axes was imported into the graphing software SigmaPlot to create a contour plot that was colour coded by predatory mode—extant ambushers in blue, extant pounce-pursuit predators in green and extant pursuit predators in red. Colour degradation within each predatory group was calibrated to the frequency of the specimens within each predatory group. After we performed this contour plot from the scores of the living taxa, we repeated the same graphs but now including the fossils according to their predicted predatory category obtained in CVA. We recovered a RGB colour code for each fossil specimen according to its position in the contour plot living morphospace.

To test whether carnivores with different predatory behaviour differed in shape irrespective of phylogeny, we used the aov.phylo function in the R package ‘geiger’62 using both canonical axes. We used Brownian motion as a model for evolutionary change and ran 1,000 simulations to create an empirical null distribution to compare with our sample.

Ancestral elbow-shape reconstruction

All taxa excluding canids were pruned from the phylogenetic tree shown in Supplementary Fig. S1 using Mesquite42. The ancestral states, or internal phylogenetic nodes, were reconstructed by squared-change parsimony26 weighted by branch lengths using MorphoJ41. The main purpose to reconstruct ancestral states was to explore patterns of elbow-joint shape evolution and by extension the evolution of predatory behaviour in canids. Also, we would be able to assess the most probable basal hunting technique for the three subfamilies and when the specialized traits of the elbow are acquired through the evolution of canids.

Body mass estimation in extinct taxa

To explore if elbow-shape changes experienced by the evolution of the family could be in part explained by a previously reported trend towards size increase through canid evolution37, we explored the issue of whether the appearance of canids with morphology indicative of more cursorial behaviour (pounce-pursuit or pursuit) is the result of their body mass reaching or exceeding 21.5 kg, which is a threshold value in carnivore biology: only carnivores of this size or greater take prey as large as, or larger than, themselves63.

We estimated the body masses (BMs) of fossil canids by means of a least-square bivariate regression of log-transformed BM on Cs for the extant taxa (log (BM)=a+bLog[Cs]) using SPSS v.19. Although both body mass and elbow centroid size are variables, each with its own errors and distributions under the control of the investigator, we chose the least-square regression model (instead of a type II) because our purpose here was to explore the body mass dependency of elbow centroid size. Note that the main objective here was to obtain a model in which body mass can be predicted from elbow centroid size. For these reasons, we assumed dependence between both variables64.

The body mass for each extant species was taken from recent studies on body mass estimation in extinct taxa64,66 and references therein. Their Cs were then computed as the Cs average for all the specimens belonging to the same species. The accuracy of the bivariate regression function was evaluated by means of the F-statistic and Pearson’s correlation coefficient (r). However, as the correlation coefficient can be high even with high residuals67, it is not a highly appropriate statistic to evaluate the predictive power of a regression equation. Thus, to evaluate the predictive power of the function, we computed the per cent prediction error (%PE) and the per cent standard error of the estimate (%SEE)67. Once we had confirmed the significance of the Log-BM on Log-Cs regression equation and its predictive power, we predicted body mass values for the extinct canids.

Enamel microstructure specialization in extinct canids

The association between environmental change and canid predatory behaviour was investigated using the arrangement of the enamel prism bands—known as bands of Hunter–Schreger (HSB)—as a proxy. The reason is because carnivores that live in open environments consume more intrinsic (bones from carcasses) and extrinsic (grit) hard food items, which relates to heavily folded HSB28,51,68,69,70. Studies of carnivorans demonstrate a link between bone and carcass consumption and increased kleptoparasitism in open habitats because of increased carcass detectability30,31,32. Therefore, we utilize HSB analysis as proxy of habitat openness, in terms of increased enamel resistance to wear from ingested grit, as well as hard vertebrate tissue consumption associated with increased carcass availability in open environments.

To quantify the degree of specialization in the enamel prism bands arrangement (bands of Hunter–Schreger; HSB), we used a stereomicroscope at 5–40 × magnification. We coded for each tooth in the dentition, the degree of HSB specialization in each of the three regions on the tooth crown (base, tip, and the mid region), each occupying approximately one-third of total crown height. We quantified the degree of HSB specialization as follows: the less-specialized undulating pattern; then the intermediate acute-angled undulating pattern; and followed by the most specialized zig-zag HSB pattern70. The final data set included both newly collected data as well as published data in Tseng69 (Supplementary Table 5).

The coded HSB pattern for each tooth was given as a score from 0 (for undulating pattern in all three regions of the tooth crown) to 1 (for zig-zag pattern in all three regions of the crown). We studied lower cheek dentitions (premolars and molars) of fossil canid species. However, we used the upper dentitions for those species without preserved dentaries. Furthermore, where both left and right dentitions were available, averages of scores were taken across the corresponding tooth and treated as a single entry in the data analysis. Because not all fossil specimens contained complete data for the premolar–molar sequence, 10% of the data set represent composite HSB data from multiple specimens of the same species. The final HSB score for each specimen data point was calculated per tooth (total HSB score divided by number of teeth coded per specimen) to account for differences in completeness and tooth count of sampled species (Supplementary Table 5).