Data Acquisition

To gather three-dimensional (3D) morphological data for the taxa studied via FEA, specimens were CT-scanned via high-resolution micro-computed tomography (μCT) using the GE v|tome|x s 180/240 kV dual-tube HRX μCT system at the Microscopy and Imaging Facility (MIF) of the AMNH. Scans were made at a voltage of 150–170 kV, tube current of 55–180 mA, and voxel size of 37.69–136.00 μm, depending on skull size (see Table S2 for additional model details). Other species models were derived from previous FEA studies. All specimens analyzed as a part of this study are housed in museum collections in South Africa, Spain, and the U.S.A. All leptarctine specimens in this study are from the AMNH Division of Paleontology Fossil Mammal collections, including Leptarctus primus AMNH FM 25385 (72 microns voxel size) and Leptarctus sp. AMNH FM 54198 (76 microns voxel size), from the Frick Collection; and Hypsoparia bozemanensis (cast) AMNH FM 48829 (84 microns voxel size), L. primus (neotype) AMNH FM 18241 (cranium; 52 microns voxel size), and AMNH FM 18270 (mandible, not scanned) from the Fossil Mammal (FM) collection. The Leptarctus model used in subsequent FEAs was based on AMNH FM 25385; AMNH FM 18241 was used as a visual reference but was not combined in any way with AMNH FM 25385. Canines from AMNH FM 54198 were transposed digitally onto AMNH FM 25385 to restore that part of the skull. The cast specimen AMNH FM 48829 was digitally grafted onto AMNH FM 25385, replacing the zygomae to represent a hypothetical, Hypsoparia-like proxy model for additional analysis. Newly sampled extant non-leptarctine carnivoran specimens, including Bassariscus astutus (AMNH 135964; 56 microns voxel size), Gulo gulo (AMNH 182936; 74 microns voxel size), Potos flavus (AMNH 239990; 72 microns voxel size), Spilogale putorius (AMNH 35207; 58 microns voxel size), Taxidea taxus (AMNH 120577; 82 microns voxel size), and Urocyon cinereoargenteus (AMNH 147213; 51 microns voxel size), are from the AMNH Division of Mammalogy collections. Non-leptarctine taxa (all extant, unless noted) FEA results from previous studies include Canis lupus (Tseng, 2009); Canis mesomelas (Slater et al., 2009); Ursus maritimus and Ursus arctos (Slater et al., 2010); Lycaon pictus and Crocuta crocuta (Tseng and Stynder, 2011); Parahyaena brunnea (Tseng, 2013); Ailurus fulgens (Figueirido et al., 2014); Panthera pardus, Herpestes javanicus, Mephitis mephitis, and Procyon lotor, as well as the extinct creodont Thinocyon velox and carnivoramorphan Oodectes herpestoides (Tseng and Flynn, 2015b).

As far as we are aware, there have been no prior quantitative tests of the biomechanical properties of leptarctines other than bite force estimates using the dry skull method, which estimates muscle physiological cross-section areas using 2D photos (Calede et al., 2017). Therefore, one goal of this study is to establish a 3D anatomical and biomechanical framework in which to quantitatively compare L. primus and H. bozemanensis with the diverse array of previously listed extant and extinct Carnivora taxa. We surveyed 18 extant specimens to ensure representation across carnivoran clades at different taxonomic levels, creating a broad phylogenetic context for examination of cranial biomechanics of Leptarctus. We categorized the trophic level of each extant species using a combination of classification provided in the PanTHERIA database (Jones et al., 2009) and published natural history accounts of the species in question. The trophic levels recognized in this study include carnivore, omnivore (including insectivores), herbivore, and frugivore (Table S1).

Finite Element Analysis

Finite element analysis in biology typically involves a simplified model of a structure of interest that is analyzed with simulations representing biological phenomena under a series of simplifying assumptions about the structure being modeled (Rayfield, 2007). Most importantly, by focusing on a comparative analysis rather than a derivation of absolute magnitudes, the FEA approach employed here simplifies the cranium in several ways. Cranial models incorporate load-bearing cortical bone structure only, and do not account for variations in the contact between bone and tooth row, in order to generate a taxonomically broad sample of models to analyze through phylogenetic comparative methods (see below).

By implementing the same approach as in previous studies for all cranial models in this analysis, we can directly compare simulation outputs across studies, thereby greatly expanding the sample sizes available for comparative analyses in a timeintensive but informative methodology that incorporates 3D osteological information. Following the protocols in Tseng and Flynn (2015a, 2015b, 2015c), developing a finite element model includes the following major steps: (1) shape reconstruction of the object of interest (in this instance, various cranial structures); (2) specification of the material properties of that object, or rather the characteristics of the cranial structures that govern the way in which the object moves under mechanical stresses and strains (such as a simulated bite); and (3) a simulation in which the object or structure is loaded with mechanical forces (the stresses and strains on the skull of a bite), otherwise understood as a parameterization of the model. This simplification of form is universal among tested crania, especially in light of varying fossil quality and diagenetic deformation, which render ex vivo experimental biomechanical assessments of fossilized specimens extremely challenging or impossible. The effect of the model simplifications made in our protocol permit comparison of relative, but not absolute, magnitudes of the biomechanical parameters. For sensitivity of model outputs to the various parameters employed here, see Bright and Rayfield (2011b), Tseng et al. (2011), Tseng and Flynn (2015a), and Walmsley et al. (2013).

Coronal CT scan images of each specimen were imported as an image sequence in ImageJ (Schneider et al., 2012) and standardized for data set size by subsampling the image stacks to represent similar numbers of slices and pixel resolution. This was achieved using the ‘import sequence’ option and setting the sampling interval to reduce the number of final images used. The coronal images were then exported from ImageJ into Mimics (Materialise NV) where segmentation of the regions of interest was conducted. In Mimics, threshold and manual edit operations were conducted to select only cortical bones and to separate the cranium from the mandible. The segmented regions were then reconstructed into 3D surface meshes and exported from Mimics into Geomagic Studio 12 (3D Systems) as STL files, where the cranium surface mesh underwent improvement operations to create ‘water-tight’ structures for solid meshing. Using Geomagic Studio 12, meshes were cleaned of overlapping and skewed triangular elements both internally and topically. For damaged or incomplete skulls, modest reconstruction took place to restore symmetry or patch gaps in the preservation. Because none of the skulls in our study had major destruction, it was possible to highlight entire sections of the skull, mirror a section, and attach the newly mirrored section to the damaged area, thus eliminating the potential errors or conjecture involved in hand-crafting a new segment of the skull. Initial appraisal of the models was undertaken by opening the ‘display’ tab, allowing the model transparency, and ‘clipping’ the model, or creating an active and moveable cross-section of the skull to identify collapsed elements and holes. The model was then smoothed using the ‘quicksmooth’ function to flatten the slice artifacts produced by the coronal image stack compilation; problematic holes, tunnels, and spikes (sharp edges) within the mesh were repaired by the ‘Mesh doctor’ function, which identifies and corrects common meshing errors. Special attention was paid to the tooth cusps at this stage in order to avoid unintentional altering of tooth cusp height or morphology. The triangular surface elements that compose the meshes were reduced via the ‘decimation’ function to a more uniform number (between 200,000 and 280,000 surface triangular elements) among all specimens. Because finite element densities have been known to affect simulation outputs, this was a necessary standardization step (Tseng and Flynn 2015a, 2015b, 2015c). At the time of triangular element ‘decimation,’ triangles were constrained to have a controlled aspect ratio (edge/edge is 8, edge/height is 6) to remove highly skewed elements within the model.

The refined reconstructions were then imported into STRAND7 (G+D Computing) as STL files. Skulls were checked for extraneous nodes, erroneous triangles in the meshes, and other errors that would interfere with solid meshing. Iterative mesh improvements were conducted in Geomagic after using the error-checking functions in STRAND7 if collapsed elements or overlapping edges were found, which occurred in roughly 80% of the meshes. Once all errors had been resolved within the surface mesh, FE meshing and model boundary conditions could be initiated. The surface was then converted into a solid mesh using the automatic meshing function in STRAND7; meshes of three different resolutions were generated for each model to account for variability of simulation outcomes as a function of mesh density (Tseng and Flynn, 2015c). The muscular attachments of the cranium were identified for the temporalis, masseter, and pterygoid groups in Geomagic, then exported to MATLAB. The BONELOAD script (Dumont et al., 2007) was used to apply vector forces evenly over the surface of the muscle attachment sites using the ‘tangential distribution’ function; all muscle forces were calculated by multiplying muscle attachment area by 0.3 N, the maximum contractile force for mammalian muscle fibers (Tseng and Flynn, 2015a, 2015c). Because only unilateral bites were simulated, muscle forces were adjusted depending on whether the muscle group in question is on the working side (100% of calculated force) or balancing side (60% of calculated force) (Table S2). Constraints to the models were applied at each of the temporomandibular joints (TMJs), and at a given tooth position (at the tip of the tallest cusp) using single nodal constraints. The TMJ constraints allowed translation only along the axis formed by the TMJs, and bite position constraints allowed translation in the plane perpendicular to the vertical axis of the tooth. All cheek dentition positions present (from canine to last molar) were simulated. All models were given the same set of homogeneous material properties of 20 GPa (Young's modulus) and 0.3 (Poisson ratio). All models were solved using the linear static solver in STRAND7.

Raw data from FEAs were standardized to produce the comparative FE data set of 21 species (Table S3). The magnitude of bite force (measured as nodal reaction force at the bite position) was scaled relative to input muscle force (the sum of all calculated muscle forces on a given model) to produce a measure of mechanical efficiency (ME) (i.e., more efficient cranial systems can produce more bite force per muscle force input). Because ME is a ratio, it is a size-free measure of biomechanical efficiency. Skull strain energy (SE) measures skull stiffness (as work put into structural deformation by the muscle forces), which is biomechanically significant, because lower strain energy values indicate more input muscle force into the output bite force as work done on food during mastication (rather than work done to deform the cranium). Therefore, low SE and high ME values would indicate a lower risk of damage to the cranium and higher bite force relative to amount of muscle force input and therefore high maximum capability to utilize harder food sources (Tseng and Flynn, 2015b).

In addition to comparing numerical FEA outputs, we examined von Mises stress generated by bite simulations for the species deemed most similar to Leptarctus in multivariate analyses (see below) to determine similarity in regional skull stresses and strains during bites in each tooth position (canine, second premolar, third premolar, fourth premolar, and first molar). These heat maps indicate the weakest areas of the individual's skull when a particular tooth locus was simulated to ‘bite.’ Within an engineering context, von Mises stress approximates locations of likely failure in a ductile material, such as bone (Irons and Ahmad, 1980). In the context of skull deformation, von Mises stress indicates which areas of the skull are weakest under load: cooler colors (blues and greens) indicate less stress in that area, and warmer colors (reds and yellows) indicate greater stress. Areas that are white exceed the threshold of the index scale (0–50 MPa). Skulls with similar color codes (i.e., similar von Mises stresses) indicate comparable skull deformation during a ‘bite,’ which in turn suggests similar feeding capabilities.

Hierarchical Cluster Analysis with Resampling and k-means Analysis

The first set of multivariate analyses employed a series of clustering methods for classifying the data gathered through FEA (i.e., values of cranial strain energy [SE] and mechanical efficiency [ME] for each of the tooth loci when simulating a masticatory motion in STRAND7), with statistical analyses performed in R version 3.2.5 (GNU GPL version 2). The data were first run through hierarchical cluster analyses using the ‘hclust’ function in the STATS package (R Core Team, 2017). Eight different cluster methods available in the R STATS package (average, centroid, complete, mcquitty, median, single, Ward D, Ward D2) were employed using Euclidean distances between taxa's FEA output. The goal of employing cluster analysis was to identify the most similar taxon (or taxa) to Leptarctus in all of the biomechanical simulation values compared.

Bootstrapping is frequently implemented in quantitative paleontological analyses where data are pooled together and resampled via replacement (Kowaleski and Novack-Gottshall, 2010). The probability of observing the clusters generated using hierarchical cluster methods was therefore assessed using bootstrap resampling, with two types of bootstrap metrics to assess support of clusters: regular BP (bootstrap probability) and approximately unbiased (AU) BP, which additionally resamples the data at different sample sizes (Shimodaira, 2004). These analyses were performed on all cluster methods at 1000 iterations each. The bootstrap percentages generated by the resampling methods assess the significance of the percentage of times a specific cluster appears out of those one thousand iterations. Clusters with AU P-values greater than 95% were considered strongly supported by the data. The masticatory capabilities of Leptarctus primus then were assessed by the similarity in its strain energy and mechanical efficiency values relative to extant species with known dietary preferences.

In addition to hierarchical cluster analyses, k-means analyses were applied as another measure of cluster support (Tseng et al., 2016). Given that the data in each of our clusters are approximately Gaussian, it is logical to use the k-means algorithm to order numeric values (Huang, 1998). Given a cluster number (i.e., an a priori grouping into which specimens can sort), k-means calculates the distributional means of all the clusters (i.e., a mean vector for comparing the groups into which specimens are assigned) for classification. These data in turn facilitate two operations: (1) determining clustering association of specimens, then (2) permitting aggregate-level analyses on a specific cluster. Classification of groups at different preassigned k-means values provides a simple way of evaluating similarity and another means of measuring support for a particular group membership. We examined the 21-species data set through 19 iterations of k-means clustering (from k = 2 to k = 20; k = 1 would form a single cluster, and k = 21 would form 21 singletaxon clusters, so both were excluded) to assess the consistency of the clusters containing L. primus.

Phylogenetic Signal Test and Phylogenetic Principal Components Analysis

Phylogenetic signal can be defined as the tendency for nonrandom distributions of measured characteristics to be derived solely from evolutionary relatedness (Felsenstein, 1985). To test for phylogenetic signal in this multivariate data set, we used (1) a multivariate version of Blomberg's K statistic (Blomberg and Garland, 2002; Adams, 2014) and (2) a Mantel test (Harmon and Glor, 2010). Blomberg's K is a quantitative measure of phylogenetic signal strength as a function of (1) the ratio between data variance and variance with consideration of phylogenetic covariance and (2) the ratio between expected variance under a Brownian motion model of evolution given the phylogeny in question (Blomberg et al., 2003). The Mantel test assesses whether two distance matrices are significantly correlated. The adjusted biomechanical attribute–inferred simulations (SE and ME) were tested both separately and together. Significant phylogenetic signal in any of the data partitions would indicate that the clustering patterns observed can be explained as either partially or totally derived from evolutionary interrelatedness, whereas lack of significant phylogenetic signal would suggest that other factors (such as feeding adaptations) play a more important role in generating the data distributions observed. Presence of phylogenetic signal does not preclude an initial covariance between a functional shift and clade origination, but once that function was established within that clade, it does not require additional special explanation as a functional adaptation. Blomberg's K was assessed using a modified version of the function provided in the R package GEOMORPH (Adams et al., 2017), and the Mantel test was assessed using the ‘mantel.test’ function in the R package APE (Paradis et al., 2004) using Euclidean distance matrices for FEA data and branch length distance matrices for each of the six branch length configurations.

In addition to testing for phylogenetic signal, we conducted multivariate ordination using phylogenetic principal components analysis (pPCA) to visualize the main axes of variation within the FEA output data using the PHYTOOLS package in R (Revell, 2012). For both phylogenetic signal tests and pPCA, a composite phylogenetic tree was constructed in Mesquite based on the topology recovered in Flynn et al. (2005); placements for mustelid taxa not analyzed in Flynn et al. (2005) were based on Koepfli et al. (2008) (Fig. 2). We tested six branch length configurations to assess their effect on the calculation of phylogenetic signal and phylogenetic PCA(Table S4). These configurations represent either different fossil occurrence dates or molecular divergence time estimates from the literature (see Supplemental Data). Because Leptarctus has not yet been phylogenetically placed using cladistic methods within the framework used in the current study, we conservatively place Leptarctus as an unresolved branch of Arctoidea, with Ailuridae and other Musteloidea (Fig. 2).

Linkage between Dietary Category and Biomechanical Simulation Output

To assess the degree to which biomechanical data from FEA are related to dietary categories among extant species studied, we used (1) comparisons of ME-SE profiles that have been previously shown to distinguish hypercarnivores from omnivores (Tseng and Flynn, 2015b) and (2) univariate and multivariate Kruskal-Wallis tests between FEA data and dietary categories. The ME-SE profiles were constructed from the numerical output data obtained in the FE simulations (see Finite Element Analysis, above), and the profiles of ME and SE values across the tooth row were examined for presence of multiple inflection points of low SE values along the tooth row (indicating hypercarnivory) or either a smoothly sloping profile or presence of a single inflection point (indicating omnivory). Kruskal-Wallis tests are the nonparametric equivalents of analyses of variance: ANOVA (univariate) and MANOVA (multivariate); we used the ‘kruskal. test’ function from the CORE STATS package in R for univariate tests and the ‘nonpartest’ function (with n = 10,000 iterations) from the npmv package for multivariate tests (Burchett et al., 2017). The dietary categories are as previously defined using a combination of the PanTHERIA database and individual natural history accounts (Table S1). Because of the small samples of herbivores (n = 1) and frugivores (n = 1) in the data set, only carnivores and omnivores were tested as categories. Assessments of similarity of Leptarctus FEMs to herbivore and frugivore data were done using the ME-SE profiles only, without Kruskal-Wallis tests.

Functional Implications of Morphological Variation between Leptarctus primus and Hypsoparia bozemanensis

As noted above, Hypsoparia bozemanensis is a close relative or possible conspecific of Leptarctus primus, but it is not represented by complete material. To determine the possible range of variation in biomechanical characteristics and feeding capabilities of these two organisms, we assessed the effect on FE output values of the range of possible cranial variation by including distinguishing H. bozemanensis characteristics cited by prior authors (unicuspid p4, canine groove, and deeper zygomae) in analyses of a proxy model of the L. primus specimen (Dorr, 1954; Webb, 1969; McKenna and Bell, 1997; Baskin, 1998; Qiu and Schmidt-Kittler, 1982; Lim and Miao, 2000). It quickly became apparent that there was only one viable characteristic that differentiated the two. The widely referenced depiction of the p4 of L. primus had been incorrectly drawn as having a bicuspid p4. Instead, upon direct inspection of the specimen (AMNH FM 18270), we discovered that its p4 instead was unicuspid, with only a small cuspule rather than a distinct second cusp. Because neither L. primus nor H. bozemanensis possesses a bicuspid p4, this character was not analyzed further. Second, the canine groove is a mandibular trait and cannot be tested in this study of overall cranial biomechanics (rather than dental morphology). Therefore, the only remaining potentially variable character to include in the model was the deep zygomae. In the absence of a complete skull (i.e., there is only one verifiable H. bozemanensis skull known, and it is mostly fragmented, making it a poor candidate for CT scanning and FEA of a full skull), a cast of H. bozemanensis (CMNH 9674) was scanned and the deeper zygoma in that specimen was size-standardized and grafted onto the L. primus model. The combined proxy model was then analyzed with the same protocol as for the original Leptarctus model. The goal of this secondary set of analyses was to determine whether or not the individual models for H. bozemanensis (proxy model reconstruction) and L. primus performed similarly in FEA, even with differences in their zygomae.

ME-SE Profiles from Finite Element Analyses

Unilateral bite forces, summarized as the force of the average of both the left and right tooth loci, with 95% confidence intervals (CIs) calculated from multiple-resolution models analyzed per species, are presented in Figure 3. The FEAs of unilateral bite forces yield biomechanical profiles of ME versus standardized SE values plotted over all tooth locus simulations. These range in ME from 9.3% (the canine in the fossil creodont Thinocyon velox) to 77% (the first molar in Taxidea taxus). Within the 0.028–0.23 J range for SE, Potos flavus had the highest overall SE (lowest skull stiffness) and Spilogale putorius the lowest SE (highest skull stiffness). Leptarctus exhibits some of the highest SE values, with its unilateral bite force being just below that of Lycaon pictus, comparable to those of the fossil stem-carnivoramorphan Oodectes herpestoides and the extant canid Canis mesomelas in its low stiffness. However, in the progression of ME and SE values from anterior to posterior tooth loci, the profile of Leptarctus is similar to those of Lycaon pictus, Taxidea taxus, and Thinocyon velox in having two inflections along the curve (but note that 95% CIs encompass the whole range of SE values across tooth loci) (Fig. 3; see also Table S3). Leptarctus primus has the second highest ME value, below only Taxidea taxus. Leptarctus primus and T. taxus also have more similar ME-SE bite profile shapes across the tooth row than to Po. flavus, which has a single, significantly lower SE inflection point relative to other tooth positions in its profile. Gulo gulo, although exhibiting a much lower SE value, has the third highest ME value, just below those of T. taxus and L. primus. The similarity observed between the ME-SE bite profiles of L. primus and T. taxus is examined further via analysis of von Mises stress distributions on skull models.

In terms of general distribution of stresses throughout the skull, Taxidea taxus and Leptarctus primus exhibit highly similar von Mises stresses in nearly every tooth bite position (Fig. 4). Significant stress levels and distributions around the postorbital constriction, the posterior and anterior sections of the zygomae, and the top of the cranium of the individual are repeatedly and similarly observed in both taxa. Stress distribution in Po. flavus is more concentrated around the working side zygoma and postorbital constriction compared with L. primus and T. taxus, and the temporomandibular joints (TMJs) of Po. flavus exhibit higher stress as well. Lycaon pictus has very low overall von Mises stress, with slightly elevated stress levels around the TMJs and the orbits. Gulo gulo shows the highest level of stress throughout the skull, with the most widespread elevated stress (green to red shades on the heat maps), especially around the postorbital constriction, TMJs, orbits, and midline at nasal bones (Fig. 4).

Hierarchical Cluster Analysis with Resampling and k-means Analysis

In addition to inferring the feeding capabilities of L. primus via analyses of unilateral bite forces and their visualizations through von Mises stress, cluster analyses were used to determine similarities in numerical FE values. Clusters containing L. primus signify biomechanical similarity between those taxa that Leptarctus pairs with consistently. Of the eight clustering methods employed, the cluster that contains L. primus with the highest level of bootstrap support was the one using the centroid cluster method, containing a Ly. pictus–Po. flavus pair, which is in turn paired to L. primus (Fig. 5). A similar cluster occurred using the average cluster method. Leptarctus primus did not pair with other species in two of the analyses, using the medianand single-cluster methods, respectively (Fig. S1). In one of the cluster methods (complete), L. primus paired with Po. flavus (84% AU, 36% BP), and in three of the cluster methods (Ward D, Ward D2, and Mcquitty), L. primus paired with T. taxus (44–62% AU, 41–52% BP) (Fig. S1).

TABLE 1.

Results of k-means analyses.

FIGURE 3.

A, plot of mechanical efficiency (ME) versus adjusted strain energy (SE) of unilateral bite simulations at all tooth loci present on all 21 species FE models. Darker-shade (red) solid curves are ecological carnivores; thin dotted (green) curves are omnivores; light (blue) solid curve represents the frugivore Po. flavus; thick dashed dark curve represents the herbivore A. fulgens. The extinct musteloid Leptarctus primus is represented by the thick dark solid curve. B, comparison between extant T. taxus (darker-shaded red curve) and extinct L. primus (thick black curve) shows the similarity in SE patterns and higher ME values in T. taxus across the tooth row. C, the comparison between L. primus (solid black curve) and a modified L. primus model (dotted curve; as a proxy for the likely conspecific Hypsoparia bozemanensis) with deeper zygomae indicates that although adjusted strain energy values are similar, the model with deep zygomae exhibits increased mechanical efficiency at every tooth locus compared with the original L. primus model. Error bars represent 95% confidence intervals within species for which multiple resolution FE models were available.

Among the 19 k-means analyses (from K = 2 to K = 20), L. primus was classified as a standalone ‘cluster’ (of one species) 57.89% of the time (Table 1); 31.58% of the time L. primus clustered with T. taxus, 21.05% of the time with Po. flavus, and 10.53% of the time with Ly. pictus. The rest of the taxa sampled did not cluster with L. primus at any of the K values tested.

FIGURE 4.

von Mises stress distributions in crania of extinct Leptarctus primus and extant Taxidea taxus, Potos flavus, Lycaon pictus, and Gulo gulo. Unilateral bites are shown for the right upper tooth row, with C (canine), P2 (second premolar), P3 (third premolar), P4 (fourth premolar), and M1 (first molar). Cooler (darker) shades indicate lower stress; warmer (lighter) shades indicate higher stress.

FIGURE 5.

Hierarchical and resampling cluster analysis. Dendrogram based on ME and SE values using the centroid cluster method shows a highly supported cluster of L. primus clustering with Ly. pictus and Po. flavus. Node values indicate approximately unbiased (AU) bootstrap probability. For dendrograms of analyses employing other cluster methods and their bootstrap support, see Figure S1. Superscripted asterisks indicate fossil species. Abbreviation: U ., Urocyon.

Phylogenetic Signal and Phylogenetic Principal Components Analyses

For combined ME and SE values of the whole data set, and ME and SE outputs separately, Blomberg's K values ranged from 0.33 to 0.78 among the six branch length configurations, with P-values from 0.14 to 0.69 (Table 2). Mantel tests were run on the same data partitions (ME+SE, ME, SE). The Z statistic values ranged from 3,111 to 67,387 and P-values from 0.07 to 0.99 across all six branch length configurations (Table 3). Based on these results, the FEA data do not include significant amounts of phylogenetic signal among the specimens tested.

Phylogenetic PCA using the six different branch length configurations returned PC axes that accounted for 54.8–65.8% of total variation along PC1 and 32.7–43.9% of total variation along PC2 (Table 4). The first two PC axes together accounted for more than 98% of total variation in all six configurations. Principal Component 1 in all cases is heavily weighted by SE values, and in two of the configurations ME values contribute equal weight to PC1 values. Negative PC2 values are strongly influenced by ME (higher ME values on more negative PC2 values), and positive PC2 values are less strongly influenced by SE values (higher SE values on more positive PC2 values) (Figs. 6, S2). In five out of six branch length configuration pPCAs, Po. flavus falls on the extreme end of the first PC axis with the highest SE values or a combination of higher SE and high ME values. In all of the configurations, T. taxus falls on the most negative extreme of the PC2 axis by having the highest ME values.

Linkage between Dietary Category and Biomechanical Simulation Output

Both univariate and multivariate Kruskal-Wallis tests returned statistically nonsignificant results for the ME and SE values between the carnivore and omnivore dietary categories (Table 5). The ME-SE profiles of carnivorous and omnivorous species have overlapping SE ranges, but the highest SE values are observed in the omnivorous U. cinereoargenteus. The ME value ranges also partially overlap between species of carnivores and omnivores, but the highest ME values (>0.40) are observed only in carnivorous species (Ly. pictus, T. taxus, G. gulo), plus the frugivorous Po. flavus. The herbivorous A. fulgens has ME-SE values deeply nested within a group of carnivores and omnivores, characterized by relatively low ME and low SE values (Fig. 3).

Morphological Variation of H. bozemanensis Proxy and L. primus Models

Overall, the composite proxy model representing H. bozemanensis has a higher ME value than the model for L. primus for every tooth locus sampled (C, P2, P3, P4, M1) and a marginally higher SE value in the canine (C) (Fig. 3C). However, only the ME differences are significantly different when the 95% CIs are considered. The SE values have ranges of uncertainty that overlap substantially between the corresponding tooth loci tested (Fig. 3C).

The bite profiles between the two individual models are similar. In both models, the canine retains the highest SE values. Hypsoparia bozemanensis lacks the noticeable dip in P4 SE values characteristic of hypercarnivores (Tseng and Flynn, 2015b) but mirrors the increasing biomechanical outputs from P4 to M1 in Leptarctus primus. The only significantly different SE values are between the canine and P2 bite positions in both models; the premolar and molar bite positions show an overall, non–statistically significant increase (as evidenced by overlapping 95% CIs of SE values) from front to back of the tooth row.

All model files used or generated in this study are freely available online at the URLs listed in Table S5.

Inferring the Diet of Leptarctus

As detailed in Introduction, Leptarctus primus had been variably inferred to be a carnivore, frugivore, herbivore, insectivore, or omnivore. However, based on its placement within the MESE plot for the 21 species sampled (Fig. 3), its clustering tendencies from hierarchical and k-means analyses (Figs. 5, S1; Table 1), and pPCA (Figs. 6, S2), our results suggest that the primary feeding capability in Leptarctus primus most likely was carnivory, with a secondary capability of omnivory/opportunistic omnivory. These results are consistent with the only previous quantitative analysis of the potential feeding adaptation of leptarctines (Calede et al., 2017), which used different, simpler methodologies to arrive at a similar interpretation. Prior hypotheses of insectivory, frugivory, and herbivory are not supported as the primary dietary paleoecology of Leptarctus primus.

The results from the ME-SE plot indicate that L. primus has similar biomechanical capability to some extant carnivores, with a bite profile most similar to the carnivorous badger Taxidea taxus: the canine bite returns the greatest strain energy (J), the P2 bite exhibits a drop in SE values, and there is an increase in SE (making the first molar comparable to the canine in strain energy) toward the back of the tooth row. The high mechanical efficiency of the entire bite profile relative to other taxa modeled suggests that L. primus and T. taxus both had the capability of producing higher bite forces relative to size compared with other species, and this is corroborated by the k-means analyses (Table 1). The diet of Taxidea taxus is composed mostly of a variety of small vertebrates such as rodents (Table S1). Although the bite profiles are nearly identical in shape, the significantly greater M1 mechanical efficiency of T. taxus suggests that it has a higher molar crushing ability than did L. primus, which may be associated with its omnivorous secondary dietary capability. Furthermore, whereas T. taxus has a much lower adjusted strain energy, SE values for L. primus (Table S3) are much closer to those of other tested caniforms, such as Canis and Ursus (Tseng and Flynn, 2015b). This suggests that L. primus had a relatively low cranial stiffness, more like some caniforms rather than the stiffened, feliform-like cranium of T. taxus. von Mises stress distributions indicate that among the four extant species most commonly associated with L. primus in the various multivariate analyses, the L. primus model is most similar to T. taxus in each tooth locus. Potos flavus has a skull with more concentrated stresses around the working (biting) side orbit and both TMJs, Ly. pictus has very low stresses overall, and G. gulo has relatively high stresses throughout the skull (Fig. 4). Insectivory as a primary dietary capability is highly unlikely when comparing the biomechanical output of L. primus with that of the three primary omnivore/insectivores, Spilogale putorius, Herpestes javanicus, and Mephitis mephitis (see Table S1). The more bunodont dentition of Leptarctus relative to insectivorous carnivorans is consistent with this inference from the biomechanical data. Low but posteriorly increasing ME and SE values for Spilogale indicate a small but steady increase in biomechanical efficiency across the tooth row (C–M1). In comparison, the values for Leptarctus differ in that they have much higher ME and SE values, thereby more closely resembling primary carnivores. This does not preclude L. primus from having sustained itself partially on insects (see ‘omnivory/opportunistic omnivory,’ as defined in Table S1), but given Leptarctus’ significantly higher mechanical efficiency, it is highly unlikely that it was limited to or had a primary dietary specialization for insectivory as characterized by the species tested in this study. This statement operates on the assumption that individuals with higher ME and SE capabilities would rely on food sources that require this increase in ME-SE values, although it is possible that, like some extant, highly muscled carnivorans (see Ursus), it would consume prey or food items that do not require the full cranial biomechanical capability of the predator.

TABLE 2.

Results of multivariate Blomberg's K statistic for testing phylogenetic signal.

In sum, ME and SE data indicate that L. primus was primarily carnivorous and secondarily omnivorous. It is important to note that the FE models were scaled so that muscle force is proportional to muscle attachment area, and as such, differential body size is partially corrected for in this approach. So, unless the larger mustelids also tend to show muscle attachment areas that are larger than already expected by this scaling relationship, size alone would not be an explanatory factor for the similarities exhibited here.

Dendrograms generated by hierarchical cluster analyses indicate that L. primus variably groups with T. taxus, Ly. pictus, or Po. flavus, reflecting the distributions observed in the pPCA (Figs. 5, 6, S1, S2). Because both the ME-SE profile and kmeans analyses of L. primus further support a close association with T. taxus, we conclude that the American badger is the best extant biomechanical analog for L. primus, with Lycaon pictus, Po. flavus, and G. gulo as secondary analogs. These results corroborate the conclusions of Olsen (1957) more than 60 years ago, who suggested that T. taxus could serve as an extant analog for Leptarctus, and the more recent interpretation of carnivory by Calede et al. (2017) based on dry-skull bite force estimates. Taxidea taxus is a primary carnivore, the diets of G. gulo and Ly. pictus consist mostly of hunted or scavenged vertebrate animal matter, and Po. flavus is a frugivore (Table S1). This seemingly divergent connection to both carnivores and frugivores is driven by the intermediate ME and SE values of L. primus between the similarly high-SE and lower-ME bite profiles of Po. flavus and Ly. pictus on one end, and the similarly lower-SE and higher-ME bite profiles of T. taxus and G. gulo on the other (Table S3). Both a fruit-based diet (as in Potos) and a hypercarnivorous diet with low bone consumption (as in group-hunting Lycaon) are considered to be dietary specializations with low mechanical demand compared with obligately plant-based (e.g., tough fiber diets such as bamboo in Ailuropoda and Ailurus) and bone-based hypercarnivore diets like in Crocuta (Figueirido et al., 2013). On the other hand, the mustelids T. taxus and G. gulo are solitary opportunistic hunters of a variety of food items that may vary seasonally and which may require great bite force by a single animal for killing and completely consuming small mammal prey or scavenging tough carcasses (Long, 1973; Pasitschniak-Arts and Larivière, 1995). The intermediate position of L. primus suggests a distinctive dietary specialization and may indicate either that it was capable of a combination of those behaviors to some extent, or that it was not as specialized or generalized in diet as the two most similar biomechanical analog groups identified by our analyses. Although it clusters with both the carnivorous and frugivorous species mentioned above, L. primus shows similarities in SE values to both of the less mechanically efficient ursids (Ur. arctos, Ur. maritimus) as well as the blackbacked jackal (Ca. mesomelas). Its deeper zygomatic arches and the associated increased area for masseter muscle attachment would enable Leptarctus to bite down quickly and hold on tightly, potentially similar to behavior observed in extant small mustelids such as species of Mustela (Ewer, 1973). The laterally displaced arrangement of parasagittal crests in L. primus also might explain the disparity in skull stiffness (measured by strain energy) between L. primus and T. taxus. Taxidea taxus has a narrow, low-vaulted sagittal crest with more muscle attachment area over the surface of the cranial vault, whereas L. primus has a more dorsoventrally oriented lateral cranial vault capped by strong parasagittal crests and exhibits elevated stresses on the dorsal cranium, most likely from tensile forces pulling the dorsal cranium laterally during temporalis muscle contraction (Fig. 4).

TABLE 3.

Results of Mantel tests for phylogenetic covariation of biomechanical simulation data.

TABLE 4.

Percentages of total variation accounted for by the first two principal component axes from phylogenetic principal components analyses (pPCAs).

In hierarchical cluster analyses, L. primus is widely separated from the herbivorous Ailurus fulgens (which falls within a distant cluster). Based on this result, it is unlikely that Ailuruslike herbivory was the main feeding capability for L. primus. The association between L. primus and Po. flavus in some cluster analyses and the pPCA suggests that the secondary omnivory capability of L. primus may have included some frugivory. However, the similarities in ME-SE profiles between Po. flavus and Ly. pictus, the latter a hypercarnivorous canid, suggest either (1) that ME-SE profiles alone may not be sufficient to distinguish frugivores from other types of feeders or (2) that frugivores do not differ significantly in cranial biomechanical capability from some carnivores as measured by the FE method employed in this study. Consistent with the potential inability of tooth row ME-SE profile values to distinguish dietary categories in larger taxonomic samples, the nonparametric Kruskal-Wallis test returned no significant differences in ME, SE, or ME plus SE values between carnivorous and omnivorous species (Table 5).

These results raise another question: with so much prior attention paid to the dental morphology and craniodental characters of the L. primus neotype (AMNH FM 18241, AMNH FM 18270), why have most previous studies not considered carnivory as the most likely dietary specialization for this taxon? The L. primus neotype—upon which most of the prior L. primus dietary inferences were made—has considerable wear on its entire tooth row. A lack of apical wear on all tooth cusps of the fully erupted dentition in the L. primus (AMNH FM 25385) specimen analyzed here suggests that this specimen was ontogenetically younger than most or all other described Leptarctus specimens. This contrasts greatly with the L. primus neotype, whose cusps are heavily worn, with substantial dentine exposure. It is probable that the extensive tooth wear present on most Leptarctus specimens led early workers to interpret a noncarnivorous diet, because a reduction in (or total lack of) shearing and piercing facets on the premolars would render carnivory highly inefficient (e.g., trenchant cusps in carnivorous felids versus bunodont cusps in omnivorous ursids). However, based on the young adult specimen of L. primus analyzed in this study (AMNH FM 25385), a biomechanical capability for carnivory in L. primus not only is possible but also is well supported by the results of our broad comparisons across Carnivora.

FIGURE 6.

Bivariate plot of the first two principal components in the phylogenetic principal components analysis using uniform branch length configuration. 1, Ca. lupus; 2, L. primus; 3, G. gulo; 4, S. putorius; 5, T. taxus; 6, H. javanicus; 7, M. mephitis; 8, O. herpestoides; 9, Pa. pardus; 10, P. lotor; 11, Th. velox; 12, A. fulgens; 13, C. crocuta; 14, Ly. pictus; 15, Ur. maritimus; 16, Ur. arctos; 17, Ca. mesomelas; 18, Pr. brunnea; 19, B. astutus; 20, U. cinereoargenteus; 21, Po. flavus.

TABLE 5.

Results of nonparametric statistical test (Kruskal-Wallis test) for difference between carnivorous and omnivorous species in the data set.

In sum, we found little evidence via FEA that L. primus was solely specialized for herbivory, frugivory, or insectivory. Instead, this analysis provides compelling evidence that L. primus was capable of eating a diet resembling that of the American badger: a largely carnivorous diet of small vertebrates, snakes, and birds, supplemented secondarily by opportunistic omnivorous feeding on insects, eggs, seeds, fungi, etc. Given the close association between L. primus, T. taxus, and Gulo gulo in the pPCA plots (Figs. 6, S2), further studies comparing L. primus with additional badger ecomorphs within Mustelidae would be worthwhile in further testing this dietary paleoecological inference.

Phylogenetic Signal

Vertebrate paleontologists have long used morphology and noncomputational assessments of functional biomechanics as proxies for determining paleodiet. This study expanded on these approaches by applying computational biomechanics methods and statistical analyses across broad taxonomic and ecological samples of taxa to infer cranial properties and dietary preferences, after testing for or accommodating phylogenetic signal as a potential primary explanator of shared characteristics. Phylogenetic signal in this paper determines the dependency ofMEand SE values on a specimen's phylogenetic relationship to other specimens in the study (see Table 2, Blomberg's K statistic values and P-values). Every Blomberg's K statistic value is less than 1, indicating that there is less phylogenetic signal than expected within the sample sets. Based on this, we conclude that specimens are not clustering based on significant phylogenetic influence on their cranial biomechanical properties. This conclusion is further supported by results of Mantel tests of distance matrices indicating that phylogenetic relatedness (e.g., between L. primus and some of its closest biomechanical analogs [mustelids, procyonids, or canids]) did not significantly influence the associations observed in the data.

The results from these phylogenetic signal analyses, although implying that the association between Po. flavus, T. taxus, Ly. pictus, G. gulo, and L. primus is not phylogenetically significant, do not necessarily ensure that they are more closely linked because of functional evolution related to masticatory mechanics. Carnivoran cranial evolution has many driving factors, with a species' ability to efficiently accommodate stresses as they bite into prey items being only one of them (e.g., Tseng and Flynn, 2018). This is a topic that should be explored further through evaluation of other potential correlates with cranial biomechanics, as well as direct tests for covariance between biomechanical attributes and diet, and reexamined for past studies that proceeded on this same assumption that craniodental form is directly—in some cases, solely—caused by masticatory pressures for specific diets or feeding strategies.

Ecological Signal in FE Models of Extant Species

The clustering of taxa with mixed feeding ecologies within the L. primus cluster also is observed in the rest of the dendrograms (Figs. 5, S1). For example, the durophagous hypercarnivores C. crocuta and Pr. brunnea and durophagous obligate herbivore A. fulgens cluster with the small-bodied carnivore/omnivore B. astutus and the large-bodied carnivore Pa. pardus. Similarly, although the highly carnivorous, large-bodied Ca. lupus and Ur. maritimus cluster together, they are joined by the smallerbodied omnivores M. mephitis and P. lotor (Fig. 5). This mix of biomechanical attributes between taxa of different body sizes and feeding ecologies (e.g., large hypercarnivores with smaller omnivores) suggests that body size and its allometric effects may be additional important factors in distinguishing feeding ecologies using biomechanical traits. Smaller-bodied omnivores may face similarly stringent biomechanical demands when feeding on prey that are close to their own body size as do large-bodied carnivores feeding on larger prey. Because simulated biomechanical attributes analyzed in this study are sizestandardized, and comparative muscle data are lacking to enable determination of size-specific model parameters, any additional diagnostic biomechanical traits that are size-related were not included in these results. That all extant species closely clustering with L. primus have larger skulls suggests that L. primus possessed biomechanical capabilities expected of carnivorans much larger than its actual size would otherwise indicate. This interpretation further supports a potentially unique ecological niche for L. primus that combined the capabilities of its closest extant analogs or represented a distinctive feeding ecology that is not represented by any of the extant species included in this analysis. Additional quantitative analyses are required to test and distinguish between these alternative explanations, potentially testing muscle-specific biomechanical capability and accounting for influence of nonfeeding ecological traits on cranial biomechanical characteristics (Tseng and Flynn, 2018).

Comparison of Hypsoparia and Leptarctus Cranial Biomechanics

The biomechanical analyses on the skull of Leptarctus primus and the Hypsoparia bozemanensis proxy model showed consistent increases in mechanical efficiency when the zygomae are deepened in the Hypsoparia proxy model, and no statistical differences in strain energy values at analogous tooth loci (Fig. 3). Given the bite profiles of each of those models, it is possible that adaptations within this species (if the two were conspecific), or in the closely related taxa if they were distinct species, varied to some degree, with morphologies indicating diets ranging from more carnivorous to more omnivorous. These observed differences may represent geographical variation in diet and skull biomechanics within a single species, or they may indicate an ecomorphological transition from more carnivory to less carnivory among closely related species within the genus, as proposed by Calede et al. (2017) for leptarctines overall. Although Ly. pictus and Po. flavus both represent specialized dietary ecomorphs (hypercarnivory and frugivory, respectively), dietary flexibility is observed in Gulo gulo (a primary carnivore and secondary opportunistic omnivore, also known as an opportunistic scavenger, its name translating to ‘Glutton glutton’) as well as in the carnivore Taxidea taxus, whose primary diet is comprised of vertebrate prey, with some invertebrate prey (minor omnivory or insectivory) depending on ecological factors (Potter, 1924; Errington, 1937; Grinnell et al., 1937; Hamilton, 1939; Snead and Hendrickson, 1942; Drake and Presnall, 1950; Rausch, 1959; Jackson, 1961; Rausch and Pearson, 1972). This carnivore/frugivore specialist versus omnivore generalist dichotomy corresponds to high-SE and low-ME versus low-SE and high-ME bite profiles, respectively (Table S3). Thus, a flexible, generalized omnivorous diet appears to be associated with higher mechanical efficiency and a stiffer skull compared with specialized carnivorous or durophagous herbivore diets, potentially to enable consumption of tough or hard foods that may be part of the diet. The intermediate ME and SE values suggest that L. primus was not as specialized as the extant frugivores or hypercarnivores, nor did it have the biomechanical properties in its skull that would permit processing extremely tough and/or hard foods to the extent shown by the carnivores G. gulo or T. taxus. Using the dietary preferences of these extant carnivorans as guides gives further insights into leptarctine dietary adaptations, within species or across closely related species, and permits more refined paleoecological reconstruction of these distinctive and previously enigmatic fossil musteloid carnivorans.