Main text

The Cretaceous-Paleogene (K-Pg) mass extinctions accelerated the formation of modern-day global biota. In particular, the more than 6,000 species of living placental mammals trace their origins to the diversification of major orders around or after the K-Pg boundary ^1–3. Different hypotheses (for example: long-fuse, explosive, and short-fuse) about the timing of their adaptive radiation have been variably supported by analyses of molecules, fossils, and model-based approaches ^1,4–6. Regardless of the diversification scenarios favored by the competing explanations, the Paleocene epoch (66-56 Myr ago) has been highlighted as the most critical time interval for establishing phylogenetic and rate parameters for the placental adaptive radiation ^5,7. However, substantial asymmetry exists in the quantity and quality of terrestrial fossil records that document the first 10 Myr of the ‘Age of Mammals’, with ∼80% of known terrestrial K-Pg boundary sections occurring in North America ⁸. Earliest Paleocene mammals are not known from Europe ^9,10, rendering detailed reconstructions of post-K-Pg recovery in that continent difficult. The post-K-Pg fossil assemblages from Asia have hitherto not been considered in analyses of post-K-Pg recovery dynamics ¹¹.

The North American fossil record suggests that post-K-Pg mammal taxonomic diversity recovery was relatively rapid at higher latitudes; most occurred within the first ∼0.3 to 1 Myr of the Paleocene epoch ^12–17. This initial placental mammal diversification was driven by archaic groups (i.e., stem lineages and those with no living relatives), followed by the first appearance of many modern orders during two peak hyperthermal events in the past 66 Myr, at the end-Paleocene and early Eocene climatic maxima, respectively ^18,19. Distinct from these taxonomic recovery patterns, high selectivity of mammalian ecomorphological extinction across the K-Pg boundary indicates a primary productivity filter at the start of the Cenozoic Era (66 Myr ago to the present), disproportionately impacting large herbivores and their predators ²⁰. Additionally, the first 10 Myr of mammalian brain evolution after the K-Pg was marked by niche partitioning along a size gradient; endocranial traits reflecting more complex sensory processing did not appear until the Eocene (56 Myr ago). This phenomenon has been termed the ‘brawn before brains’ hypothesis ²¹. A similar pattern in initial size-driven diversification followed by expansion of ecomorphological disparity is also observed in mammal jaws, indicating a general dynamic across phenotypic systems ²². Whether this mode of evolution, whereby size disparity increase precedes other ecometric traits, should be understood as a global phenomenon during the post-K-Pg placental radiation remains untested, in large part because no such analyses have centered on non-North American continental fossil mammal records.

We investigated dental topography-performance shifts and timing of ecomorphological diversification by developing and leveraging the largest dataset to date of Paleocene Asian mammal assemblages. Our analyses focused on three of the most fossiliferous and biogeographically isolated Paleocene sedimentary sequences in paleotropical Asia: The Nanxiong, Qianshan, and Chijiang Basins in present-day south China ^23–27 (Fig. S1). We generated a new phenotypic dataset of 200 Asian Paleocene mammal teeth using high-resolution microcomputed tomography and laser scanning, capturing 37 species endemic to low-latitude east Asia and which are brought to bear on K-Pg recovery dynamics for the first time (Data S1-S2). We used dentitions as ecological indicators ²⁸ and examined temporal shifts in tooth crown complexity, curvature, and height and their association with tooth deformation resistance using topographic and simulation analyses. Our dataset spans the Paleocene, the first 10 Myr of the Cenozoic, enabling us to test the hypothesis that dental topography and tooth puncturing and shearing performance linkages showed delayed niche expansion relative to body size disparity increase during this initial period of post-K-Pg mammal recovery in Asia.

Dental traits paralleled Paleocene global and regional environmental conditions

We found that dental topographic trait variability in Paleocene mammals in south China tracked global and regional climatic changes despite stasis in high-level taxonomic composition over the course of the first 10 Myr after the K-Pg transition (Fig. 1). Dental height and sharpness variability were low in the beginning and end of the time interval, with a spike in the middle Paleocene corresponding to a short-lived negative excursion in global temperature. In contrast, elevated low-level taxonomic turnover of genera and species within the Paleocene indicates cladogenetic shifts, rather than anagenetic adaptation, underpin this dental topographic evolution (Data S1)¹¹. These findings suggest that in addition to its impact on crown group mammals, the K-Pg extinctions and subsequent climatic fluctuations also filtered out archaic placental mammals with lower speciation rates in favor of rapidly speciating taxa ⁴.

Fig. 1.

By contrast, there is no support for significant shifts in dental topographic trait mean values from the early to the middle Paleocene for the majority of analytical iterations of clade and dental data partitions (Figs 1A, S4; Table 1; Data S7). However, in most analyses we observed a significant shift in at least one dental topographic metric from the middle to the late Paleocene (Table 1; Data S7). The larger-bodied Chinese endemic pantodont (an ungulate-like archaic placental group; “CEP”) mammals tend to increase dental complexity (OPCR, orientation patch count rotated) and curvature (DNE, Dirichlet normal energy) whereas smaller, non-pantodonts (Chinese arctostylopids and anagalids; rodent-like archaic placental groups) in the dataset exhibited no significant dental traits shifts from the middle to late Paleocene (Fig. S4; Data S7). The transformation by CEPs in complexity and curvature indices relates to the capacity of teeth to resist wear and coincides with temperature and aridity increase towards the end-Paleocene thermal maximum event (Figs. 1-2). Multiple geological and geochemical proxies suggest that paleoclimate in and around the Nanxiong Basin K-Pg section in south China reflects a latitudinally much broader global tropical zone during the Paleocene ²⁹ relative to present-day Earth, as well as rapid shifts between more versus less humid intervals during the first 10 Myr of the Cenozoic ³⁰. Results from feldspar-quartz (F:Q) ratios, clay mineral composition analysis, diffuse reflectance spectroscopy (DRS), stable carbon isotope analysis, and total organic carbon (TOC) analysis all support this regional paleoclimate profile ^30,31. Local climate reconstructions for the latest Cretaceous indicate relatively warm and dry intervals, followed by warm and humid climates during the earlier Paleocene, then a return to less humid but still warm conditions in the later Paleocene ³¹. Chinese endemic pantodont dentitions tracked shifts in this paleoenvironmental progression.

Fig. 2.

Table 1.

The overall increase in dental complexity and curvature also coincided with an increase in drought-tolerant flora in south China (Fig. 2) and specifically with paleoenvironmental reconstructions in the Nanxiong Basin ³². Palynological evidence suggests that in addition to a predominance of broad-leaved, deciduous plants mixed with a smaller percentage of conifers than in fossil localities further north ³³, south China also recorded drought resistant taxa such as the maidenhair fern Pterisisporites in the early Paleocene Shanghu formation ³⁴ and again in late Paleocene samples ³⁵. Expanded comparisons across time and space suggest fossil pollen taxa that are indicators of seasonal aridity increased from 20.3% of late Paleocene pollen samples to 34.3% of early Eocene pollen samples ³⁵. These shifts appeared to have stabilized by the end of the Eocene, when taxonomically modern floras became established in Asia³⁶. In addition to tracking Paleocene temperature trends (Fig. 1), Asian mammal dental topography (OPCR and DNE in particular) also mimicked trends in the marine realm, where global planktonic foraminifera records demonstrate a similar pattern in species richness curve, as well as with global δ¹³C ^37,38 and south Pacific CO₂ levels ³⁹. These broader associations offer strong indications that endemic Paleocene mammals in south China composed a dynamic assemblage evolving in step with regional and global environmental changes during the first 10 Myr after the K-Pg mass extinctions.

We detected a shift and expansion of mammal dental morphospace in the late Paleocene (Fig. 2). An overall shift towards increased dental topographic trait magnitudes in late Paleocene samples is driven mainly by CEPs (Table 1), despite the fact that they constitute a minority of the overall data (41% of teeth) as well as late Paleocene partition (21% of teeth). Additionally, dental metric disparity is significantly higher in the late Paleocene partition than in the two preceding time bins (Table 1). This pattern is driven both by increased dental curvature and complexity (DNE & OPCR) in larger-bodied CEPs and crown height (RFI) in smaller-bodied non-pantodonts, in addition to differences in disparity among the clades (Data S7; bootstrapped variance for DNE, OPCR, and RFI are at least twice as large in non-pantodonts compared to CEPs). This suggests that an expansion of dental disparity in the late Paleocene occurred across the size spectrum of endemic Asian mammals. Over the same time interval examined, body size disparity appeared to have reached maximum saturation early on, in the early Paleocene (Fig. 1B, S6A, Table 1, Data S5), indicating that substantial increases in the disparity of dental complexity, curvature, and height lagged behind maximum size disparity during the Paleocene. This marked dental evolutionary pattern of ‘brawn before bite’ in endemic Asian mammals mirrors the endocranial and jaw functional morphology patterns identified in their North American and European counterparts ^21,22, underscoring that an initial size-driven post-K-Pg recovery followed by ecomorphological radiation was a global phenomenon, even as regional tectonic events such as the initial collision of the Indian subcontinent with Asia and Deccan Traps volcanism influenced local mammal evolution ^40,41.

Topography-performance integration underlies mammal dental shifts

Within an overall pattern of increasing integration between dental topographic traits and bite performance traits across the Paleocene time intervals (Fig. 3; Table S1), topographic versus performance trait variability shifted markedly both from early to middle and from middle to late Paleocene time bins (Figs S3-S4;). Early Paleocene to middle Paleocene DTA-FEA intra– and inter-correlations remained stable (Fig. 3B), but DTA-FEA inter-correlations strengthened while intra-DTA correlations weakened in the Middle to Late Paleocene transition (Fig. 3C). This transition pattern is in part due to the divergence of shape-performance linkages in CEP versus non-pantodont mammals. Dental topography variability tracked bite performance variability in non-pantodonts through time, peaking in the middle Paleocene (Fig S4D). By contrast, both intra– and inter-partition comparisons of topography and performance trends in CEPs showed a complete decoupling to shifts through time (Fig. S4G). These findings indicate that form-function malleability, the coexistence of distinct topography-performance relationships in each time and taxon partition while overall integration between the two trait groups increases between time bins, was present throughout the Paleocene. Complex form-function linkages generally promote evolutionary redundancy and can enhance optimization of phenotypic traits when selective trade-offs are present ^42–44. The presence of functional redundancies underlies the high levels of dental topographic variability in Asian Paleocene mammals. Alternatively, varying degrees of independence between the two performance traits and dental topographic traits analyzed could allow the two aspects of the dentition to evolve in a decoupled manner (Fig. 3).

Fig. 3.

South Asia as a Paleogene ‘Garden of Eden’ for placental mammals?

Among the most consequential implications of accurately interpreting post-K-Pg mammalian recovery dynamics in Asia is the ability to reconstruct the evolutionary conditions during the biogeographic origination of modern placental orders. Current knowledge and an overwhelming majority of data about the first post-dinosaur mammal ecological communities center on North American localities ^12,45, and secondarily from European records ⁴⁶. Yet, the fossil records of those continents pinpoint a high turnover of mammal faunas at the Paleocene-Eocene boundary, driven by the Paleocene origination of major modern mammal clades elsewhere and subsequent dispersal of those lineages to the European and North American continents, respectively. At least five living orders (Primates, Rodentia, Lagomorpha, Perissodactyla, Eulipotyphla), representing over half of mammalian species richness today, trace their evolutionary origins to Asia ^18,47–49. There is also mounting evidence that other organisms, including fish and plant lineages, followed a similar biogeographic pattern ⁵⁰.

The abrupt appearance of early representatives of modern mammal lineages in North America has been formally operationalized by the ‘East of Eden’ model, wherein Asia is the ‘Garden of Eden’ and a biodiversity pump of living orders of placental mammals ⁵¹ and other organisms. Biogeographic modeling analysis of modern mammalian diversity also strongly identifies tropical south Asia, including the geographic regions occupied today by southern China, India, and southeast Asian countries, as the cradle of mammalian diversification since the Paleocene ⁵². As geographically proximate faunal assemblages to the epicenter of the ‘Garden of Eden’ hypothesis, the endemic Paleocene mammals of the Asian paleotropics analyzed in this study exhibited a high degree of dental topographic variability while strengthening topography-performance integration during a period of biogeographic isolation and regional and global climate changes. This finding favors the scenario whereby many living orders of placental mammals were borne out of phenotypically and functionally plastic ancestral assemblages, including those that lived in south China, during the Paleocene ⁵³.

We further infer from these findings that episodes of fluctuating global warming during the first 10 Myr of the Cenozoic promoted the evolution of ‘all-purpose’ mammalian dentitions rather than those with specialized functions (Fig. 3). This prolonged dental topographic variability tracks extended post-K-Pg floral recovery times ^35,54, and suggests a ‘brawn before brains’ mode of placental mammal evolution may in part represent a series of correlated evolutionary shifts in sync with the steadily increasing complexity of Paleocene primary producers ^20,36. That a global signal would be detected in the splendidly isolated south Asian Paleocene assemblage is unexpected, and indicates that effects of climate forcing on post-K-Pg mammal recovery may be ubiquitous, as has been predicted for Earth’s ongoing rapid environmental shifts ⁵⁵. In response, Paleocene mammal clades in south China reached maximum tooth size disparity early and maximum integration between dental topography and bite performance later, all the while maintaining high levels of variability in dental complexity and convexity (Fig. 1). These preconditions set the stage for the subsequent taxonomic turnover between archaic and crown lineages during the Paleogene to Neogene modernization of mammal communities ^19,24. As a primary interface against shifting food resources and the environment, the placental dentition was poised to play an outsized role in the explosive diversification following the K-Pg extinctions, with the masticatory complex having already completed the macroevolutionary transition to a mechanically stiff but phenotypically flexible jaw during the Mesozoic Era ^56,57. The end-Cretaceous extinctions and the climatic volatility that followed then set in motion the ecological release that allowed mammals to explore dental form-function to a far fuller extent in 66 Myr of the Cenozoic Era than during the preceding 150 Myr under the reign of dinosaurs. The transformational evolutionary sequence of ‘brawn before bite’ in tooth size, topographic, and performance evolution was likely central to this placental success story.

Method details

We focused sampling on three clades: Chinese endemic Pantodonta, Chinese Arctostylopida, and Anagalida. Additional data were collected on other clades (e.g., Chinese Tillodontia) opportunistically when well-preserved specimens were available (Data S1). These three main clades together represent >50% of the species found in Paleocene faunas across China ¹¹. They also individually represent the most diverse clades across the three Asian Land Mammal Ages of the Paleocene: Shanghuan, Nongshanian, and Gashatan. By contrast, these clades are reduced to <25% of species diversity in Eocene assemblages in China, finally disappearing altogether in the late Oligocene. Therefore, we take the three clades to be representative of mammal assemblage dynamics during the Paleocene time interval in China.

Historically, Paleocene fossil mammal faunas in China have been defined based on biostratigraphic criteria and supplemented by magnetostratigraphic correlations where available. Here we follow the assessment of ⁵⁹ in generally correlating the Shanghuan Asian Land Mammal Age (ALMA) with the Puercan and Torrejonian North American Land Mammal Ages (NALMAs), the Nongshanian ALMA with early to middle Tiffanian NALMAs, and the Gashantan ALMA with late Tiffanian and Clarkforkian NALMAs, respectively (see ¹¹ for an alternative interpretation). Additionally, the Shanghuan and Nongshanian boundary was interpreted by ²³ to be at the top of Chron C27N, which is dated at 60.920 Myr ago in the Geomagnetic Polarity Time Scale ⁶⁰ but indicated as closer to ∼62.3 Myr ago in ⁵⁹. Similarly, the Nongshanian-Gashatan boundary has been variably defined at 59.24 Myr ago ⁵⁹ or the base of Chron C26N, which is 57.911 Myr ago ⁶⁰. In contrast, the end of the Gashatan coincides with the Paleocene-Eocene boundary at 56 Myr ago and has been consistently defined as such⁵⁹. Given these existing uncertainties, we use the terms ‘early Paleocene’, ‘middle Paleocene’, and ‘late Paleocene’ to refer generally to the Shanghuan, Nongshanian, and Gashatan ALMAs, respectively, in this study.

We analyzed 200 individual teeth from 48 specimens, representing 37 species (Data S1-S2). The specimens represent 2 upper first premolars, 4 upper second premolars, 15 upper third premolars, 19 upper fourth premolars, 22 upper first molars, 23 upper second molars, 16 upper third molars, 3 lower first premolars, 5 lower second premolars, 13 lower third premolars, 17 lower fourth premolars, 20 lower first molars, 22 lower second molars, and 19 lower third molars. The assigned geologic age of the studied specimens spans the entire Paleocene, binned into early (n = 91 teeth), middle (n = 46 teeth), and late Paleocene (n = 63 teeth) data partitions. To maximize sample size while minimizing disturbance to more delicate specimens, a combination of original specimens and previously produced high-fidelity specimen casts were used in our sampling of dental crown morphology.

Tooth shape estimates using Dental Topographic Analysis

We digitized dental crown morphology via microCT either using a GE v|tome|x m 300/180 kV micro-computed-tomography system (GE Measurement & Control Solutions, Wuntsdorf, Germany), housed at the Institute of Vertebrate Paleontology and Paleoanthropology (IVPP), Chinese Academy of Sciences (CAS), or a GE Phoenix Nanotom M in the Functional Anatomy and Vertebrate Evolution Laboratory, University of California, Berkeley. Projection images (1,000 to 1,500 images depending on specimen size) were acquired with an isotropic voxel size of 10 to 40 μm, at an energy range of 120-150 kV, current of 100-150 μA. Additionally, surface 3D scans (generated using IVPP’s Artec Space Spider 3D scanner, with a precision of 0.05 mm and resolution of 0.1 mm) were used for larger specimens to efficiently obtain surface morphology data. The enamel caps (the crown portion of each tooth above the enamel-dentine junction) were extracted from the specimen models using Geomagic Wrap v2021; all crown surface holes (generated from scanning imperfections or digitally removed sediment) were patched, the models were remeshed with triangles that have aspect ratios (base:height or edge-edge) of <10 and decimated to ∼10,000 triangles. Finally, we standardized the spatial orientation of all tooth crowns before exporting them as .ply files for dental topographic analysis (DTA). Individual tooth crown (enamel cap) models have the Z axis normal to the occlusal plane, the X axis parallel to the mesial-distal axis, and the Y axis parallel to the labial-lingual axis.

Exported meshes were further vetted and prepared for DTA within the R programming environment. We used functions implemented in the molaR R package ⁶¹ for all of the steps described next. First, all .ply files are further cleaned using ‘molarR_Clean’ to remove floating points/vertices and any triangular faces with zero area. Then, the batch function ‘molarR_Batch’ was used to calculate four dental topographic metrics: DNE, OPCR, Slope, and RFI.

Dirichlet Normal Energy (DNE) is a measure of occlusal sharpness, using a quantification of surface energy in tooth crowns relative to a gently curving or flat mesh surface. Steep, high, and/or shearing cusps tend to produce higher DNE values in DTA, whereas bulbous cusps tend to produce lower DNE values ^62,63. DNE has been used to successfully distinguish folivores, omnivores, and frugivores in euarchontan (primates, colugos, tree shrews) mammals ^62,64.

Orientation Patch Count Rotated (OPCR) is a measure of tooth crown complexity, using quantification of distinct patches on the crown face that possess unique orientations. Teeth with larger number of cusps and crenulations tend to produce higher OPCR values, whereas teeth with fewer cusps and simpler ridges tend to produce lower OPCR values. This metric has been used effectively to distinguish convergently evolved carnivore versus herbivore species across multiple mammalian orders ⁶⁵. DeMers and Hunter ⁶⁶ suggest that OPCR does not currently provide a reliable approximation for the degree of adaptation to herbivory, and that any interpretation of ecomorphological differences should be made within a taxon-specific context. Given the narrow set of clades targeted in our dataset, we interpret increasing OPCR as an indication of increased functional capability to process vegetation. A related measure, slope, quantifies the mean tilt of the tooth surface relative to the occlusal plane. Teeth with lower and more gently curving cusps tend to have lower slopes, whereas teeth with sharp and/or abruptly ridged crests tend to have higher slopes.

Relief Index (RFI) is a measure of the relative height and complexity of the tooth crown, using the values of the 3D surface and 2D ‘footprint’ areas of the tooth crown. RFI has been shown to be effective in distinguishing frugivores from insect and leaf specialists in euarchontan mammals ⁶². Furthermore, RFI is also able to distinguish carnivoran species at different trophic levels and dietary breadths ⁶³.

All tooth crown models are provided as ply files and deposited in FigShare (10.6084/m9.figshare.28611854).

Tooth performance estimates using Finite Element Analysis

The common inference that dental morphology can reflect dietary adaptations makes an underlying connection between the two traits through biomechanical performance; the mechanical performance imbued by a particular morphological configuration accomplishes one or more food acquisition and/or processing tasks. Such mechanical performance links can be tested using experimental and simulation approaches on both theoretical grounds and actual tooth shapes ⁶⁷. Here we applied a simulation approach to estimate two general performance traits commonly interpreted for mammalian tribosphenic teeth: puncturing/compressing and shearing/grinding ⁶⁸. The ancestral therian tribosphenic stroke, or the chewing movement underlying mammals with unfused mandibular symphysis and tribosphenic molars, has been reconstructed to involve significant components of (1) long-axis rotation in each hemimandible and (2) mortar-and-pestle grinding between upper and lower cusps ⁶⁹. These movements are associated with the ‘lock-and-key’ occlusion of a stereotypical mammalian tribosphenic dentition, which exhibit tall crest-like cusps in the trigonid and lower basin-like valleys in the talonid of lower molars that articulate with their corresponding upper teeth (Fig. S3).

We used finite element analysis (FEA) to estimate the work efficiency of different tooth crowns subjected to compressive/puncturing and shearing/grinding forces, the two major masticatory actions of tribosphenic teeth ^68,70,71. Adequate compressive and puncturing forces applied through individual or combined action of cusps are needed to overcome the fracture strength of brittle and tough food such as seeds and invertebrate exoskeletons in order for individuals to access soft and/or more nutritious tissues within. Sufficient shearing and grinding forces applied through surfaces or cusp edges are necessary to further break down food boli into smaller pieces to aid in the digestive process. The combination of these two major functions of tribosphenic teeth is thought to be a foundational adaptation that enabled the mammalian radiation ^68,69,72. We focused on the dental enamel portion of the tooth crown in our biting simulations because it is the most mineralized vertebrate tissue and typically the best-preserved component of fossil teeth. The fracture mechanics of mammalian dental enamel is thought to be intimately related to dietary function ⁷³, including the evolution of different enamel thickness in different feeding ecomorphologies ²⁸. It is important to note that because the dental data available for this study included not only CT scans, but also surface scans and specimen casts, enamel thickness was not quantifiable for most specimens in our dataset (see below for a description of enamel thickness scaling procedure applied in the simulations).

The two masticatory scenarios tested, compressing/puncturing versus shearing/grinding, have previously been shown to be broadly connected to tooth topology. Taller cusps (approximated in our study by Relief Index and partially by Slope) and more convex crowns (approximated in our study by DNE) tend to exhibit higher strain in finite element models of hypothetical tooth shapes ⁷⁴. However, there is by no means a consensus on the presence of a tight form-function linkage across all tooth types and/or clades. Whereas a dental topography to food mechanical property linkage can be detected in some extant carnivores ^75,76, bats ⁷⁷, and primate-like tooth models ⁷⁸, a form-function relationship between dental topography and tooth/food mechanical properties is absent in other model systems ^64,79. Therefore, we incorporated biomechanical simulations as an explicit test of the mechanical performance associations implied by the DTA-flora correlations established in the dental topography and paleoenvironmental portion of our analyses. The main question addressed by these simulations is not necessarily whether a form-function linkage exists in the Paleocene mammal dataset, but whether there is a consistent relationship between DTA and FEA traits across the three time bins of the Paleocene. If such consistent linkages are observed, it would suggest a stasis in dental topography and compressive/shearing performance through the first 10 million years of the Cenozoic in our dataset; alternatively, any significant changes in the topography-performance relationship would indicate substantial evolutionary change through the same time period.

We imported the same tooth crown models generated for DTA into the Strand7 finite element software (Strand7 Pty Ltd, Australia) to create finite element models of the individual teeth. The tooth meshes are composed of two-dimensional, three-noded triangular elements ranging between 9,000-10,000 elements, with assigned thickness parameter (see below) so they behave similarly to single-layer 3-D elements. A convergence test was done on two of the specimens (IVPP V5228, V5231) and indicated that the <10,000 three-noded triangular plate elements (tri3) returned simulation output values with less than 10% deviation from higher resolution meshes of the same models (Fig. S3). Therefore, we deemed the default resolution of the tooth meshes used in DTA to be sufficient for finite element simulations. Furthermore, because the dataset contained digital models built from both CT and surface scan data, not all specimens have enamel thickness information. Although increased enamel thickness has been associated with durophagy in mammals, there appears to be taxon-specific patterns of enamel thickness to tooth size and dietary ecology, in addition to extensive variation in intra-tooth enamel thickness distributions ^80–82. In absence of information on enamel thickness differences among the taxa studied, we standardized the thickness of the tooth crown models isometrically to be 10% of the total surface area of a given tooth model. As such, the simulation results should be treated as enamel cap shape-derived performance traits, rather than those based on a fully parameterized enamel model incorporating localized thickness and inter-specific allometric differences. We scaled applied forces on the teeth, simulating compressive or shearing loads, to be a value equivalent to total model surface area (but in units of Newtons instead of mm²). Compressive loads were represented on each tooth by a single nodal force vector directed towards the base of the crown along the height axis of each tooth, on the tallest cusp or cuspid (which is typically the protocone for upper molars and protoconid for lower molars, or the main cusp in premolars). This configuration simulated a compressive load directed into the tallest cusp (towards the root) on a given tooth from a food item (Fig S3), as reconstructed for the mammalian tribisphenic bite during a crushing or puncturing movement ^68,69. Shearing loads were represented on each tooth by two nodal force vectors, one on each of the two tallest cusps or cuspids; the two nodal forces are of equal magnitude but opposing directions, perpendicular to the long axis of the tooth and in the occlusal plane. This configuration portrays a shearing motion along the short axis of the tooth (Fig. S3), and simulates tooth-to-tooth and tooth-to-food contact during the ‘mortar and pestle’ grinding rotation movement reconstructed for the ancestral mammal tribosphenic bite ⁶⁹. This force magnitude standardization procedure effectively generates identical force to area ratios across the models, ensuring that the magnitudes of mechanical stress placed on the model are the same across models of different sizes; furthermore, we adjusted the strain energy values collected for each tooth model simulation using the volume and input force ratio-based correction equation provided in ⁸³. Given our objective to assess the association between relative performance traits and DTA metrics, with the latter being size-free variables, only the relative magnitudes of bite performance, rather than absolute values, were collected from the performance analyses.

After each tooth crown model was defined using the criteria outlined above, homogenous material properties were assigned. All models were defined as plate models (with thickness scaled to surface area) with a Young’s (elastic) modulus of 80 GPa ⁸⁴ and Poisson Ratio of 0.3. Next, each tooth was constrained from translation and rotation using four nodal constraints distributed in the four corners of each tooth, respectively (Fig. S3).

All bite scenarios were solved using the linear static solver implemented in Strand7. A total number of 400 analyses were performed (one compressive bite and one shearing bite simulation for each of the 200 tooth specimens in the dataset). We extracted two traits from the models: compressive bite strain energy and shear bite strain energy. Strain energy is defined as the area under a stress-strain curve of an object under load; operationally it measures the amount of work done by the load to deform the object (analogous to experimental work-to-fracture measurements)⁸⁵. An object that is more resistant to deformation would have lower strain energy than an object that easily deforms under load. The use of work-to-fracture measures to assess tooth performance is consistent with the understanding of mammalian tooth enamel as a fracture-prone, yet fracture-resistant, biological tissue ⁸⁶. In such a framework, fracture resistance is expected be directly related to the amount of force an individual is able to exert through the tooth crown during mastication, and by extension the hardest or toughest food item that can be processed without catastrophic damage to the tooth crown itself. Fracture resistance (as approximated by strain energy in this study) is also likely to be a stronger target of selection in mammals compared to those of other toothed vertebrates, given the former’s diphyodonty (having only two sets of tooth instead of continuous dental replacement) and thus necessity to prolong usable tooth lifespan⁸⁰. Therefore, we used strain energy values under compressive and shear bite simulation scenarios as a proxy for the effectiveness of each tooth crown model at resisting deformation from the respective biting forces. A tooth that is more resistant (i.e., have lower strain energy values) from compressive forces would be able to more effectively crush/puncture harder food items, and a tooth more resistant from shear forces would be able to cut/grind tougher food items. We modify Irschick et al.’s⁸⁷ definition of whole-organism performance and define bite performance in the context of this study as the capacity of individual tooth models to resist simulated compressive and shearing forces, which represent ecologically relevant factors that influence masticatory efficacy.

All finite element models are provided as NASTRAN/text files and deposited in FigShare (10.6084/m9.figshare.28611854).

Sensitivity validation of original versus cast specimen models

The tooth dataset developed in this study contains a mixture of CT scans of original and casts of fossil specimens. In order to assess the potential discrepancies in DTA and FEA trait values collected from cast versus original specimen derived models, either due to phenotypic details not captured in cast specimens or deterioration of epoxy or resin-based casts over time, we randomly sampled two specimens for which both original and cast based models were analyzed. Results from DTA and FEA of the CT-derived models of lower first molars from IVPP V5228 (Altilambda pactus) and V5231 (A. tenuis) were compared to those conducted on digital models built from CT scans of casts of those specimens.

We used the ‘n-point’ and ‘global’ alignment functions in Geomagic Wrap to first align the models using arbitrary homologous landmarks visually identified on both tooth models and then the automatic global alignment function to maximize overlap between the two models in 3D coordinate space. The ‘deviation’ function in Geomagic Wrap was then used to calculate summary statistics for linear deviations orthogonal to the surface of the original model. The V5228 Lm1 model showed a maximum deviation of 0.58 mm (5.8% of crown height) and average deviation of 0.04-0.09 mm (<1% of crown height), with a standard deviation of 0.11 mm (1% of crown height). The V5231 Lm1 model showed a maximum deviation of 0.38 mm (6% of crown height) and average of 0.03-0.04 mm (<1% of crown height), with a standard deviation of 0.06 mm (1% of crown height)(Fig. S3).

Quantification of uncertainty ranges for bootstrap sampling

We then subjected the cast-based models to DTA and FEA. In both cases, the DTA values between cast and specimen models are the most different for DNE (30-40% difference), followed by OPCR (18-32% difference), then RFI (2-13% difference), and finally Slope (1-5% difference) (Data S4). Also in both cases, adjusted compressive SE values differed by 14-41% and adjusted shear values differed by 20% (Data S3). Based on these specimen-cast validation tests, we set a conservative +/-40% uncertainty range for all DTA and FEA values obtained from all cast-based tooth models. DTA and FEA trait mean and variance estimates used for downstream statistical analyses were then calculated using a bootstrap sampling scheme where 1,000 replicates of the DTA and FEA datasets were passed through the statistical tests (see Quantification and statistical analysis section, below). Overall, the early and middle Paleocene time bins contain around 50% cast data, whereas the late Paleocene time bin contains 66% cast data. If data uncertainty from original versus specimen models present a major signal, we would expect the late Paleocene time interval to always show higher variance given an abundance of cast data in that time bin. We did not observe any consistent trends of high variance in the late Paleocene in the 1,000 bootstrap samples, suggesting that the results are not significantly biased by data quality differences between cast and original specimen derived models (Fig. S4).

Time bin duration correction and tooth position ratios

According to Ni et al. [6], the Shanghuan Asian Land Mammal Age (ALMA) faunas in Guangdong and Anhui range from 66-61.6 Ma (or 4.4 Myr in duration), the Nongshanian ALMA faunas in Guangdong, Jiangxi, and Anhui range from 61.6-59.2 Ma (or 2.4 Myr in duration), and the Gashatan ALMA faunas in Anhui rangesfrom 59.2 to 56 Ma (or 3.2 Myr in duration). Given the different durations represented in each of the three time bins used in our analyses, we corrected all variance estimates by dividing the variance calculated in each of the bootstrapped samples by the time duration of the respective time bin. In this manner, the variance values reported in the study represent per-million-year values.

We used pie chart analysis to verify that the proportions of tooth positions in each time bin are not substantially different from each other (Fig. S1). However, we caution that time bin comparisons of individual tooth positions are unlikely to be statistically robust because of small sample sizes. We only use aggregates of dental positions (all teeth, molars, or premolars) in our data partition analyses and interpretation.

Quantification and statistical analysis

We used the following features of the dental topography and performance data as the basis for assessing faunal assemblage trait shifts through the Paleocene: trait mean, variance, trait-to-trait correlation, and partition-to-partition correlation. We used ANOVA (Analysis of Variance) and pairwise t test to compare trait means, F test to compare trait variance, linear regression analysis and Kendall’s τ to compare trait-to-trait correlation, and two-block partial least squares (2B-PLS) analysis to compare partition-to-partition correlation. There are a total of six traits forming two main data partitions: topographic partition (DNE, OPCR, Slope, RFI) and performance partition (Compressive SE, and Shear SE). To assess the sensitivity of the results to subsets of the data, all of the trait comparisons were done iteratively using the total dataset, by time bin (early, middle, late Paleocene), taxon (Chinese endemic pantodonts versus non-pantodonts), and/or by dental position (all teeth, molar teeth, premolar teeth). All statistical tests were performed on 1,000 bootstrap resampled datasets that were parameterized using results of sensitivity tests on original versus cast specimen models, finite element mesh convergence tests, and time bin duration comparisons (see previous section). In addition, we quantified tooth size using two measures: total 2D surface area and square root of surface area. The R script for all statistical analyses and plots are included in Data S9.

Statistical resampling and tests

We used bootstrap resampled variance of dental topographic metrics and dental performance variables as a measure of tooth form-function disparity. 1,000 replicates of the 200-specimen DTA-FEA trait datasets were generated by sampling each DTA and FEA trait value from a uniform distribution defined using a +/-40% range of the original model-derived values. The bootstrap resampling was done with replacement, so it was possible for a given DTA/FEA trait value in a sampled tooth model to be repeated in a different replicate, but very unlikely for entire datasets to share similar values with another replicate. We used this resampling scheme to account for the uncertainty in trait estimates introduced by the totality of modeling, specimen preservation, cast versus original, and other potential sources of uncertainty in trait values. For each replicate, variance was calculated for each metric, time interval, for Chinese endemic pantodont (CEP) versus non-pantodont, and molar versus premolar data partitions. To assess shifts through time, statistical differences between the variance of each dental topographic metric sample in adjacent time intervals (A, early Paleocene; B, middle Paleocene; C, late Paleocene) were evaluated using F tests. A null hypothesis of a variance ratio of 1 between pairwise comparisons was tested using the var.test() function in R ⁸⁸. The outputs of the variance tests on a give bootstrap sample were then compiled for all 1,000 replicates, and the overall variance mean and variance test output used as the statistical basis for detecting significant differences between time, taxon, and tooth data partitions. Unless indicated otherwise, all tests described below were also done using the 1,000 bootstrap samples.

We evaluated differences in central tendency (mean value) of dental topographic metrics and dental performance variables through time using analysis of variance (ANOVA) and pairwise t tests. Each dental topographic metric was tested against the three time intervals defined above, using the aov() function in R. For statistically significant (at the p < 0.05 level) results, we additionally evaluated the pairwise intervals that contribute significant differences in mean dental topographic values. We assessed pairwise differences using pairwise t tests implemented with a Holm correction for multiple comparisons. The t test was conducted using pairwise.t.test() in R.

We constructed a tooth morphospace using all four dental topographic metrics analyzed by principal components (PC) analysis. The first and second PC axes were chosen to visualize the two-dimensional morphospace. Additionally, we quantified the degree of morphological disparity and statistical differences in disparity between adjacent time interval data partitions. All PC scores generated from the PCA were included in the disparity analysis. The prcomp() function in R was used for PCA, and the morphol.disparity() function implemented in the Geomorph R package ⁸⁹ was used for morphological disparity significance tests. Input data for the PCA were from the original dataset values, not from the bootstrap samples.

Linear regression analysis between individual dental topographic metric (DNE, OPCR, Slope, RFI) and dental performance metric (compressive bite strain energy, shear bite strain energy) was performed to quantify the correlation between dental form and function. Adjusted R² and p values generated from the lm() function in R were used to evaluate the goodness of fit and statistical significance of the form-function relationships. The distribution of R² values and their corresponding p values are reported in Fig. S6.

In addition to pairwise form-function linear regression analyses, we also evaluated the degree of correlation between the DTA and FEA data blocks as a means to measure integration between dental topography and deformation resistance. We used two-block partial least squares analysis⁹⁰ coupled with bootstrap resampling to generate distributions of correlation coefficients (r-pls) for each of the three time bins, and then tested for significant differences in the magnitude of DTA-FEA correlation between adjacent time bins using Welch’s two-sample t tests on the correlation coefficient distributions. DTA-FEA correlation differences at p < 0.05 are interpreted as a significant shift in the degree of integration of the two traits from one time bin to the next time bin.

All R scripts used in the analyses are included as supplementary files (Data S8, S9).

Resource availability

• Requests for further information and resources should be directed to and will be fulfilled by the lead contact, Jack Tseng (zjt@berkeley.edu).

• Code associated with analyses is included with this article.

• Raw data are available in the supplementary materials.

• All original fossil specimens are accessioned in the Institute of Vertebrate Paleontology and Paleoanthropology, Chinese Academy of Sciences.

Acknowledgements

We thank I. Ruiz for assisting with generating 3D tooth models from CT data; X. Tseng for assisting with data collection; P. Holroyd for assisting with locating specimen casts in the collections of the University of California Museum of Paleontology; J. Liu for helpful discussion and comments on earlier versions of the study. Z. Luo provided earnest and constructive feedback on an earlier version of the study.

Additional information

Funding

This study was funded in part by grant 213109 from the Nanjing Institute of Geology and Paleontology, Chinese Academy of Sciences (to QL and ZJT) and National Key Research and Development Project of China (2024YFF0807603; to QL).

Author contributions

Conceptualization: All authors. Methodology: ZJT.

Formal analyses: ZJT.

Investigation: All authors. Results acquisition: ZJT.

Raw data acquisition (surface and CT scan data): ZJT, QL. Writing, original draft: ZJT.

Writing, review, and editing: All authors. Visualization: ZJT, QL.

Additional files

Document S1. Figures S1–S6, Tables S1 and S2, and supplemental references.

Data S1.

Data S2.

Data S3.

Data S4.

Data S5.

Data S6.

Data S7.

Data S8.

Data S9.

Data S10.