Clarifying relationships between cranial form and function in tapirs, with implications for the dietary ecology of early hominins

Paleontologists and paleoanthropologists have long debated relationships between cranial morphology and diet in a broad diversity of organisms. While the presence of larger temporalis muscle attachment area (via the presence of sagittal crests) in carnivorans is correlated with durophagy (i.e. hard-object feeding), many primates with similar morphologies consume an array of tough and hard foods—complicating dietary inferences of early hominins. We posit that tapirs, large herbivorous mammals showing variable sagittal crest development across species, are ideal models for examining correlations between textural properties of food and sagittal crest morphology. Here, we integrate dietary data, dental microwear texture analysis, and finite element analysis to clarify the functional significance of the sagittal crest in tapirs. Most notably, pronounced sagittal crests are negatively correlated with hard-object feeding in extant, and several extinct, tapirs and can actually increase stress and strain energy. Collectively, these data suggest that musculature associated with pronounced sagittal crests—and accompanied increases in muscle volume—assists with the processing of tough food items in tapirs and may yield similar benefits in other mammals including early hominins.


Finite Element Model Construction
All CT data was collected from medical CT scanners and exported in DICOM format. The resulting images were imported into Avizo v9.0 (Visage Imaging, Inc.), where the cranium and mandible were selected using a combination of automated thresholding based on gray values, and they were manually edited to separate the cranium from the mandible and generate 3D surface models. The surface reconstructions were then smoothed and edited to improve the quality of the mesh, including testing the aspect ratio and the dihedral angles of the surface triangles to minimize distorted elements which can cause errors. The aspect ratios of the triangles were adjusted to below 10 (aiming for around 5) and the dihedral angles were set at above 10 degrees to ensure a good quality mesh. The surface meshes were then converted to solid 3D finite element meshes composed of 4-node tetrahedral elements above 1 million elements (supplemental table 2). Finally, each model was exported as an Abaqus input file (*.inp) for easy importation to the FEA software package Abaqus CAE v6.14 (Simulia). The cranium and mandible were imported separately for each species, so that the mandible could be used for muscle force alignment, but excluded from the analysis.
The jaw adductor musculature was modeled as three main components; the temporalis, masseter and pterygoid muscles. Given that there are no measured recordings of bite force or muscle force in tapirs from dissections or in vivo experiments, muscle input forces were calculated by multiplying the physiological cross-sectional area (PCSA), estimated using the dry skull method, by 0.3 N, the maximum force of contraction for mammalian muscle fiber [1]. The dry skull method estimates the PCSA of the combined masseter and pterygoid system, so approximately 10-15% of the total muscle force was assigned to the pterygoid for each species. Muscle forces were applied to the skulls by distributing the load for each muscle over the entire surface of the muscle origin. Muscle force orientations were determined by creating a local coordinate system with a vector between the origin and the corresponding insertion on the mandible. A coupling constrain was used to represent the fan shaped nature of the muscle fibers of the temporalis.
Each model was then imported to Abaqus CAE v6.14 where material properties and boundary conditions were applied. To prevent the models from rigid body motion, a single node was constrained at both temporomandibular joints (TMJ) and at each premolar (P2-P4) or molar (M1-M3) on the left, right or both sides for unilateral or bilateral biting, respectively. The left TMJ was fully constrained against displacement in the x-(lateral), y-(vertical), and z-direction (anterior-posterior), and the node on the right TMJ was constrained in the y-and z-axis to allow lateral displacement of the skull. The bite point(s) were constrained in the axis perpendicular to the occlusal plane, so as not to over constrain the model. A number of biting scenarios were modeled to simulate bilateral biting at each individual tooth along the tooth row, and unilateral biting at each tooth with the balancing side muscles adjusted to 60% of the working side. Both left and right unilateral biting was modelled. Even though single node constraints can create artificially high stress and strain in the elements surrounding the constraint, overall patterns of stress and strain distribution are not affected [2,3]. As our goal was to examine overall stress and strain distributions in the skull, with a focus on the function of the sagittal crest not the region around the constraints, these single node constraints are acceptable and minimize the complexity of the model. However, to account for artifacts at these constraints, we examined both the maximum stress and an adjusted maximum in which we removed the top 2% of stress values.
In the absence of material properties data for tapirs, and due to incorporating a fossil specimen with which material properties cannot be estimated from the CT scans, each model was assigned homogeneous and isotropic average values of Young's modulus (E = 20 GPa) and Poisson's ratio (ν = 0.3) for mammalian bone [4]. This may produce stiffer models than models with multiple material properties, however, this methodology was considered suitable for the present study, which compares only relative stress and strain values.
Due to the variation in size among the tapir species, the models were standardized to the same size for comparison of shape alone. Stress and strain will depend on the size of the structure. For example, a smaller structure will be stiffer than a larger one, whereas a larger structure will be stronger.
All analyses were performed on the original size models and these results were scaled so that the surface area or volume of each model was equivalent to T. terrestris. To compare the stress (strength) between models, the models were scaled to the same ratio of muscle force to surface area. When comparing strain energy, which is the amount of energy stored per unit volume, the models were scaled to the same force: volume ratio, following Dumont and colleagues [5].

Enamel Thickness Estimates.
Avizo was used to non-orthogonally reslice teeth along the protoloph and metaloph orientations (note, protoloph and metalophs were not parallel); measurements were taken on the resliced scans. All measurements were taken on the left upper second molar (M2), except for T. pinchaque, where the right M2 was measured (the left being damaged). Measurements were taken on the protoloph and the metaloph (see supplemental figure 1). As CT data for T. polkensis showed no contrast between dentine and enamel, and also represented an old individual with visibly worn teeth, it was not measured. All other specimens had the M3 in, or partially in a crypt (the reason for measuring the M2). Since estimating enamel thickness was not a primary aim of our study we did not collect microCT data of the teeth in addition to the lower resolution CT scans used for FEA, therefore these measurements should only be considered estimates. Enamel thickness is reported as raw values, and two metrics were used to remove the effect of body size (as estimated from the width of the protoloph and metaloph of M2s). Specifically, enamel thickness on the protoloph and metaloph were divided by the width of the protoloph or metaloph, respectively. Additionally, enamel thickness was scaled to the size of the largest tooth width (T. terrestris in both cases) by dividing smaller teeth by the tooth width of T. terrestris, subsequently multiplying the product by the raw enamel thickness values.

Finite element analysis
Unilateral biting. The performance for each model at each bite location for unilateral biting is presented in supplemental table 2. For unilateral biting the mean stress, adjusted maximum stress, and strain energy for each model decreased along the tooth row from anterior to posterior until approximately half way at which point it increased for all models. Tapirus bairdii, however, increased along the entire tooth row for stress and strain, indicating T. bairdii is stronger and more energy efficient at more anterior bite points.
The highest stress was recorded in T. pinchaque for the most anterior bite points, but was higher in T.
bairdii for the posterior bite points. Mechanical efficiency increased for all models as the bite point moved posteriorly. T. bairdii has the highest mechanical efficiency and T. pinchaque the lowest.
Bilateral biting. The performance of each model at each bite location for bilateral biting is presented in supplemental table S3. Variation between unilateral and bilateral biting was most notable when comparing stress. For bilateral biting the highest stress was recorded in T. polkensis and the lowest was recorded in T. bairdii for adjusted maximum and mean stress, but was lowest in T. terrestris for maximum stress. However, for unilateral biting T. polkensis has the lowest adjusted maximum stress.
Strain energy did not vary much between bilateral and unilateral biting and neither did mechanical efficiency. (supplemental table 1 ); Asfc, area-scale fractal complexity; epLsar, anisotropy; Tfv, textural fill volume; HAsfc 3x3 , HAsfc 9x9 , heterogeneity of complexity in a 3x3 and 9x9 grid, respectively. *= low parasagittal ridges that do not unite to form a sagittal crest in adults, or **=sagittal crest present in T. polkensis at the Gray Fossil site.