The trunk replaces the longer mandible as the main feeding organ in elephant evolution

Abstract

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. Our research provides insights into how unique environmental pressures shape the extreme evolution of organs, particularly in large mammals that developed various peculiar adaptations during the late Cenozoic global cooling trends.

Introduction

Proboscideans are known for their exceptionally elongated and versatile trunks (Shoshani, 1998). However, unlike modern elephants, proboscideans underwent a prolonged evolutionary phase characterized by the presence of greatly elongated mandibular symphysis and mandibular tusks. This elongation can be traced back to the Late Oligocene species Palaeomastodon and Phiomia, which are among the earliest elephantiforms (Andrews, 1906), and continued through to the Late Miocene Stegotetrabelodon, a stem elephantid (Shoshani, 1996; Tassy, 1996). Extreme longirostriny, a feature observed in fossil and modern fishes, reptiles, and birds, was relatively rare among terrestrial mammals and its occurrence in large-bodied proboscideans is particularly intriguing. Particularly, during the Early and Middle Miocene (approximately 20–11 Ma), the morphology of mandibular symphysis and tusks exhibited remarkable diversity, with over 20 genera from six families (Deinotheriidae, Mammutidae, Stegodontidae, ‘Gomphotheriidae’, Amebelodontidae, and Choerolophodontidae) displaying variations (Shoshani, 1996; Tassy, 1996). Why did proboscideans have evolved such a long mandible of so diversified morphology? How did fossil proboscideans use their strange mandibular symphysis and tusks, and what was the role of trunk in their feeding behavior? Finally, what was the environmental background for the co-evolution of their mandible and trunks, and why did proboscideans finally lose their long mandible? These important issuers on proboscideans evolution and adaptation remain poorly understood. Addressing these significant aspects of proboscidean evolution and adaptation requires comprehensive investigations into the functional and eco-morphology of longirostrine proboscideans.

Results

Mandible morphology of bunodont elephantiforms

During the Early and Middle Miocene, bunodont elephantiforms, as the ancestral group of living elephants, flourished (Osborn, 1936; Gheerbrant and Tassy, 2009; Cantalapiedra et al., 2021). Bunodont elephantiforms include Amebelodontidae, Choerolophodontidae, and ‘Gomphotheriidae’; all possess a greatly elongated mandibular symphysis (Gheerbrant and Tassy, 2009; Cantalapiedra et al., 2021). Our comprehensive phylogenetic reconstructions contained the majority of longirostrine elephantiform taxa at the species level and strongly supported this taxonomic scheme (Figure 1A; Figure 1—figure supplement 1, Figure 1—source data 1). The three families were characterized by their distinctive mandible and mandibular tusk morphology (Shoshani, 1996; Tassy, 1996). The paraphyletic ‘Gomphotheriidae’ have clubbed lower tusks (Figure 1B and K; Figure 1—figure supplement 2B); their mandibular symphysis is relatively narrow. This morphology is rather unspecialized, and the extant elephants are derived from ‘Gomphotheriidae’ (Tassy, 1996; Cantalapiedra et al., 2021). The mandibular symphysis of Amebelodontidae is generally shovel-like and the mandibular tusks are usually flattened and wide. Platybelodon is the most specialized genus within this family; it possesses extremely flattened and widened mandibular tusks with a sharp distal cutting edge (Wang et al., 2013: Figure 1C, I, and J; Figure 1—figure supplement 2A). Choerolophodontidae is unique because it completely lacks mandibular tusks and it has a long trough-like mandibular symphysis (Figure 1D, L, and M, Figure 1—figure supplement 2C–F). A very deep slit is present on each side of the distal alveolar crest (distal mandibular trough edge: Figure 1L), which is presumably for holding a keratinous cutting plate (Figure 1M), similar to the slits for holding large claws in felids and some burrowing mammals (e.g. anteaters). The anterior mandibular foramen of Choerolophodon is extremely large and tube-like (Figure 1—figure supplement 2D), which indicates a very developed mental nerve and eponymous artery (Eales, 1926), for the nutrition of tissues (i.e. keratins) grown from symphyseal dermis. The three lineages exhibit different evolutionary states of food acquisition organs (PC1 scores from mandible and tusk characters; see Materials and methods: Figure 1—figure supplement 4; Figure 1—figure supplement 5B). Different food acquisition organs morphology strongly indicate different methods of food acquisition among the three gomphothere families.

Figure 1.

Morphology of the narial region and mandible of three gomphothere families compared with an extant elephant, and the elephantiformes phylogeny.

(A) Phylogenetic reconstruction of major longirostrine elephantiforms at the species level based on the Bayesian tip-dating method. The node support (the number at each node) is the posterior probability, and the bars represent chronologic ranges of each taxon. (B–D) Representative cranium and mandible specimens of the three gomphothere families, including IVPP V22780, cranium, and IVPP V22781, mandible, of Gomphotherium tassyi [B], ‘Gomphotheriidae’, from Heijiagou Fauna, Tongxin region (TX4); HMV 0930, cranium and associated mandible of Platybelodon grangeri [C], from Zengjia Fauna, Linxia Basin (LX3); and IVPP V23457, cranium and associated mandible of Choerolophodon chioticus [D], from Middle Miaoerling Fauna, Linxia Basin (LX2). (E–H) Narial morphology of bunodont elephantiforms and elephantids in dorsal view, including IVPP OV733, Elephas maximus [E], a living elephantid; HMV 0930, P. grangeri [F]; IVPP V22780, G. tassyi [G]; and IVPP V23457, C. chioticus [H]. Mandibular morphology of bunodont elephantiforms. (I–J) Mandibular symphysis and tusks of HMV 0930, P. grangeri, in dorsal [I] and distal [J] views. (K) Mandibular symphysis and tusks of IVPP V22781, G. tassyi, in dorsal view. (L) Mandibular symphysis of IVPP V25397, C. chioticus, showing the deep slits at both sides of the distal alveolar crests in dorsal view. (M) Reconstruction of keratinous cutting plates in the slits, in dorsolateral view. Anatomic abbreviations: ce, cutting edge of the distal mandibular tusk in Platybelodon; kcp, reconstructed keratinous cutting plates in Choerolophodon; nb, nasal process of nasal bone; mc, slit or groove for mesethmoid cartilage insertion (white in color); pf, perinasal fossa; ps, prenasal slope in Platybelodon; s, slit for holding kcp in Choerolophodon.

Figure 1—figure supplement 1.

Phylogenetic reconstruction of longirostrine bunodont elephantiforms and mammutids using the maximum parsimony method (tree length, 253; CI, 0.470; RI, 0.779).

The node support (the number at each node) was calculated by symmetric resampling with 0.33 change probability (1000 replicates).

Figure 1—figure supplement 2.

Mandible of longirostrine bunodont elephantiforms.

(A) P. grangeri, HMV 0939, from Zengjia Fauna, Linxia Basin (LX3). (B) G. tassyi, IVPP V22781, from Hejiagou Fauna, Tongxin region (TX4). (C) Choerolophodon connexus, IVPP V8567, from Halamagai Fauna, Junggar Basin (JG2). (D–F) C. chioticus, from Lower and Middle Miaoerling Faunas, Tongxin region (TX2), including IVPP V23457, male, lateral view [D], showing the large, tube-like anterior mental foramen, and dorsal view [E] and IVPP V23213 [F], female, dorsal view. Anatomic abbreviations: mf, mental foramen; s, slit for holding keratinous cutting plates in Choerolophodon.

Figure 1—figure supplement 3.

Narial region of longirostrine bunodont elephantiforms.

(A and B) Choerolophodon guangheensis showing the large nasal bone processes and insertion for mesethmoid cartilage, lacking perinasal fossae; IVPP V17685 [A], type specimen, male, from Dalanggou Fauna, Linxia Basin (LX1); IVPP V31268 [B], female, from Lower Miaoerling Fauna, Tongxin region (TX1). (C) C. chioticus, IVPP V23457, male, from Middle Miaoerling Fauna, Tongxin region (TX2), showing a large nasal bone process and mesethmoid cartilage insertion, with very small perinasal fossae. (D and E) G. tassyi showing a moderately developed nasal bone processes and mesethmoid cartilage insertion, and well-developed perinasal fossae; HMV 0028 [D], type specimen, from Zengjia fauna, Linxia (LX3); IVPP V22780 [E], from Hejiagou fauna, Tongxin region (TX4). (F) Aphanobelodon zhaoi, HMV 1880, type specimen, from Upper Miaoerling fauna, Tongxin (TX3), showing a small nasal bone process and mesethmoid cartilage insertion, well-developed perinasal fossae, and initially developed prenasal slope. (G–J) P. grangeri, showing the very small nasal bone process and mesethmoid cartilage insertion, well-developed perinasal fossae and prenasal slope, including HMV 1841 [G], HMV 1840 [H], HMV 0940 [I], and HMV 0939 [J], from Zengjia Fauna, Linxia Basin (LX3). (K) Platybelodon tongxinensis, from Lintong region showing the small nasal bone process and mesethmoid cartilage insertion, well-developed perinasal fossae, and initially developed prenasal slope. (L) Platybelodon tetralophus, AMNH 26460, from Tamuqin Fauna, Tunggur region (TG3), showing a small nasal bone process, and well-developed perinasal fossae and prenasal slope. Anatomic abbreviations: nb, nasal process of nasal bone; mc, slit or groove for mesethmoid cartilage insertion (white in color); pf, perinasal fossa; ps, prenasal slope.

Figure 1—figure supplement 4.

Evolutionary levels of narial region (A) and of characters in relation to horizontal cutting (B).

Value in A was PC1 of characters 54–57 (see Appendix, Figure 1—figure supplement 5C, and Figure 1—source data 1); and that in B was PC2 of characters 5, 9, 11, 72, and 77 (see Figure 1—figure supplement 5D). Note that the clade of Platybelodon (marked by an asterisk) shows high evolutionary levels and that of Choerolophodontidae (marked by a circle) shows low evolutionary levels in both character combines, strongly suggesting the co-evolution of narial morphology and horizontal cutting behavior.

Figure 1—figure supplement 5.

Character combines in relation to feeding behavior in trilophodont longirostrine bunodont elephantiforms.

(A) PC1–PC3 plan of character combine of tusks and mandibular symphysis; (B) PC1 scores of A that are mapped on the phylogenetic tree (mammutids and stem elephantimorphs excluded); (C) PC1–PC3 plan of character combine of narial region; (D) PC2–PC3 plan of character combine in relation to horizontal cutting behavior. In A, C, and D, character loads were shown as line segments. Notes: character load vectors were marked in purple, narial region; red, mandible; blue, upper tusk; and green, mandibular tusk, colors, respectively. Taxa abbreviation: Afki, Afrochoerodon kisumuensis; Amfr, Amebelodon fricki; Apzh, A. zhaoi; Arfi, Archoerobelodon filholi; Chch, C. chioticus; Chco, Choerolophodon corrugatus; Chcn, C. connexus; Chgu, C. guangheensis; Chng, Choerolophodon ngorora; Eumo, Eubelodon morrilli; Gnth, Gnathabelodon thorpei; Goan, Gomphotherium angustidens; Gobr, Gomphotherium browni; Goin, Gomphotherium inopinatum; Gost, Gomphotherium steinheimense; Gopr, Gomphotherium productum; Gosy, Gomphotherium sylvaticum; Gota, G. tassyi; Melu, Megabelodon lulli; Phse, Phiomia serridens; Plba, Platybelodon barnumbrowni; Plda, Platybelodon danovi; Plgr, P. grangeri; Plte, P. tetralophus; Plto, P. tongxinensis; Prch, Protanancus chinjiensis; Prbr, Protanancus brevirostris; Prma, Protanancus macinnesi; Prwi, Protanancus wimani; Rhtl, Rhynchotherium tlascalae; Seba, Serbelodon barbourensis; Telo, Tetralophodon longirostris; Tobr, Torynobelodon britti.

Evolutionary dynamics and ecological niches of three longirostrine bunodont elephantiforms families

The three longirostrine bunodont elephantiforms families were widely distributed in the Early–Middle Miocene in northern China, ranging from ~19 to 11.5 Ma (Figure 2, Figure 2—figure supplement 1) (most of the Shanwangian and Tunggurian stages of the Chinese Land Mammal Age) (Wu et al., 2018). However, the relative abundances of these three families were very different and varied through time. We investigated fossil bunodont elephantiforms from four regions: Linxia Basin (LX), Tongxin region (TX), Junggar Basin (JG), and Tunggur region (TG) (Figure 2, Figure 2—figure supplement 1, Figure 2—source data 1). We focused on different fossil assemblages in different ages. We evaluated the amounts of fossils of each elephantiform taxon from major fossil assemblages in the above regions from three museums (IVPP, HPM, and AMNH, for full names, please see the Materials and methods) and calculated their proportions among all bunodont elephantiforms fossils (Figure 2A; Figure 2—figure supplement 2A). During ~19–17 Ma, the quantities of both Choerolophodon (Choerolophodontidae) and Protanancus (Amebelodontidae) were larger than those of Gomphotherium (‘Gomphotheriidae’) despite the high diversity of Gomphotherium in species level. During ~17–15 Ma (the Mid–Miocene Climate Optimum [MMCO]), the relative fossil abundances of the three gomphothere families were similar. However, within Amebelodontidae, Platybelodon appeared and rapidly replaced the primitive Protanancus. After ~14.5 Ma (i.e. the beginning of the Mid–Miocene Climate Transition [MMCT]) (Westerhold et al., 2020), Choerolophodon suddenly experienced regional extinction, and Gomphotherium also gradually declined. Within Amebelodontidae, Aphanobelodon, which is morphologically similar to Platybelodon, only occurred in one fossil assemblage, and the Platybelodon population greatly expanded. After ~13 Ma, Gomphotherium were regionally extinct, and Platybelodon became the only genus among bunodont elephantiforms.

Figure 2.

Relative abundance, tooth enamel δ¹³C, and stratigraphic ranges of the three gomphothere families in the Shanwangian and Tunggurian stages (~20–11.5 Ma) of northern China.

(A) Relative abundance of the three gomphothere families, including Choerolophodontidae, only represented by Choerolophodon (pink); Amebelodontidae, represented by Protanancus, Aphanobelodon (dark blue), and Platybelodon (light blue); and ‘Gomphotheriidae’, only represented by Gomphotherium (green). Horizontal bars indicate the average ages of the fossil assemblages, which are shown in C–F. The ages were determined by paleomagnetism (Supplementary file 1). The red curve shows the global reference benthic foraminifer oxygen isotope curve, which represents the global temperature (after Westerhold et al., 2020). (B) Tooth enamel stable carbon isotopic compositions of various gomphothere taxa. Each circle represents the bulk enamel δ¹³C values of a single tooth. The data of LX4 and TG2,3 are from previous publications (Wang and Deng, 2005; Zhang et al., 2009). Solid and dashed lines represent the mean and maximum enamel δ¹³C values for C₃ diets that have been corrected for Miocene atmospheric CO₂δ¹³C (after Tipple et al., 2010). (C–F) Synthetic stratigraphic columns of typical fossil bearing regions of northern China during ~19–11.5 Ma, which incorporated different fossil assemblages with different ages, from the Linxia Basin [C], Tongxin region [D], Junggar Basin [E], and Tunggur region [F]. Vertical bars represent the temporal ranges of different gomphothere taxa. Abbreviations for gomphothere taxa: Apzh, A. zhaoi; Chch, C. chioticus; Chco, C. connexus; Chgu, C. guangheensis; Goan, G. cf. angustidens; Goco, Gomphotherium cooperi; Goin, G. inopinatum; Gost, G. steinheimense; Gota, G. tassyi; Plda, Platybelodon dangheensis; Plgr, P. grangeri; Plte, P. tetralophus; Plto, P. tongxinensis; Prbr, P. brevirostris; Prsp, Protanancus sp.; Prwi, P. wimani.

Figure 2—figure supplement 1.

Relative abundance and tooth enamel δ¹⁸O values of fossil elephantiforms from the Shanwangian and Tunggurian stages (~19–11.5 Ma) in northern China.

(A) Relative abundances of the four elephantiform families, including Choerolophodontidae, only represented by Choerolophodon (pink); Amebelodontidae, represented by Protanancus, Aphanobelodon (dark blue), and Platybelodon (light blue); ‘Gomphotheriidae’, only represented by Gomphotherium (green); and Mammutidae (red). Horizontal bars represent the average ages of the fossil assemblages. The ages were determined by palaeomagnetism (Supplementary file 1). The red curve shows the global benthic foraminifer oxygen isotope curve, which represents the global temperature (after Westerhold et al., 2020). (B) Tooth enamel stable oxygen isotopic compositions of various gomphothere taxa. Each circle represents a single data point. (C–I) Plots of gomphothere tooth enamel stable carbon and oxygen isotopes in different fossil assemblages, including Dalanggou [C] and Laogou [D], Linxia Basin; lower Miaoerling [E] and upper Miaoerling [F], Tongxin Region; Halamagai [G] and Kekemaideng [H], Junggar Basin; and Moergen+Tamuqin [I], Tunggur region. (J) Isotopic data statistics from all gomphothere samples. The circles represent the average values, and the error bars represent the standard deviations. Abbreviations for fossil assemblages: JG1, Top Suosuoquan Fauna; JG2, Halamagai Fauna; JG3, Kekemaideng Fauna; JG4, Dingshanyanchi I Fauna; LX1, Dalanggou Fauna; LX2, Shinanu Fauna; LX3, Zengjia Fauna; LX4, Laogou Fauna; TG1, Tairum Fauna; TG2, Moergen Fauna; TG3, Tamuqin Fauna; TX1, Lower Miaoerling (Miaoerling A/B) Fauna; TX2, Middle Miaoerling (Miaoerling C) Fauna; TX3, Upper Miaoerling (Miaoerling D/E) Fauna; TX4, Heijiagou Fauna.

Figure 2—figure supplement 2.

Geographic distribution of Platybelodon worldwide.

The red circles represent the four regions examined in the present article, the light purple circles show the other regions with Platybelodon in northern China, and the dark blue circles denote the Platybelodon localities outside of China. This geographic map was generated by GeoMapApp (V3.6.15) (https://www.geomapapp.org/).

We analyzed stable isotope values of the tooth enamel for each gomphothere taxon in different fossil assemblages (Figure 2—source data 2). Overall, the δ¹³C values indicate a relatively open environment that consisted of a diverse range of habitats (including grasslands, wooded grasslands, and forests) dominated by C₃ plants in northern China from ~19 to 11.5 Ma (Figure 2B, Figure 2—figure supplement 1). In LX1 (Figure 2B; Figure 2—source data 2), at approximately ~18.5 Ma (see Supplementary file 1), both δ¹³C and δ¹³O values of Choerolophodon and Protanancus are similar and have a wide range of overlap, and those of Gomphotherium are in the middle of the range.

This indicates that the niches of these three groups overlapped without obvious differentiation. In TX1,2 (~17.3 Ma) (Deng et al., 2019), during which Platybelodon first appeared, the δ¹³C values of Platybelodon are lower than those of Protanancus. The δ¹³C values of Choerolophodon and Gomphotherium, however, are all within the δ¹³C range of Amebelodontidae. In JG2 (~16 Ma), the δ¹³C value of Platybelodon shows a distinct positive shift, which indicates expansion into more open habitats (Figure 2B). However, Choerolophodon may have persisted in a relatively closed environment (Figure 2B). Moreover, the ecological niche of Gomphotherium appeared to be in between those of Platybelodon and Choerolophodon, and they potentially lived in the boundary area of open and closed habitats (Figure 2B). This isotopic niche pattern is also observed in JG3 (~14.3 Ma), with the rise of Platybelodon and decline of Choerolophodon and Gomphotherium. In TX3, LX4, and TG2,3 (after ~14 Ma), only Platybelodon were examined (as Gomphotherium specimens are too rare to sample); its ecological niche was similar to that previously occupied by Choerolophodon after the latter went extinct.

Specialized feeding behaviors of three gomphothere families

To reconstruct the feeding behaviors of the three gomphothere families, we performed finite element (FE) analyses on three models as representatives of each family: Choerolophodon, Gomphotherium, and Platybelodon (for model settings, see Figure 3—figure supplements 1–4; Supplementary files 1 and 2). We conducted two kinds of tests: the distal forces test and the twig-cutting test. In the distal forces test, after full muscle forces were exerted, a 5000 N vertical force is loaded on the distal end of each mandible (Figure 3A–C, Videos 1–3). The mandible strain energy curve (MSEC) of Platybelodon suddenly reaches a very high value, while the MSEC increases of Gomphotherium and Choerolophodon are far less than that of Platybelodon (Figure 3D). Then, keeping the magnitude of the external force unchanged, with change of the external force from vertical to horizontal direction, the MSEC of Platybelodon decreases even lower than that of Gomphotherium and was similar to that of Choerolophodon. This result indicates that the mechanical performance of the Platybelodon mandible is disadvantageous under distal vertical external forces, but greatly improved under horizontal external forces (Figure 3A–C, Videos 1–3).

Figure 3.

Finite element analyses of feeding behaviors among three longirostrine gomphothere families and reconstruction of their feeding ecology.

(A–C) von Mises stress color maps of Choerolophodon, Gomphotherium, and Platybelodon models, with the full muscle forces exerted, and an additional 5000 N external vertical force applied on the distal end of the mandibular symphysis. (D) Strain energy curves of the three mandibles under the following three steps: (1) muscle forces linearly exerted; (2) a 5000 N external vertical force suddenly applied on the distal end; and (3) the 5000 N external force gradually changed from vertical to horizontal. (E) Sum of equivalent plastic strain from total elements (SEPS) of twigs cut by mandible models in three different directions (i.e. twig horizontal, 45º oblique, and vertical). Larger SEPS values indicate higher efficiency of twig cutting. (F) Scenery reconstruction of feeding behaviors of the three longirostrine gomphothere families (by X Guo), represented by Choerolophodon (Choerolophodontidae), feeding in a closed forest, Gomphotherium (‘Gomphotheriidae’), feeding at the margin between the closed forest and open grassland, and Platybelodon (Amebelodontidae), feeding on open grassland. (G) Detailed 3D reconstruction of Platybelodon feeding by grasping the grass blades using their flexible trunk and cutting the grass blades using the distal edge of their mandibular tusks.

Figure 3—figure supplement 1.

Geometric skull models for finite element (FE) analysis of three representative bunodont elephantiforms, including Gomphotherium (A), Choerolophodon (B), and Platybelodon (C), generated by Materialise 3-matic Research (V12.0).

For comparative purposes, the three bony parts of mandibles (excluding the tusk or cutting plate) were scaled to the same volume as those of Choerolophodon, and the other parts (cranium and tusk or cutting plate) were also scaled to the bony mandible of each model. The green areas indicate the surfaces for attaching temporalis (At), and superficial masseter+zygomaticomandibularis+ pterygoideus internus (AA), of which the areas were calculated for muscle force estimation. Note that the muscle forces of Choerolophodon were used in FE analyses, and those of Gomphotherium and Platybelodon were assigned the same magnitude as the muscle forces of Choerolophodon for comparative purposes.

Figure 3—figure supplement 2.

Mechanical settings of Platybelodon horizontal twig-cutting modeling; the blue cylinder represents the twig model.

(A) Constraints and boundary conditions of the model. (B) Modeling of the temporalis muscle; the muscle forces were exerted on the red axes. (C) Modeling of the superficial masseter muscle. (D) Modeling of the zygomaticomandibularis muscle. (E) Modeling of the pterygoideus internus muscle.

Figure 3—figure supplement 3.

Mechanical settings of Gomphotherium vertical twig-cutting modeling.

(A) Constraints and boundary conditions of the model. (B) Modeling of the temporalis muscle. (C) Modeling of the superficial masseter muscle. (D) Modeling of the zygomaticomandibularis muscle. (E) Modeling of the pterygoideus internus muscle.

Figure 3—figure supplement 4.

Mechanical settings of Choerolophodon oblique twig-cutting modeling.

(A) Constraints and boundary conditions of the model. (B) Modeling of the temporalis muscle. (C) Modeling of the superficial masseter muscle. (D) Modeling of the zygomaticomandibularis muscle. (E) Modeling of the pterygoideus internus muscle.

Video 1.

Finite element (FE) modeling of Platybelodon distal force test, color map showing the von Mises stress.

Video 2.

Finite element (FE) modeling of Gomphotherium distal force test, color map showing the von Mises stress.

Video 3.

Finite element (FE) modeling of Choerolophodon distal force test, color map showing the von Mises stress.

In twig-cutting tests, a cylindrical twig model of orthotropic elastoplasticity was posed in three directions to the distal end of the mandibular tusks of Platybelodon and Gomphotherium, and to the keratinous cutting plate of Choerolophodon. The sum of the equivalent plastic strain (SEPS) from total twig elements was calculated (equivalent plastic strain represents the irreversible deformation of an element, and the sum from all twig elements can reflect the cutting effects) (Figure 3E), and the cutting videos are provided (Videos 4–9; Video 10; Video 11; Video 12; Video 13; Video 14; Video 15; Video 16; Video 17; Video 18 Video 19). When the twig was placed horizontally, the SEPS of the Choerolophodon model is the largest, which means that Choerolophodon has the highest twig-cutting efficiency, followed by Gomphotherium; while Platybelodon exhibits the much lower efficiency in cutting horizontal twigs than the other two models. When the twig was placed obliquely (45° orientation), the Platybelodon model still shows the smallest SEPS value, although it was nearly one order of magnitude higher than that of the horizontal twig. Finally, when the twig itself was in a vertical direction, the growth direction of the keratinous cutting plate determines that Choerolophodon cannot cut in this condition. The SEPS of the Gomphotherium model also decreases and shows lower cutting efficiency. In contrast, the SEPS of the Platybelodon model increases another order of magnitude, which is substantially larger than that of any other taxa in any cutting state. These data strongly indicate that Platybelodon mandible is specialized for cutting vertically growing plants. However, the Choerolophodon mandible is specialized for cutting horizontally or obliquely growing plants, this explains the absence of mandibular tusks, but they are likely not able to feed on vertically growing plants. The cutting effect of the Gomphotherium mandible is relatively even for all directions.

Video 4.

Finite element (FE) modeling of Platybelodon vertical twig-cutting test, the total model, color map showing the von Mises stress.

Video 5.

Finite element (FE) modeling of Platybelodon vertical twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 6.

Finite element (FE) modeling of Platybelodon oblique twig-cutting test, the total model, color map showing the von Mises stress.

Video 7.

Finite element (FE) modeling of Platybelodon oblique twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 8.

Finite element (FE) modeling of Platybelodon horizontal twig-cutting test, the total model, color map showing the von Mises stress.

Video 9.

Finite element (FE) modeling of Platybelodon horizontal twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 10.

Finite element (FE) modeling of Gomphotherium vertical twig-cutting test, the total model, color map showing the von Mises stress.

Video 11.

Finite element (FE) modeling of Gomphotherium vertical twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 12.

Finite element (FE) modeling of Gomphotherium oblique twig-cutting test, the total model, color map showing the von Mises stress.

Video 13.

Finite element (FE) modeling of Gomphotherium oblique twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 14.

Finite element (FE) modeling of Gomphotherium horizontal twig-cutting test, the total model, color map showing the von Mises stress.

Video 15.

Finite element (FE) modeling of Gomphotherium horizontal twig-cutting test, the twig model (facing to the surface that is in contact with the tusk), color map showing the equivalent plastic strain.

Video 16.

Finite element (FE) modeling of Choerolophodon oblique twig-cutting test, the total model, color map showing the von Mises stress.

Video 17.

Finite element (FE) modeling of Choerolophodon oblique twig-cutting test, the twig model (facing to the surface that is in contact with the cutting plate), color map showing the equivalent plastic strain.

Video 18.

Finite element (FE) modeling of Choerolophodon horizontal twig-cutting test, the total model, color map showing the von Mises stress.

Video 19.

Finite element (FE) modeling of Choerolophodon horizontal twig-cutting test, the twig model (facing to the surface that is in contact with the cutting plate), color map showing the equivalent plastic strain.

Co-evolution of narial morphology and characters of horizontal cutting among gomphothere families

The evolutionary level of the trunk can be completely inferred from the morphology of the narial region (Tassy, 1994). Here, we first showed the narial region of a living elephant (E. maximus, IVPP OV733) (Figure 1E), focusing on the following four morphological factors. (1) The dorsal border of the narial aperture is slightly caudal to the postorbital process. This part is for the attachment of maxilla-labialis, which is the key muscle for manipulating the entire trunk (Boas and Paulli, 1908; Eales, 1926). (2) The narial aperture is wide, showing a pair of deep and sub-circular perinasal fossae. This part related to the insertion of lateralis nasi, and functions in enlarging the nostril cavity to suck up water in the trunk (Boas and Paulli, 1908; Eales, 1926). (3) The nasal process of the nasal bone is moderately developed. (4) The insertion slit for the mesethmoid cartilage is narrow and small, and deeply concealed in the narial aperture. The latter two points might be related to trunk flexibility. The smaller the nasal bone process and mesethmoid cartilage, the more flexible the trunk. These points were carefully discussed by Tassy, 1994.

The three lineages with different mandibular morphology also exhibited different stages of trunk evolution, which can be inferred from the narial region morphology (Tassy, 1994). Among the three groups, Gomphotherium has similar narial morphology to living elephantids (Figure 1E and G, Figure 1—figure supplement 3D, E), with comparable size and morphology of the nasal bone process, insertion for mesethmoid cartilage, and perinasal fossae. Alternatively, the narial region of Choerolophodon shows a relatively primitive evolutionary stage (Tassy, 1994; Figure 1H, Figure 1—figure supplement 3A–C). It has a very wide and large nasal bone process for an elephantiform, a fairly wide and long groove for mesethmoid cartilage insertion, and a pair of somewhat incipient (or even absent) perinasal fossae. Amebelodontids usually possess common narial morphology similar to that of Gomphotherium (Tassy, 1994). However, Platybelodon has noteworthy differences in narial morphology (Wang and Li, 2022; Figure 1F, Figure 1—figure supplement 3G–L). Its nasal aperture is greatly enlarged, which results in a very broad area for attaching maxillo-labialis (Boas and Paulli, 1908; Eales, 1926), the ‘core’ muscle of the proboscis. The nasal bone process is very short and stout, with a slight dorsal bulge. The slit for mesethmoid cartilage insertion is very tiny. Beside the well-developed perinasal fossa, a vast inclined region is positioned rostral to the perinasal fossa, this is hereafter referred to as the prenasal slope (Wang and Li, 2022), and it potentially provides additional attachment for the nasialis.

The narial morphology of Platybelodon and the closely related genus Aphanobelodon (Wang et al., 2017a: Figure 1—figure supplement 3F), is unique in elephantiforms, and shows even more derived characters than living elephants. Principal components analysis (PCA) was respectively performed for the characters of the narial region and the food acquisition organs (mandible and tusks) to extract their synthetic characters. On the phylogenetic tree (Figure 1—figure supplement 4A), the Platybelodon lineage (the dark purple star) displays the most derived narial morphological combination (PC1 scores from characters in relation to the narial region), while Choerolophodontidae (the dark purple circle) show the least specialized narial morphology, which is close to the stem taxon Phiomia. Interestedly, the evolutionary level of the character combine in relation to horizontal cutting (in terms of PC2) is highly correlated with that of narial region (Figure 1—figure supplement 4B). In the Platybelodon clade, this character combine also shows high evolutionary level; while in the Choerolophodontidae clade, the evolutionary level of horizontal cutting is rather low, comparable with that of narial region. The result strongly suggests a highly co-evolution between narial region and horizontal cutting behavior in the trilophodont longirostrine bunodont elephantiforms.

Discussion

In several fossil and living terrestrial mammalian groups, including living Tapirus (Moyano and Giannini, 2017) and extinct Proboscidipparion (Ma et al., 2023), Astrapotherium (Kramarz et al., 2019), and Macrauchenia (Blanco et al., 2021), elongated noses have evolved as food procuring organs. However, none of these groups possess a long and dexterous trunk like elephants, as they have not lost their mandibular incisors. Only living elephants are capable of solely accomplishing food procurement with their trunks. Coincidentally, longirostrine proboscideans, characterized by their extremely elongated mandibular symphysis and tusks, are the only group where these highly developed organs may have co-evolve (Tassy, 1996). It is evident that the great elongation of the nose is well matched with the extreme elongation of the mandibular symphysis. However, different lineages exhibit different evolutionary strategies strongly influenced by their ecological adaptations.

Our focus is on the ancestral elephantids, longirostrine bunodont elephantiforms, which include three families (Amebelodontidae, Choerolophodontidae, and ‘Gomphotheriidae’). These groups flourished during the MMCO, a period of global warmth from 17 to 15 Ma (Zachos et al., 2001; Westerhold et al., 2020). With the climatic shift during the MMCT (after ~14.5 Ma), Choerolophodon sharply declined and experienced regional extinction in northern China. Gomphotherium also declined but persisted until around 13 Ma, while Platybelodon became predominant until the end of the Middle Miocene (Figure 2A, Figure 2—figure supplement 1A). Worldwide, Choerolophodon and Gomphotherium continued to flourish during the MMCT and survived to the early Tortonian in other regions (Lambert and Shoshani, 1998; Göhlich, 1999; Sanders et al., 2010). However, Platybelodon was mostly restricted to Central Asia, especially northern China, with only a few records found in other regions (Wang et al., 2013; Wang and Li, 2022: Figure 2—figure supplement 2). These uneven distributions were likely influenced by climatic and ecological factors, such as the relatively drier climate and more open ecosystems in northern China, and Central Asia, strongly affected by the elevation of the Tibetan Plateau (Zachos et al., 2001; Wu et al., 2022). Extensive studies, using evidence from various disciplines including geochemical and geomagnetic proxies, and fossil records, have been conducted to explore this issue (Guo et al., 2002; Miao et al., 2012; Tang and Ding, 2013; Wang et al., 2022).

As discussed, the three groups can be distinguished by different mandibular symphysis and tusk morphologies (Gheerbrant and Tassy, 2009), which indicate different feeding behaviors. Their distinct narial regions further reflect different evolutionary stages of their trunks (Tassy, 1994). Choerolophodon possesses a highly specialized mandibular symphysis (Figure 1—figure supplement 5A and B) for cutting horizontally growing plants and was confined to relatively close habitats. The low evolutionary level of the narial region suggests a relatively primitive or clumsy trunk in Choerolophodon. The feeding strategy of Gomphotherium was unspecialized and flexible, relying on the coordination of the enamel band of an upper tusk and the corresponding lower tusk. Previous research also suggested that G. steinheimense, the sister taxon of elephantoids (Figure 1A, Figure 1—figure supplement 1) from the Halamagai Fauna, fed on grasses (Wu et al., 2018). While FE analysis could not provide a clear suggestion for the trunk function of Gomphotherium, the narial evolutionary stage, which is close to living elephants, suggests a relatively flexible trunk in Gomphotherium, as it is phylogenetically closer to living elephants than the other two groups.

Enamel isotope results also support the idea that Platybelodon expanded its living habitats to more open environments, such as grasslands, more than any other bunodont elephantiforms, likely due to its distinct feeding strategy. Our FE analyses strongly indicate that the Platybelodon mandible is specifically suited for cutting vertically growing plants. In open environments, there are vertically growing plants, such as soft-stemmed herbs. Platybelodon did not survive to the Late Miocene in northern China, possibly due to a mass extinction caused by the global climatic shift known as the Tortonian Thermal Maximum (Westerhold et al., 2020). However, from another perspective, if the Platybelodon trunk functioned similarly to that of living elephants, such as pulling out herbs from the earth (McKay, 1973), the greatly enlarged mandibular symphysis and tusks would have become redundant.

The reduction of the mandibular symphysis and tusks in the Late Miocene occurred in every lineage of elephantiform, coinciding with the large-scale expansion of C₄ grasses in the middle and low latitudes (Cerling et al., 1997). This may reflect the functional evolution of trunk grasping and manipulation in all elephantiforms lineages. From around 8 to 5 Ma in Africa, some derived members of ‘Gomphotheriidae’ (e.g. Anancus) and stem taxa of Elephantidae (e.g. Stegotetrabelodon, Primelephas) showed a strong inclination toward grazing on C₄ grasses, even though their cheek tooth morphology remained relatively primitive (Lister, 2013). Thus, open environments might be a key factor in both trunk development and the evolution of modern elephants.

As we have discussed, mandibular elongation was a prerequisite for the extremely long trunk of proboscideans, and open-land grazing further promoted the evolution of trunks with complex manipulative functions. This may explain why tapirs never developed a trunk as dexterous as that of elephants, as tapirs never shifted their adaptation zones to open lands. Living in dense forests, foods were easily accessible and procured through the mouth. Furthermore, the living mammals of the Late Cenozoic, in various open areas, have undergone specific organ evolution, such as the elongated necks of giraffes (Wang et al., 2022), the extravagate saber-tooth evolution in carnivores (Jiangzuo et al., 2023), and the development of various strange cranial appendages in different ruminants (Janis, 1982). It is possible that the highly evolved trunks of elephants evolved somewhat accidentally, under the pressure of ecological changes from closed to open environments.

Conclusion

In this study, we have examined the functional and eco-morphology, as well as the feeding behaviors, of longirostrine bunodont elephantiforms. Our findings demonstrate that multiple eco-adaptations have contributed to the diverse mandibular morphology observed in proboscideans, while open-land grazing has driven the development of trunk coiling and grasping functions and ultimately led to the loss of the long mandible. Specifically, the longirostrine elephantiform, Platybelodon, represents the first known proboscidean to have evolved both grazing behavior and trunk coiling and grasping functions. We have arrived at this conclusion through three lines of evidence, including the palaeoecological reconstruction based on tooth enamel stable isotope data, the reconstruction of feeding behaviors through FE analyses, and the examination of mandibular and narial region morphology correlated with characteristics associated with horizontal cutting behavior. The coiling and grasping ability of the trunk in Platybelodon evolved within the ecological context of Central Asia, which experienced regional drying and the expansion of open ecosystems following the MMCT (Miao et al., 2012; Tang and Ding, 2013). As a result, Platybelodon outcompeted other longirostrine bunodont elephantiforms and flourished in the open environment of northern China until the end of the MMCT. This scenario sheds light on how proboscideans overcame an evolutionary bottleneck. Initially, the elongation of mandibular symphysis and tusks served as the primary feeding organs, with the trunk being used as an auxiliary tool. However, through some necessary modifications such as tactile specialization and water intake adaptations (Shoshani, 1998; Purkart et al., 2022), the feeding function gradually shifted entirely to the trunk, which offered advantages in terms of flexibility and lighter weight of the feeding organs. Consequently, elephantiforms rapidly reduced the length of their mandibular symphysis and tusks (Shoshani and Tassy, 2005; Gheerbrant and Tassy, 2009). Similar stories of open-land adaptation and the acquisition of iconic characteristics have been observed in various mega-mammalian lineages, highlighting the crucial role of open-land adaptation for successful survival in modern ecosystems (Janis, 1982).

Materials and methods

Materials

The materials examined in this work are from three longirostrine gomphothere families, i.e., Choerolophodontidae, Amebelodontidae, and ‘Gomphotheriidae’ (Gheerbrant and Tassy, 2009). These materials include complete crania, mandibles, and teeth of different species, including C. chioticus, C. connexus, and C. guangheensis (Choerolophodontidae); P. dangheensis, P. grangeri, P. tetralophus, P. tongxinensis, P. brevirostris, Protanancus sp., P. wimani, and A. zhaoi (Amebelodontidae); G. cf. angustidens, G. cooperi, G. inopinatum, G. steinheimense, G. tassyi (‘Gomphotheriidae’). All were housed in three museums: the Institute of Vertebrate Paleontology and Paleoanthropology (IVPP), American Museum of Natural History (AMNH), and Hezheng Paleozoological Museum (HPM; HMV is the specimen prefix). For the detailed specimen list, please see Figure 2—source data 1. These materials were discovered from four regions, the Linxia Basin, Tongxin region, Junggar Basin, and Tunggur region (Figure 2; Figure 2—figure supplement 2); these regions are fossil-rich, especially during the Shanwangian and Tunggurian stages (~20–11 Ma), with different fossil assemblages in different ages (Deng et al., 2013; Wang et al., 2016; Wang et al., 2022; Qiu et al., 2013). For the age information of each fossil assemblage (Wang et al., 2022; Qiu et al., 2013; Wang, 2021), please see Supplementary file 1.

Cladistic analysis

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 P. serridens, an Oligocene basal elephantiform, was selected as outgroup. Additionally, a basal elephantoid, T. longirostris, 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 Wang et al., 2017b. 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 (Gavryushkina et al., 2014; Ronquist et al., 2012; Zhang et al., 2016). The Lewis Mkv model (Lewis, 2001), with gamma rate variation across characters (Mkv+G) (Yang, 1994), was initially used; subsequently, the timetree was modeled by the fossilized birth death process (Heath et al., 2014; Stadler, 2010). 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 (Lepage et al., 2007), 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 larger than 200 for all parameters and average standard deviations of split frequencies smaller than 0.01 (Geyer, 1992). The BTD analysis was performed in MrBayes 3.2.7 (Ronquist et al., 2012).

MP reconstruction was performed by TNT1.1 (Goloboff et al., 2008). 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>j (i to j in the matrix), and assigning the value ‘j–i’ to the blanks where i<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 major consensus tree was then calibrated using the time PaleoPhy function in R package paleotree 3.3.25 (Bapst, 2012), and is shown in Figure 1—figure supplement 1.

Stable carbon and oxygen isotope analysis

In this study, 83 tooth enamel samples from three gomphothere families were collected from gomphothere specimens (Figure 2—source data 2) for stable carbon and oxygen isotope analysis. We also compiled previously published isotope data of Platybelodon from the Tunggur region and Laogou, Linxia Basin (Wang and Deng, 2005; Zhang et al., 2009). The fossil teeth used in this study are well preserved and showed no visible signs of alteration. Tooth enamel samples were obtained by cutting a small patch of enamel from a tooth or drilling along the entire length using a rotary drill; then, the samples were ground into fine powder. All samples (2–3 mg enamel powder) were then pre-treated with 5% sodium hypochlorite (NaOCl) overnight to remove any possible organic contaminants and then cleaned with distilled water. These samples were then treated with 1 mol acetic acid overnight to remove non-structural carbonates and subsequently cleaned with distilled water. The treated samples were then freeze-dried.

The dried enamel samples were reacted with 100% phosphoric acid (H₃PO₄) at 25°C for approximately 72 hr. Carbon and oxygen isotope data were measured at Florida State University using a Finnigan MAT Scientific Delta Plus XP stable isotope ratio mass spectrometer coupled with a Thermo Scientific GasBench II. The lab standards that we used include MERK, MBCC, ROY-CC, and PDA. Results are reported in the standard delta (δ) notation as δ¹³C and δ¹⁸O values in reference to the international carbonate standard VPDB (Vienna Pee Dee Belemnite). We reconstructed the diet δ¹³C values of proboscideans from enamel δ¹³C values using an enrichment factor (ε*) of 13‰ for non-ruminants, like proboscideans (Cerling and Harris, 1999; Passey et al., 2005). Detailed results and specimen information are shown in Figure 2—source data 2.

FE analysis

We investigated the feeding behaviors of the three longirostrine gomphothere families, i.e., Choerolophodontidae, Amebelodontidae, and ‘Gomphotheriidae’, using FE stimulation. Three species were selected to represent each family, C. chioticus, IVPP V23457 (cranium and mandible); P. grangeri, HMV 0930 (cranium and mandible); and G. tassyi, IVPP V22780 (cranium) and IVPP V22781 (mandible), respectively. Note that, in Gomphotherium, we were unable to get access to cranium and mandible that belonged to one individual. We used a handheld Artec Spider 3D 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 Choerolophodon, 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 Platybelodon and Gomphotherium, and the keratinous cutting plate in Choerolophodon). 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 x–z plane, and the x-positive direction was set along the rostral direction. For comparison between different taxa, the three mandibular models were scaled to the same volume as Choerolophodon (6,563,708.0146 mm³) (Figure 3—figure supplement 1). 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 (V6.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 (V6.14) (Figure 3—figure supplements 2–4), 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 millimeter (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 (Drake et al., 2016; Bo and Quanshui, 2000). 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 (Risbrudt et al., 2010).

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 (df) between the two points (Figure 3—figure supplements 2–4). 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 df of the two extreme points, and allows exerting force along the connector (Figure 3—figure supplements 2–4). Four jaw-closing muscles were considered, including temporalis, superficial masseter, zygomaticomandibularis, and pterygoideus internus. Ten axis connectors were assigned to the temporalis, 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 (At) and ascending ramus (AA) were measured in Materialise 3-matic Research (V12.0) (Figure 3—figure supplement 1). We estimated the muscle force of temporalis as follows Tseng et al., 2017:

At (mm²)×0.3

this force was equally distributed to the 10-axis connectors for temporalis.

Alternatively,

AA (mm²)×0.3

was considered the gross force for superficial masseter, zygomaticomandibularis, and pterygoideus internus. This force was also equally distributed to the other 10-axis connectors.

Note that the At and AA in the Platybelodon and Gomphotherium models were not true. These models were scaled to the same volume as that of Choerolophodon. In the simulation, we uniformly assigned muscle forces to make it easy to compare models (Supplementary file 1).

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 dfs (simulating the rotation of the mandibular articulation). The dfs 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 (dft) and twig-cutting test (tct). The dft includes two steps: (1) applying the muscle force; and (2) exerting a distal 5000 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 tct, 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 (Figure 3—figure supplement 2), 45° obliquely (Figure 3—figure supplement 3), and vertically (Figure 3—figure supplement 4). One extremity of the ‘twig’ was fixed, and contact properties were assigned (hard contact normally and 0.3 frictional coefficient tangentially). The tct also includes two steps: (1) applying the muscle force as in dft; and (2) displacing the cranium and mandible 10 mm toward the ‘twig’ to stimulate cutting action of proboscideans, by which we determined the cutting efficiency of directions for each taxon. In the results, the SEPS from total twig elements was calculated and reported for each model. The plastic strain represents the irreversible deformation of an element, and the SEPS from all twig elements can reflect the cutting effects in each model. The videos for von Mises stress contour color maps were also generated in the dft (Videos 1–3) and tct (Video 4, Video 6, Video 8, Video 10, Video 12, Video 14, Video 16, Video 18) modeling, and for the EPS contour color maps of twigs in tct (Video 5, Video 7, Video 9, Video 11, Video 13, Video 15, Video 17, Video 19) modeling.

PCA and PC scores mapping on the tree

In PCA, only taxa of the three gomphothere families—members of Choerolophodontidae, Amebelodontidae, and ‘Gomphotheriidae’, in addition to P. serridens—were retained; the mammutid taxa (Losodokodon, Eozygodon, and Zygolophodon) and the stem elephantimorphs (Eritreum and Gomphotherium annectens) were excluded from the analyses.

A character combine, including characters 1–14 of upper and lower tusks, as well as 71–73, 76, 77, and 80 of the mandibles (Supplementary appendix and Figure 1—source data 1), was generated using PCA, which represents the synthetic character states of food acquisition organs. Besides, another two character combines, four characters concerning the narial region (i.e. characters 54–57) and five characters (characters 5, 9, 11, 72, and 77) in relation to the horizontal cutting behavior, were also generated using PCA. PCA was performed in Past 4.04, in which the missing values were predicted using mean value imputation. The data on PC1 vs. PC3 or PC2 vs PC3 plans were plotted (Figure 1—figure supplement 5A, C, D).

The PC1 of the food acquisition organ combination and that of the narial region, as well as the PC2 of the character combine in relation to the horizontal cutting behavior were selected to represent the synthetic evolutionary state of each (Figure 1—figure supplement 4; Figure 1—figure supplement 5B). We use PC2 rather than PC1 in the last test because the characters in relation to the horizontal cutting behavior were not evolved in one way. For example, the character 9 indicates the width of the mandibular tusks, which includes three states: 0, wide; 1, very wide; 2, narrow. In this case, the larger value (2=narrow mandibular tusker) does not mean the stronger horizontal cutting effect. Finally, these PCs were respectively mapped on the BTD tree (mammutids and stem elephantimorphs removed), using the contMap function in R package phytools 0.7-90 (Bapst, 2012).