DOI: 10.5061/dryad.d2547d86g
title: The trunk replaces the longer mandible as the main feeding organ in elephant evolution--Supplementary codes

description: The long-trunked elephantids underwent a significant evolutionary stage
 characterized by an exceptionally elongated mandible. The initial
 elongation and subsequent regression of the long mandible, along with its
 co-evolution with the trunk, present an intriguing issue that remains
 incompletely understood. Through comparative functional and
 eco-morphological investigations, as well as feeding preference analysis,
 we reconstructed the feeding behavior of major groups of longirostrine
 elephantiforms. In the Platybelodon clade, the rapid evolutionary changes
 observed in the narial region, strongly correlated with mandible and tusk
 characteristics, suggest a crucial evolutionary transition where feeding
 function shifted from the mandible to the trunk, allowing proboscideans to
 expand their niches to more open regions. This functional shift further
 resulted in elephantids relying solely on their trunks for feeding.
description: <strong>Cladistic analysis</strong>
 Cladistic analyses were performed to evaluate the phylogenetic
 hypothesis of trilophodont longirostrine proboscideans. The data matrix
 contained 37 taxa, including most of the known trilophodont longirostrine
 taxa at the species level, and <em>Phiomia serridens</em>, an
 Oligocene basal elephantiform, was selected as outgroup. Additionally, a
 basal elephantoid, <em>Tetralophodon longirostris</em> was
 also included to assess which clade the true elephantids originated from.
 The morphological characters included 5 characters from upper tusks, 9
 from mandibular tusks, 37 from cheek teeth, 19 from the cranium, and 10
 from the mandible, mainly following Tassy, 1996; Shoshani, 1996; and Wang
 et al., 2017. For a description of the characters and states, please see
 Appendix; for the data matrix, please see Figure 1–Source Data 1. Two
 methods, Bayesian tip-dating (BTD) and maximum parsimony (MP) analyses
 were performed. In BTD analysis, the fossil ages were
 incorporated as tip calibrations (<em>Gavryushkina et al.,
 2014</em>; <em>Ronquist et al., 2012</em>;
 <em>Zhang et al., 2016</em>). The Lewis Mkv model
 (<em>Lewis, 2001</em>), with gamma rate variation across
 characters<sup> </sup>(Mkv + G) (<em>Yang,
 1994</em>), was initially used; subsequently, the timetree was
 modelled by the fossilized birth death process<sup>
 </sup>(<em>Heath et al., 2014</em>; <em>Stadler,
 2010</em>). The process was conducted using the time of the most
 recent common ancestor (root age) and included hyperparameters of
 speciation rate, extinction rate, fossil-sampling rate, and
 extant-sampling probability. The root age was first assigned an
 offset-exponential, with mean age of 37 Ma and minimum age of 34 Ma,
 referring to the oldest fossil. The fossil ages were fixed to their first
 occurrence. The extant-sampling probability was fixed to 1 because no
 living genera were specified. Apart from the timetree, the other key
 component was the relaxed clock model, which models the evolutionary rate
 variation along the branches in the tree. We used the independent gamma
 rate clock model (<em>Lepage et al., 2007</em>), in which the
 mean clock rate was initially assigned a lognormal prior (–6, 1) and the
 variance parameter of the clock rate was exponential (10).
 We executed two independent runs and four chains per run (one
 cold chain and three hot chains) using Markov chain Monte Carlo. Each run
 was executed for 1 million generations and sampled every 5000 generations.
 The first 25% of the samples were discarded as burn-in and the rest of the
 two runs were combined. Good convergence and mixing were determined by
 effective sample sizes<sup> </sup>larger than 200 for all
 parameters and average standard deviations of split frequencies smaller
 than 0.01 (<em>Geyer, 1992</em>). The BTD analysis was
 performed in MrBayes 3.2.7 (<em>Ronquist et al.,
 2012</em>)<sup> </sup>(<em>Supplementary file
 3</em>). MP reconstruction was performed by
 TNT1.1 (<em>Goloboff et al., 2008</em>). In MP analysis, all
 characters were equally weighted. Characters 20–25 pertained to the
 loph/lophid numbers of cheek teeth (character numbers begin from “1” in
 our numeration; however, in MP analysis, TNT1.1 automatically numbered
 characters from “0”; therefore, in the TNT1.1 program, these characters
 were numbered “19–24”). These were treated as ordered and irreversible,
 which was performed by setting “step-matrix → of costs” under the menu
 “Data → Character settings”, assigning the value “9” to the blanks where i
 &gt; j (i to j in the matrix), and assigning the value “j–i” to the
 blanks where i &lt; j. The traditional search strategy was performed,
 and the results were reported based on a 50% majority consensus tree of
 the most parsimonious trees. Node supports were calculated by symmetric
 resampling with 0.33 change probability (1000 replicates). The script
 files are provided in Supplementary file 3. The major consensus tree was
 then calibrated using the time PaleoPhy function in R package paleotree
 3.3.25 (<em>Bapst, 2012</em>), and is shown in Figure 1–figure
 supplement 1. <strong>FE
 analysis</strong> We investigated the feeding
 behaviors of the three longirostrine gomphothere families, i.e.,
 Choerolophodontidae, Amebelodontidae, and “Gomphotheriidae”, using finite
 element (FE) stimulation. Three species were selected to represent each
 family, <em>Choerolophodon</em><em>
 chioticus</em>, IVPP V23457 (cranium and mandible);
 <em>Platybelodon grangeri</em>, HMV 0930 (cranium and
 mandible); and <em>Gomphotherium tassyi</em>, IVPP V22780
 (cranium) and IVPP V22781 (mandible), respectively. Note that, in
 <em>Gomphotherium</em>, we were unable to get access to
 cranium and mandible that belonged to one individual. We used a handheld
 Artec Spider 3-dimensional scanner to obtain the surface topology of these
 specimens. The surface meshes were produced using Artec Studio 14
 Professional. These meshes were first repaired with ZBrush 2021; for
 example, ZBrush 2021 was used to recover the broken edges of mandibular
 tusks, create the remaining mandibular tusks in the alveolus, create the
 keratinous cutting plate of the mandible in
 <em>Choerolophodon</em>, and for retro-deformation of the
 crushed cranium (i.e., the cranium of IVPP V23457). The
 rough surface meshes were edited using Materialise 3-matic Research
 (V12.0). The surface meshes were smoothed by removing the small knobs and
 filling the holes for volume mesh generation. Before volume meshing, the
 models were cut into symmetric halves along the median sagittal plane. To
 reduce computation, only the right halves were preserved for further
 analyses. Note that the cranium was only used to define the attachments or
 insertions of jaw-closing muscles (modeled by many draglines, see below),
 and was treated as a rigid body in simulation (therefore, the cranial and
 nasal cavities are not relevant). The mandible contains two parts, the
 bony structure, including the cheek teeth, and the food acquisition organs
 (the mandibular tusk in <em>Platybelodon</em> and
 <em>Gomphotherium</em>, and the keratinous cutting plate in
 <em>Choerolophodon</em>). These two parts were integrated
 using the command “Non-manifold Assembly” in Materialise 3-matic Research.
 Then, the cranium and the mandible were aligned based on their natural
 position, i.e., the occlusal surface of the upper and lower tooth rows
 were matched, and the mandibular condyle and glenoid fossa were fitted.
 The sagittal surfaces of the cranium and mandible were made to coincide
 with the <em>x</em>-<em>z</em> plane, and the
 <em>x</em>-positive direction was set along the rostral
 direction. For comparison between different taxa, the three mandibular
 models were scaled to the same volume as
 <em>Choerolophodon</em> (6,563,708.0146
 mm<sup>3</sup>) (<em>Figure 3–figure
 supplement</em><em> 1</em>). Finally, volume meshes were
 generated in the cranium and mandible across the three models, and
 exported as .inp files that could be loaded into Abaqus CAE (V 6.14), the
 engineering software for FE analysis. The parameters for volume mesh
 generation are listed in Supplementary file 1. The
 volume meshes, representing the geometric models, were imported into
 Abaqus CAE (V 6.14) (<em>Figure 3–figure
 supplement</em><em> 2–4</em>), and included two parts,
 cranium and mandible; additionally, a third part was created by Abaqus
 CAE, a long cylinder (300 mm in length, 50 mm in diameter), to model twigs
 that were cut by mandibular tusk or keratinous cutting plate. Note that
 the millimetre (mm)–ton (t)–second (s) unit system was adopted; other
 units included Newton (N) and million Pascal (MPa). Different materials,
 including bone, dentine, keratin, and wood were assigned to the
 corresponding parts. The materials of bone, dentine, and keratin were
 treated as isotropic linear elastic materials. For the detailed
 parameters, see Supplementary file 1 (<em>Drake et al., 2016; Huo et
 al., 2000</em>). However, twigs could not be treated as an isotropic
 linear elastic material, because the purpose of this simulation was to
 evaluate the food procuring efficiency of different taxa in different
 working conditions. Here the twigs were assigned as an orthotropic
 elastoplastic material, and the parameters (of wet red pine tree) were
 obtained from a wood handbook (<em>Risbrudt et al.,
 2010</em>). The occlusal surfaces were coupled by
 two arbitrary points on the upper and lower teeth. These two points were
 connected by a “beam connector”, which constrained all the degrees of
 freedom (<em>df</em>) between the two points (<em>Figure
 3–figure supplement</em><em> 2–4</em>). In this way, we
 simulated the occlusal surfaces of the upper and lower teeth. Jaw-closing
 muscles were simulated by several groups of “axis connectors”. This type
 of connector does not constrain any <em>df</em> of the two
 extreme points, and allows exerting force along the connector
 (<em>Figure 3–figure supplement</em><em>
 2–4</em>). Four jaw-closing muscles were considered, including
 <em>temporalis</em>, <em>superficial
 masseter</em>, <em>zygomaticomandibularis</em>, and
 <em>pterygoideus internus</em>. Ten axis connectors were
 assigned to the <em>temporalis</em>, and four, three, and
 three were assigned to the latter three, respectively. These connectors
 were uniformly arranged along their insertion areas based on their natural
 anatomy. The areas of the temporal fossa (<em>At</em>) and
 ascending ramus (<em>AA</em>) were measured in 3-matic
 Research (V12.0) (<em>Figure 3–figure
 supplement</em><em> 1</em>). We estimated the muscle
 force of <em>temporalis</em> as follows (<em>Tseng et
 al., 2017</em>): <em>At</em>
 (mm<sup>2</sup>) × 0.3 this force was
 equally distributed to the 10 axis connectors for temporalis.
 Alternatively, <em>AA</em>
 (mm<sup>2</sup>) × 0.3 was considered the
 gross force for <em>superficial masseter</em>,
 <em>zygomaticomandibularis</em>, and <em>pterygoideus
 internus</em>. This force was also equally distributed to the other
 10 axis connectors. Note that the <em>At
 </em>and <em>AA</em> in the
 <em>Platybelodon</em> and<em> Gomphotherium
 </em>models were not true. These models were scaled to the same
 volume as that of <em>Choerolophodon</em>. In the simulation,
 we uniformly assigned muscle forces to make it easy to compare models
 (<em>Supplementary file 1</em>). The
 cranium was treated as a rigid body and was fixed. Another boundary
 condition was assigned to a node of the mandibular condyle, which only
 allowed the y-direction rotation and constrained any other
 <em>df</em>s (simulating the rotation of the mandibular
 articulation). The <em>df</em>s of the x- and z- rotations of
 the mid-symphysis were also constrained by considering the connection to
 the other half. Two tests were carried out on the
 composite models: the distal forces test (<em>dft</em>) and
 twig-cutting test (<em>tct</em>). The <em>dft</em>
 includes two steps: 1, applying the muscle force; and 2, exerting a distal
 5,000 N force that gradually changes from horizontally to vertically, by
 which we assessed the optimum direction of the external force for the
 mandible of each taxon. In this test, the “twig” was not included in the
 model. In the <em>tct</em>, the middle
 point of the “twig” was set in close contact with the distal edge of the
 mandibular tusk and keratinous cutting plate, and was placed horizontally
 (<em>Figure 3–figure supplement 2</em>), 45° obliquely
 (<em>Figure 3–figure supplement 3</em>), and vertically
 (<em>Figure 3–figure supplement 4</em>). One extremity of the
 “twig” was fixed, and contact properties were assigned (hard contact
 normally and 0.3 frictional coefficient tangentially). The <em>tct
 </em>also includes two steps: 1, applying the muscle force as in
 <em>dft</em>; and 2, displacing the cranium and mandible 10 mm
 towards the “twig” to stimulate cutting action of proboscideans, by which
 we determined the cutting efficiency of directions for each taxon. In the
 results, the sum of the equivalent plastic strain (EPS) from total twig
 elements was calculated and reported for each model. The plastic strain
 represents the irreversible deformation of an element, and the sum of EPS
 from all twig elements can reflect the cutting effects in each model. The
 videos for von Mises stress contour colour maps were also generated in the
 <em>dft</em> (<em>Videos</em> 1–3) and
 <em>tct</em> (<em>Videos</em> 4, 6, 8, 10, 12, 14,
 16, 18) modelling, and for the EPS contour colour maps of twigs in
 <em>tct</em> (<em>Videos</em> 5, 7, 9, 11, 13, 15,
 17, 19) modelling.
description: # Title of Dataset: The trunk replaces the longer mandible as the main
 feeding organ in elephant evolutionSupplementary data Brief summary of
 dataset contents, contextualized in experimental procedures and results.
 ## Description of the Data and file structure **Codes information: (All
 codes are written in English)** Code S1  Script file for Bayesian
 total-evidence dating analysis (.nex file). (Open it in MrBayes Software)
 Code S2  Script auxiliary file for most parsimonious analysis, applying
 for some ordered and irreversible characters (characters 19–24 in the TNT
 program, please note that these character numbers were automatically
 assigned by TNT file, and begin from character “0” rather than “1”,
 therefore, these characters in supplementary Appendix and Data S1 were
 characters 20–25) (.tnt file). (Open it in TNT Software) Code S3  Script
 file for most parsimonious analysis (.tnt file). (Open it in TNT Software)
 Code S4  FEA setting file, for distal force tests (this file can be
 directly submitted in the Abaqus 6.14). Open it in Abaqus 6.14 Software;
 or open-source software Calculix
 ([http://www.calculix.de/](http://www.calculix.de/) and FEBio
 ([https://febio.org/](https://febio.org/) Code S5  FEA setting file, for
 distal twig cutting tests  (this file can be directly submitted in the
 Abaqus 6.14). Open it in Abaqus 6.14 Software; or open-source software
 Calculix ([http://www.calculix.de/](http://www.calculix.de/) and FEBio
 ([https://febio.org/](https://febio.org/) ## Sharing/access Information
 Links to other publicly accessible locations of the data: Was data derived
 from another source? If yes, list source(s): No