Two major ecological shifts shaped 60 million years of ungulate faunal evolution

Domain: earth_science
Published: 2025-06-05
Editorial outcome: minor_revision

Abstract
The fossil record provides direct evidence for the behavior of biological systems over millions of years, offering a vital source for studying how ecosystems evolved and responded to major environmental changes. Using network analysis on a dataset of over 3000 fossil species spanning the past 60 Myr, we find that ungulate continental assemblages exhibit prolonged ecological stability interrupted by irreversible reorganizations associated with abiotic events. During the early Cenozoic, continental assemblages are dominated by mid-sized browsers with low-crowned teeth, which show increasing functional diversity. Around 21 Ma, the formation of a land bridge between Eurasia and Africa triggers the first major global transition towards a new functional system featuring a prevalence of large browsers with mid- to high-crowned molars. Functional diversity continues to increase, peaking around 10 Ma. Shortly after, aridification and the spread of C4-dominated vegetation lead to a second tipping point towards a fauna characterized by grazers and browsers with high and low crowned teeth. A global decline in ungulate functional diversity begins 10 Ma ago and accelerates around 2.5 Ma, yet the functional structure of these faunas remains stable in the latest Cenozoic. Large mammal evolutionary history reflects two key transitions, aligning with major tectonic and climatic events.

Introduction
The ongoing biotic crisis is characterized by a rapid deterioration of most Earth’s ecosystems1, which could be irrevocably tipping several terrestrial and marine ecosystems towards new ecological regimes (i.e. functional tipping points2–4). Such shifts in functional states, which are abrupt and irreversible even over geological time, have been identified in the fossil record5,6, suggesting that past perturbations can be used as analogs to gain insights about the tempo and nature of immediate and future responses to environmental changes7–9. Faunal shifts in the fossil record have been traditionally tackled from a taxonomic perspective (i.e., changes in the taxonomic composition of communities and faunas10,11). More recently, paleobiological studies have progressively incorporated taxon-free functional approaches, utilizing changes in ecomorphological composition to evaluate ecological shifts beyond taxonomic turnover5,12–14.

In terrestrial ecosystems, most of these functional approaches have focused on patterns of collapse of mammalian assemblages during the late Quaternary12–14. A deeper temporal dimension of taxon–free approaches is needed to provide a more profound understanding of the evolution of mammalian ecosystem dynamics and how they responded to major environmental disruptions. Broadening the temporal scope should also reveal the mechanisms related to other patterns beyond collapse episodes, like those involved in modulating ecological resilience and stasis15, as well as processes underpinning dynamics of recovery and the unfolding of functional diversity.

Because environmental change has a huge impact on ecosystem primary productivity and biomass transfer16–18, investigating long-term patterns in primary consumer guilds should provide valuable insights into major ecological reconfigurations. Mammalian large herbivores have been the major vertebrate primary consumers in our planet’s ecosystems in Cenozoic times. Ungulates consume disproportionately more plant biomass per unit area than any other vertebrate group19. Many large mammalian herbivores are considered ecosystem engineers because they play a crucial role in processes involved in fire regimes control, seed dispersal, ecological succession of ecosystems, and soil compaction19–25. Crucially, because of the recent and ongoing decline of large herbivores, these functions are largely missing from modern ecosystems24. In a search for a broader perspective on the behavior of terrestrial functional systems, we ask: how did terrestrial functional faunas shift throughout the Cenozoic? Are functional tipping points a common feature in these systems? If so, which factors triggered such state transitions?

Here, we aim to identify the main shifts in global large herbivore functional structure of the orders Artiodactyla, Perissodactyla, and Proboscidea, which constitute most of the ungulate fossil record of the Cenozoic. As ecosystem engineers, terrestrial ungulates have a tight connection to physical changes and habitat disruptions, making them ideal to study current and past ecosystem functional transitions19,25. We developed an extensive database for large herbivores, containing occurrence information for 3012 species spanning the past 60 million years. For each species, we gathered information for 14 ecomorphological traits (13 dental traits and body size) and clustered species based on their pairwise functional distances. We refer to these groups as ‘functional types’ (FT), since their constituent species would fulfill comparable ecological roles. We then applied Network Analysis to identify ungulate continental functional faunas (UFFs)5, defined as unique associations of functional types. Notably, UFFs are not defined solely by the presence of FTs, but by their relevance, affinity, and fidelity in the faunas (see methods). Our network analysis ignores assemblage age (it focuses exclusively on the functional structure), and the temporal succession of UFFs is determined a posteriori (see methods). Lastly, we compared the results of our functional analyses with taxonomic faunas and the timing of major environmental shifts of the Cenozoic26–34. Our results unveil two major shifts in the UFF succession coinciding with abiotic events.

Results
Our results show a decoupling of taxonomic and functional assembly dynamics (Fig. 1A and Supplementary Fig. 1). The functional structure of large herbivore assemblages exhibits greater resilience compared to their taxonomic structure. On average, taxonomic reassembly (module change, see Supplementary Fig. 1) occurred every 1.7 ± 0.31 Myr, whereas functional reassembly took place every 6.25 ± 1.6 Myr. Sensitivity analyses validate these findings, ruling out potential artifacts stemming from sampling, aggregation or the selection of community detection algorithms (see methods). Large herbivore assemblages experienced long periods of ecological stasis, punctuated by two episodes of global and almost synchronous functional reassembly (tipping points, Fig. 1A). These tipping points occurred around 21 and 10 Ma, likely as a response to major abiotic events. Since 10 Ma, large herbivore assemblages have exhibited remarkable stability, first in Africa, followed by the Americas and Europe around 7 Ma, and later in Asia around 4.5 Ma (Fig. 1A).

Early Cenozoic stability
For most of the Cenozoic, the faunas of America, Europe, and Asia were grouped into a single module, representing a shared functional configuration (Fig. 1A), and occupied similar regions of functional space (Fig. 1B). Before 30 Ma, African faunas included only proboscideans among the orders examined, preventing a comprehensive assessment of the continent’s functional structure during this time period.

The Paleogene was a time dominated by mid-size browsers with low-crowned teeth (brachydont) (Supplementary Figs. S4 and S17). Between 60 and 40 Ma, the system experienced a long period of functional stability (ecological stasis), with a first reconfiguration in UFFs around 40 Ma, when American ungulate assemblages diverged functionally from those in Europe and Asia (Fig. 1A).

During the early Cenozoic, there is a global increase in functional diversity (Fig. 2). Functional diversity encapsulates both functional richness (the number of functional types) and functional evenness (the degree to which species are equally divided across functional types, see methods). Initially, this increase in functional diversity is primarily attributable to an increase in functional richness (Fig. 2), especially in American faunas, which exhibited a sustained increase in functional diversity between 50 and 20 Ma (Fig. 2B). Functional innovations were likely intertwined with biogeographic events, contributing to this enrichment. For example, around 37.5 Ma, a significant immigration of large mammals from Asia to North America occurred following a sea-level drop10. In contrast, ungulate faunas of Europe and particularly Asia underwent a prolonged functional diversity decline during the early Cenozoic (Fig. 2). At about 40 Ma, these diverging trends culminated in the split of American ungulate assemblages from those in Europe and Asia (Fig. 1A, B). This was primarily driven by the emergence of functional dominance of FTs in larger body size categories in the Americas, in contrast to the mid-sized herbivores that characterize the UFF of Eurasia (Supplementary Fig. 17). The geographical extent of landmasses could have significantly influenced the survival dynamics of different faunal groups. Specifically, the broader southern reach of the Eurasian landmass compared to North America may have provided a more stable environment for mid-sized herbivores. This stability persisted despite the overarching trend of Cenozoic cooling, which likely affected species differently depending on their latitudinal distribution.

Europe and Asia formed a single UFF and shared a common functional trajectory until around ca. 32 Ma, when the European UFF diverged from the shared module with Asia (Fig. 1A). This divergence coincided with a period of global cooling that ultimately led to the formation of permanent ice sheets in Antarctica, marking the transition from the greenhouse climate of the Early Paleogene to the onset of the Late Cenozoic icehouse climate. This climatic shift, known as the Eocene-Oligocene transition (EOC), triggered major faunal extinctions and turnover, primarily in Europe, in a faunal event also dubbed as the Grande Coupure26,35. Characterized by the influx of migrant faunas from Asia to Europe, the Grande Coupure replaced archaic ungulates with modern faunas35,36. The asymmetry in the replacement is explained by differences in tempo, biogeography and differential climatic regimes26,37. The EOC caused less faunal disruption in Asia, where archaic forms survived until the late Early Oligocene, retaining a faunal assemblage more similar to the ones of the Late Eocene (e.g., archaic Amynodontidae forms, such as Amynodon or Sianodon)37. Furthermore, modern forms gradually appeared in Asia prior to the EOC, including some families that had already emerged in the late Eocene (e.g., Entelodontidae), in contrast to the abrupt replacement observed in European faunas immediately after the EOC37. This abrupt faunal replacement drove a decline in the functional diversity of European faunas, primarily due to a reduction in functional richness (Fig. 2). This, combined with the abrupt loss of archaic ungulates and their replacement by modern forms, caused the European faunal assemblages to diverge from the late Eocene ones, which remained an important component of Asian faunas. The biogeographic and climatic context of the EOC transition also played a crucial role in shaping the observed functional shifts26. Europe’s restricted high-latitude position and its configuration as an archipelago during this period may have constrained faunal dispersal and climate tracking, making species more vulnerable to extinction as conditions changed. In contrast, larger and more continuous landmasses, such as those in North America and Asia, provided broader climatic gradients and refugia, which facilitated species persistence and functional stability (Fig. 1A). These geographic constraints, coupled with the rapid cooling and heightened seasonality of the EOC transition26, contribute to explaining the pronounced shifts observed in European faunal assemblages. This ultimately led to the divergence of European faunas in terms of functional structure (Fig. 1).

The African-Eurasian land bridge
The first global-scale reconfiguration of UFFs took place between 26 and 21 Ma, revealing a major first tipping point. This period was marked by a homogenization of the functional structure of European, Asian, and African faunas (Fig. 1). The functional coalescence is reflected in the drastic reduction of the functional distance between Africa and Eurasia, beginning around 25 Ma, during the Chattian stage (27.82–23.03 Ma) in Europe, and later followed by the Aquitanian (23.03–20.44 Ma) in Asia (Supplementary Fig. 3). A long process of environmental change, marked by a trend towards global aridification, began around 26 million years ago with the Alpino-Himalayan orogeny27–29. This caused the collision of the African and Eurasian plates, which by 21 Ma resulted in the complete closure of the Tethys Sea, the emergence of the Mediterranean Sea, and the formation of a land connection known as the Gomphotherium land bridge30.

The homogenization of faunas and the reassembly into a coherent, transcontinental UFF was therefore probably the result of the terrestrial connection between Eurasia and Africa at this time (Fig. 1A). This resulted in a net gain in functional diversity, with all three orders (Perissodactyla, Artiodactyla, Proboscidea) present across all three continents (Fig. 2). Notably, proboscideans dispersed out of Africa and spread around the world, triggering an acceleration of ecomorphological evolution within this order38. However, the impact of the Gomphotherium land bridge worked in both directions. Artiodactyls and perissodactyls, which had dominated the large herbivore fauna in Eurasia since the beginning of the Cenozoic, dispersed into Africa while exploring new regions of the functional space (Fig. 1B). Overall, the onset of the tectonic movements that lead to the formation of the Gomphotherium land bridge created new ecological configurations, marking a milestone event in the evolutionary history of terrestrial mammal faunas.

Simultaneously, the Americas underwent a reassembly of their functional assemblage, characterized by mid-size grazers with medium to high crowned teeth (Supplementary Figs. 4 and 5). However, the Americas continued to follow their distinct trajectory, retaining an exclusive state in terms of functional structure compared to the other continents (Fig. 1A). In contrast to the trend observed in Eurasia, around 15 Ma the American UFF began to experience a continuous decline in functional diversity (Fig. 2B). Interestingly, this decline follows an immigration event between Eurasia and North America at the end of the early Miocene, some 16.5 Ma10, suggesting that faunal dispersals do not necessarily enhance functional diversity but may instead reduce the functional richness of endemic American functional types. After 10 Ma, Africa and Eurasia would follow the Americas in this trend of functional diversity decline (Fig. 2B).

The global expansion of grasslands
Following a global sustained trend of increasing functional diversity fueled by a rise in functional richness (Fig. 2), the functional network reaches a second global tipping point ca. 10 Ma (Fig. 1A)32,33. From this moment onwards, we detect a drop in global functional diversity that has persisted to the present day (Fig. 2). This trend resulted in the loss of 65% of functional diversity in less than 5 Myr (from 10 to 5 Ma, see Fig. 2). The functional loss prompted an almost-contemporaneous global reassembly of UFFs characterized by a combination of grazers and browsers with selenodont teeth, represented by ruminants in Africa and Eurasia, and by perissodactyls with bilophodont teeth in the Americas (Supplementary Figs. 4, 5). This reassembly was first observed in African faunas, resulting in their divergence from the previous module shared with Eurasian faunas (Fig. 1). European faunas, followed by the faunas of the Americas and Asia, also underwent a reorganization of their functional structure similar to Africa’s (Fig. 1A). This reconfiguration made European and Asian functional assemblages to shortly diverge from each other for two million years (Fig. 1A and Supplementary Fig. S3). After around 4.5 Ma, the ungulate faunas of both continents would fall back under the same functional configuration, a coherent structure that lasted until current days (Fig. 1A).

Humid, forested environments had been more common prior to the Middle Miocene Climatic Optimum33, but subsequently the world transitioned towards a drier and colder climate31. Increased aridity and declining CO2 atmospheric concentrations favored the expansion of C4-dominated grasslands and drove the development of morphological adaptations to these new environments39–41. In this new context, with the need to cope with more abrasive particles in their diet (from plant tissues, soil and grit on the plants), several herbivore lineages developed new dental functional structures and traits, including increased relative tooth height (hypsodonty), larger molar occlusal areas, flatter occlusal topography or increasing number of cutting, transversallophs, among others38,42–46. Consequently, new functionalities emerged as ecological opportunities unfolded, shaping the exploration of different regions within the functional space (Fig. 1B). The novel environmental settings provided ideal conditions for certain groups to diversify (e.g. Equids47), while the decline of forest-humid environments worldwide triggered a drastic decrease in global functional diversity38 (Fig. 2). The decline of the functional diversity that started 10 Ma has continued until the present.

The quaternary and the functional configuration of modern faunal assemblages
The persistent functional diversity loss that started around 10 Ma was significantly exacerbated around 2.5 Ma, as confirmed by our segmented regression analyses (Fig. 2). Such acceleration in functional disruption was likely due to the severe environmental changes in aridity, vegetation cover, and ecosystem productivity resulting from intensified fluctuations in glacial cycles31,34,48.

In contrast with other continents, a recovery trend in functional diversity is recorded in the Americas stemming from an increase in functional evenness (Fig. 2A, B). We propose that similar to the formation of the Gomphotherium land bridge, the creation of a natural corridor between North and South America was a key factor driving the rebound in functional diversity. Around 2.6 Ma, the formation of the Panama Isthmus triggered the Great American Biotic Interchange (GABI)49–54. The colonization of these new environments by large herbivores from both continents significantly increased their functional diversity, largely through enhanced evenness—species became more evenly distributed across various functional roles within the ecosystems of North and South America (Fig. 2B)51,55.

However, these changes in functional diversity do not appear to significantly alter the essential functional scaffold of large herbivore continental functional pools (Fig. 1A). From 4.5 Ma, UFFs remained globally unaltered, representing a functional structure stasis that persists to the present day. This stasis has endured even through periods of major abiotic crisis, such as the Messinian Salinity Crisis or the Plio–Pleistocene glaciations, which appear to have affected the functional diversity (Fig. 2). Unlike in past events, the last functional state is resilient enough to withstand the sustained loss of functional diversity while maintaining its essential functional configuration, highlighting the complex interactions between functional structure and diversity (Figs. 1A and 2).

Interestingly, the loss of megaherbivores (>1000 kg) did not entail shifts towards new UFFs (Fig. 1). The decline of these large herbivores began during the later stages of the Cenozoic38,56. This functional attrition is better exemplified by the Proboscideans, a once-successful group that reached a global diversity peak of 33 species around 3.2 Ma, just before the intensification of Plio-Pleistocene glaciations (likely a key driver of their decline; 38). The extinction of the largest ungulates accelerated after the Late Pleistocene (129 kya onward), with most species exceeding 1000 kg disappearing, leaving only a few survivors in modern terrestrial ecosystems24. As a result, ungulate communities and the ecosystems they inhabit have undergone significant transformations24,57–59. Nevertheless, only 18% of the FTs are lost during this last episode of megaherbivore decline (Supplementary Figs. 21, 22). Although megaherbivores have been a common feature of Cenozoic faunas, they have represented only a small fraction of functional types over the last 10 Myrs, comprising 15.45% overall (15.68% in Africa, module 7; 18.18% in Eurasia, module 8; 12.5% in the Americas, module 9). Their contribution in terms of species richness remains below 14% on average (17.49% in Africa, 18.28% in Eurasia, and 6.41% in the Americas). Thus, despite megaherbivore losses, the essential functional configuration of ungulate assemblages has remained stable in terms of UFFs over the last 4.5 Myr (modules 7–9, Fig. 1).

Summary
Our study portrays a dynamic functional system within large herbivore assemblages throughout the Cenozoic, characterized by periods of ecological stability and reassembly. Over broad temporal and spatial scales, large herbivores’ functional dynamics appear mainly driven by major abiotic changes. Such events pushed the system beyond points of no-return (tipping points), forcing the reassembly into new stable states. We found two major moments of functional reassembly of terrestrial ungulate faunas that occurred synchronously across all continents. During these events, the large herbivore functional structure followed analogous trajectories in the functional space at the continental level, displaying recurrent functional homogeneity following reassembly events.

Our results show that functional systems could be resilient enough to overcome the loss of ecological roles, maintaining their functional essence during moments of environmental change, revealing the intricate relationships between functional structure and diversity. Over the last 4.5 million years, large herbivore assemblages experienced a long period of functional structure stability, overcoming even several abiotic crises, which has lasted until the present. However, this period was marked by an unchecked decline in the functional diversity of these groups, resulting in an impoverishment of large herbivore ecosystems, notably exemplified by the loss of megaherbivores in nearly all terrestrial ecosystems across the Earth. These species performed functions in terrestrial ecosystems that are largely missing due to their extinction. Therefore, conservation efforts to preserve or reintroduce these species are essential for restoring the functionality of our ecosystems. If this functional impoverishment persists and worsens due to ongoing human impacts on the environment, it is uncertain when it will irreversibly affect the fundamental pillars of the large herbivore functional pool. This situation could potentially lead to a major third global tipping point and a new functional reassembly of ungulate faunas.

Methods
Database
We built a database for all the species in the largest orders of large herbivores living in earth current ecosystems (Artiodactyla, Perissodactyla, and Proboscidea) during the last 60 Ma. First, we compiled information for these species through a review of the primary literature along with data available in the NOW database60 and the Paleobiology Database (PBDB)60,61. Second, following Cantalapiedra et al.38, each species occurrence was checked in detail by experts in their taxonomic groups, looking for valid taxonomic assignments and excluding synonyms. Only records identified at the species level were retained. Third, we used the International Chronostratigraphic Chart62 to build a consistent scheme for temporal ranges across our dataset. We checked each record in order to improve their temporal resolution through a deep survey of the literature. In some cases, there were occurrences with a broad temporal range assigned (i.e., old publications of sites with scarce stratigraphic information, early sites with a wide range of temporal assignment, or incorrect entries in the databases). When we could not find better age ranges in the literature we excluded these records from our database. We decided to exclude the occurrences with a temporal range exceeding 5 Myr in the Neogene and 10 Myr in the Paleogene. For the Plio-Pleistocene records, we manually selected the ones with broad ranges and improved them using the most recent literature.

Since we used a combination of NOW and PBDB data, we obtained duplicated records for the species in our dataset. We chose between these records based on manual, case-by-case evaluation criteria of occurrence number and overall occurrence temporal precision. As a result, after all of our criteria and review procedure was applied, we obtained a final dataset that contains 22,028 occurrences of 3046 species distributed in 10,680 sites worldwide.

Our aim was to characterize the succession of functional assemblages over evolutionary timescales. To do so, we gathered information for 13 dental traits and body size for the species in our dataset. Based on ref. 63, we compiled dental trait information for: tooth shape, degree of hypsodonty, cusp shape, number of buccal cusps, number of lingual cusps, number of longitudinal lophs, number of transverse lophs, horizodonty, as well as the presence of acute lophs, obtuse lophs, structural fortification cusps, occlusal topography, and coronal cementum (see Zliobaite et al.46, Žliobaitė et al.63, Fortelius et al.64 for a detailed description of these traits). These dental traits are fundamental because they capture the relationship of the species with environmental conditions, especially being a good predictor of net primary production or precipitation63. Besides, we classified species in eight body size categories modified from65,66: B (<1kg), C (1–10kg), D (10–45kg), E (45–90kg), F (90–180kg), G (180–360kg), H (360–1000kg), I (1000–10,000kg) and J (>10,000kg). We excluded the original category A (0–0.1kg), because we did not have species in this body size range in our database. Overall, these traits reflect several facets of the functional role of species in ecosystems, such as habitat use, trophic level, range size, energetic requirements, and resource use10,67.

Spatio-temporal standardization
Working with large databases in paleontology has some difficulties, especially regarding the duplication or wrong assignment of locality names and the heterogeneity of stratigraphic ages. In order to prevent this, when we built our occurrence matrix we merged the localities using two criteria. First, we aggregated the localities by continent (Africa, America, Asia, and Europe). To characterize macroevolutionary dynamics of Cenozoic faunas, a continental scale allowed us to track changes worldwide and check the correlation of these changes with crucial abiotic shifts in Earth’ history such as plate tectonics or climate change. Second, we aggregated localities into temporal stages defined by the International Chronostratigraphic Chart stages (66, 61.6, 59.2, 56, 47.8, 41.2, 37.71, 33.9, 27.82, 23.03, 20.44, 15.97, 13.82, 11.63, 7.246, 5.333, 3.6, 2.58, 1.8, 0.774, 0.129 and 0.0117 Myr). The aggregation of localities by continent and stage yielded a total of 78 continent-stage assemblages, which are the analytical basis for the study. The use of continent-stage pools helped to avoid potential errors in the age assignments of the fossil localities.

Characterization of ungulate function faunas
Since we consider body size and tooth shape fundamental traits to define the species’ role in the ecosystems, we excluded the species with no information about these traits (this was the case for only 34 species, representing less than 0.01% of the entire species pool). The final dataset contained functional information for 3012 species, which formed the basis for all the analysis. To study the role of these species in their ecosystems, we had to reduce the dimensionality of our large mammal trait dataset. We computed Gower distances between species trait dataset using the daysy function in R68. During the computation of functional distances between species, each trait was given a different importance, so that dental traits would not have a disproportional weight in the total difference among species. To do so, we weighted two of the traits to 1 (body size and tooth shape), because they were the most representative traits of the species’ relationship with their environment, and assigned a weight of 1/6 to the remaining 12 traits (hypsodonty, cusp shape, buccal cusps, lingual cusps, longitudinal lophs, transverse lophs, horizodonty, acute lophs, obtuse lophs, structural fortification cusps, occlusal topography, and coronal cementum), thus summing up to 2, which is the result of the other traits combined.

We used the resulting distance matrix in a Principal Coordinates Analysis (PCoA). Based on the functional morphospace obtained from the PCoA, we applied a k-means to define groups of species with similar ecological roles: functional types (FT). Species closer in the functional morphospace will be grouped in the same FT. We decided on a maximum of 406 FTs, which was the number of unique combinations of the traits in our dataset (also called functional entities Mouillot et al.69). Our approach resulted in a final dataset with 406 FT distributed along 78 continent-stage pools, which were the basis for our analysis.

Network analysis
Following Blanco, et al.5, we used the developed Network Analysis approach to track the dynamics in the ecological assembly of large mammal continental functional pools worldwide over the last 65 Myr. We used this approach to calculate two key characteristics: affinity and fidelity. Affinity refers to the extent to which a functional type occurs in a high proportion of assemblage within a given module. In this context, a functional type with high affinity will be strongly associated with a particular module, here referred to as the ungulate functional fauna (UFF). Fidelity, on the other hand, measures the extent to which a functional type is predominantly restricted to a single module, indicating its specificity to that UFF (see methods).

Network input data. To analyze the temporal structuring of both functional and taxonomic communities, we implemented a community detection framework from network theory70,71. This technique enables the identification of clusters of fossil assemblages that exhibit similar ecological roles (functional types) or taxonomic makeup, independent of their chronological order, thus offering a window into how community structure evolved functionally over time. We translated the fossil occurrence data into bipartite networks composed of two distinct node sets: one representing continent-stage bins and the other either functional types or taxa. Connections between these nodes indicate the presence of a particular functional type (or taxon) in a given continental context. For functional networks, we further weighted these connections by counting the number of species within each functional type occurring in that bin, providing a measure of their relative prominence. In contrast, taxonomic networks were constructed using unweighted links due to the absence of data on species-level abundance.

Community detection algorithm. After constructing the bipartite networks, we applied the Infomap algorithm to identify communities or modules within them72,73. These modules correspond to sets of continent-stage bins that consistently co-occur with specific functional types or taxa, effectively grouping localities based on shared ecological or taxonomical composition. Infomap performs a joint classification of both node types—continent-stage and functional or taxonomic elements—into coherent clusters74. This method is grounded in the information-theoretic concept of minimum description length, which posits that the most structured representation of data is also the one that enables the highest degree of compression75. In our context, the regularities are expressed as modules composed of densely connected continent-stage and functional (or taxonomic) nodes, and the goal of Infomap is to identify the configuration that yields the shortest possible description of the network73,74. The algorithm proceeds heuristically, iteratively testing node reassignments and accepting changes that reduce the overall code length. Infomap was chosen due to its superior performance relative to other available community detection approaches74,76. We executed 10,000 runs of the algorithm and selected the partition with the lowest code length as the optimal solution73. This extensive run set ensured a comprehensive exploration of possible solutions, and the resulting modules were found to be robust across repetitions and consistent with those derived from alternative detection methods (see Sensitivity analysis to community detection).

To visualize the structure and timing of these communities, we plotted the modules identified in both the taxonomic and functional networks over geological time (Fig. 1 and Supplementary Fig. 1).

Comparing functional and taxonomic assembly. To examine how community assembly unfolded from both functional and taxonomic perspectives, we constructed a taxonomic counterpart to the functional dataset. This involved using the same continent-stage bins as those applied in the functional analysis, ensuring a direct comparison between both levels of organization. The resulting taxonomic dataset comprised 3046 species distributed across 78 continent-stage intervals (Supplementary Fig. 1).

Node characterization. We evaluated the importance of each node—whether functional type, taxon, or continent-stage—in shaping the structure of its assigned module using the IndVal index74, which integrates two key metrics: affinity and fidelity (see Supplementary Data 1). Affinity (Ai) reflects how concentrated a node’s connections are within its own module. It is calculated as the ratio of links a node has with others in the same module (Xi) to the total number of nodes in that module (Z):

Ai = Xi / Z (1)

For example, a functional type or taxon that appears frequently across continent-stages within its module will exhibit a high affinity.

Fidelity (Fi), by contrast, quantifies how exclusive those connections are to the module in question. It is defined as the number of links a node has within its module (Ni) divided by its total number of links across all modules (L):

Fi = Ni / L (2)

A node confined to continent-stages within its own module would thus have high fidelity. The final IndVal score is obtained by multiplying these two components: IndVal = Ai × Fi. High values of this index indicate nodes that are both widespread within and specific to their module, highlighting their significance in defining that community structure.

Module characterization
We calculated the IndVal index to study the importance of the nodes in the definition of their modules. Here, we used this index to extract the most singular FTs in the modules. Through a 0.95 quantile, we extracted these singular FTs for each module. Then, we look inside these FTs to study the most abundant trait states in their species composition functional trait pool (Supplementary Figs. 4–17).

Functional space
Beta diversity. In order to track the functional evolution of large herbivores, we calculated the difference in functional diversity between the different continent-stage. To do so, we used the betapart package in R77. This package calculates multi-sites distances using pairwise dissimilarities (in our case using as a basis the distance matrix between species in each bin). The advantage of using this package is that it allows to calculate the turnover and nestedness, components of functional beta diversity. Here, we used turnover as the difference of functional types composition between continent-stages, and nestedness when functional types assemblages in functional type-poor continent-stage are a subset of the assemblages in the more functional type-rich continent-stage. We ran the functional.beta.pair function, which computes functional dissimilarities based on volume of convex hull intersections in a multidimensional functional space. Dissimilarity matrix accounting for functional turnover, was measured as Simpson-derived pairwise functional dissimilarity, nestedness-resultant functional dissimilarity and functional beta diversity, was measured as the nestedness-fraction and beta diversity of Sorensen-derived pair-wise functional dissimilarity. We plotted the change in these three measures over time, calculating the change from the previous temporal bin for each continent (Supplementary Fig. 2).

We applied a non-metric multidimensional scaling (NMDS) to reduce the dissimilarity variability in two axes to plot the functional trajectories of continental assemblages over time in two dimensions. In this case, we decided to use the turnover distance matrix for the main part of the discussion and figures, as our interest is to track the functional change over time, which is the key information in order to study the dynamics of these faunas. We plotted the results using a temporal color-scale and module scale colors, obtained from the network analysis, for each continent (Fig. 1B and Supplementary Fig S19). This figure shows the different trajectories followed by the continents over time and how they explore similar regions of the functional space (Fig. 1B). The contemporaneous exploration of similar regions of the functional space (Fig. 1B) by different continents corresponds to periods of functional structure homogeneity (Fig. 1A).

Functional diversity
Once we characterized the functional structure and its evolution using network analysis (see above), we looked into the changes in functional diversity. This allows us to study the patterns that affect the functional diversity of these faunas and compare them with the dynamics of functional structure assembly. To do so, we calculated the functional diversity for all the temporal bins separated by each continent (continent-stage) in our database. We computed this for each continent-stage based only on the functional types (i.e. groups of species with similar ecological roles, see above) belonging to the same module, from network analysis, as the continental bin.

We first use Shannon’s diversity index, which considers both the richness of functional types and the even distribution of species into functional types:

H = - Σ R pi ln(pi) (1)

where R, is the number of functional types (also referred to here as functional richness), and pi is the proportion of species belonging to the i functional type. High values of the Shannon index indicate high continental functional richness (number of functional types) and a more homogeneous distribution of species along the functional types found in a particular continent-stage. Thus, high Shannon values indicate high values of ecological disparity (high number of functional types) coupled with high values of ecological redundancy (species evenly distributed along the functional types).

Then, we also considered each constituent part of the Shannon index (i.e. richness and evenness) independently. To measure evenness, we used the Pielou’s index, calculated as,

J = H / ln(R) (2)

where H’ comes from Shannon’s diversity index. High values of evenness indicate high ecological redundancy, which means that species tend to be evenly spread across functional types (i.e. ecological roles).

For the calculation of these indexes we used only continent-stage assemblages with more than 20 species, following the sensitivity analysis to sampling performed. To account for the sampling effect of the differential number of occurrences in the different continent-stage, we extracted the residuals from linear models where indices were a function of sampling. We used the residuals from the best model (occurrences*continent) for the subsequent analysis (Supplementary Table 2).

Linear models of functional diversity
In order to correct for the sampling effect of the differential number of occurrences in the different continent-stage78,79, we extracted the residuals from linear models where indices were a function of sampling80. First, we fitted linear models where functional diversity is a function of the number of occurrences while allowing the relationship to change among continents (index ~ occurrences*continent). Second, we fitted linear models where functional diversity is a function of the number of occurrences and the quadratic of this number (index ~ occurrences*occurrences2). The best model for each index was picked based on AICc scores (Akaike’s information criterion corrected for small sample size): a more complex model was only preferred if the AIC score was smaller by at least two units. For the three indexes (Shannon, Evenness and Richness), the occurrences*continent model resulted in the best model, and we used the residuals from it for the subsequent analysis (Supplementary Table 2).

Moments of abrupt physical and environmental change have been linked to reorganizations of functional structure5. To explore this, we fitted linear models using the corrected functional diversity indexes (see above) as a function of the continent-stage age as an explanatory variable (corrected index ~ age). We used the function segmented in R to look for the breakpoints in the residuals of the linear model. We plotted along with the corrected index evolution through time (Fig. 2). Besides, to see how the functional diversity trends affected at geographical level and functional structure, we plotted the residuals indexes by continents through time (Fig. 2B).

Functional distance among continents over time
In order to identify functional resemblance among continents over time, we calculated the functional distances as the distance of the continent-stage in the functional dissimilarity space. We plotted the resulting measure over time using continent pairs (Supplementary Fig. 3).

Sensitivity analysis
Sensitivity analysis on the k-means groups. We tested the selection criterion of 406 FT based on the number of unique combinations of functional traits. We ran the k-means analysis selecting 100, 200, 300, 400 and 500 groups, and we repeated 100 times. We used this to build the analysis database and run infomap72. As a result, we plotted the quality code length (measure of modularity) and the number of modules obtained from the network analysis for all k-means randomization (Supplementary Fig. 18). In both cases the variable values become stable around 400 FT.

Sensitivity analysis to sampling. When working with paleontological data, we should always be aware of potential sampling biases that could drastically affect the results of our analyses (i.e. some epochs being better sampled, countries or regions with difficult access to conduct fieldwork, historical preferences…). In other words, the disparate number of species across the analyzed bins could impact the estimated functional similarity between them. To test this, we conducted a sensitivity analysis where we performed a rarefaction procedure79,80 on the species present in each continent-stage. We randomly selected 20 species from each bin, and then we ran the network and beta diversity analysis (see above for detailed explanation) on the rarefied dataset. We repeated the procedure 100 times. An averaged distance matrix from the 100 draws was used to estimate functional turnover between bins and to compare with our main results (Supplementary Fig. 19). The results from the rarefactions were highly congruent with those derived from the full dataset.

Sensitivity analysis to the effect of aggregation. To determine whether the modular patterns found in our functional networks were genuinely reflective of ecological structure—or merely an artifact of grouping species into functional types—we performed a null model analysis. Specifically, we generated randomized networks in which species were reassigned to functional types at random, while preserving the original number of species per functional type. We then compared the modularity of these null networks to that of the observed network. Modularity was assessed within the Infomap framework by calculating the relative code length: the ratio of the code length for the partitioned network (CLm) to the code length of the same network without any modular structure (CL). A relative value close to 1 indicates low modularity, whereas values near 0 suggest strong community structure. To express modularity in a more intuitive way, we used its complement:

M = 1 - CLm / CL (3)

Here, CLm refers to the code length of the modular network, and CL to that of the unpartitioned network. To assess statistical significance, we computed a p-value based on the proportion of 100 randomized networks (plus the observed one) that showed modularity equal to or greater than that of the empirical data (Supplementary Fig. 20).

Sensitivity analysis to community detection. We evaluated two main sources of uncertainty in our detection of network communities. First, since algorithms like Infomap rely on stochastic optimization to identify network partitions, multiple runs can yield slightly different solutions. This raises two important questions: how many runs are needed to reliably capture the full landscape of possible solutions, and how stable the resulting modules are across near-optimal partitions81. To assess the first, we conducted a 10-fold cross-validation, estimating the likelihood that a randomly held-out partition would share at least 75% similarity with one from the training set. Results showed that executing Infomap 10,000 times yielded an almost complete exploration of the solution space, with test partitions almost always exhibiting similarity above 0.75—both for functional and taxonomic networks.

To evaluate the second aspect—module stability—we examined whether the communities found in the best partition reappeared in the rest of the 9999 solutions. For each module, we computed the probability of retrieving a similarly composed group in other runs, using two thresholds of similarity: a standard one (0.5), implying one-to-one correspondence, and a stricter threshold of 0.7581. In both cases, the core modules remained highly reproducible across solutions, with probabilities close to 1, indicating strong consistency (Fig. 1 and Supplementary Table 1).

Sensitivity analysis to accuracy of network partition vs functional distances. We ran a Permutational multivariate analysis of variance (PERMANOVA) to test whether the module partition from network analysis reflects the functional distances between the continent-stage. We obtained that the module partition is similar to functional distances by a statistically significant p > 0.001.

Potential limitation of the spatio-temporal standardization
Our approach focuses on aggregating localities into a spatio-temporal sequence that allows us to study the evolution of the functional structure of large mammals at the continental level. We acknowledge the limitations of this approach, particularly regarding spatial standardization. A finer-scale geographic partitioning (e.g., by country, region, or grid) could potentially yield different functional structure patterns. However, these results would reflect a different spatial scale and should be interpreted accordingly. Our findings capture broad continental-level patterns, which remain consistent across all sensitivity analyses performed (see above).

One limitation of our study is the merging of North and South America in the functional analyses. Due to the limited availability of data (both in occurrences and functional traits) for South American species (only 50 taxa), analyzing the continents separately would have resulted in highly uneven datasets, potentially biasing comparisons. South America’s fossil record is less extensively studied than North America’s, mainly due to preservation biases in tropical regions, further constrained by geopolitical and economic factors. Additionally, certain endemic South American herbivore groups, such as notoungulates, have been less studied due to the scarcity and fragmentary nature of their fossil record.

While this approach allows for more robust statistical comparisons, we acknowledge that North and South America had largely independent faunal histories prior to the Great American Biotic Interchange (GABI). The merging of these datasets may obscure some biogeographic patterns, particularly those associated with faunal exchange events. Future work incorporating a more complete South American fossil record, including a broader representation of endemic taxa, will help refine our understanding of the dynamics of functional diversity across both continents.

We used R software68 for all the analyses. All the scripts and database used are at: https://github.com/f-blanco/ecoherb and refs. 82,83.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a very interesting paper, where the authors use an approach based on network analyses for studying the faunal dynamics of ungulate continental assemblages (orders Artiodactyla, Perissodactyla and Proboscidea) during the Cenozoic and their reorganization in response to major tectonic and climatic events.\nThey find (i) a prevalence of long periods of ecological stasis; (ii) a decoupling of taxonomic and functional assembly dynamics; and (iii) two major shifts (tipping points) in the functional diversity of the guilds of primary consumers at 26-21 Ma (i.e., at the Oligocene-Eocene transition) and after 10 Ma (i.e., at late Micene times), respectively.\nAccording to their interpretation, the first tipping point relates to the formation of a land bridge between Eurasia and Africa, which led to a transition towards a prevalence of large browsers with mid- to high-crowned molars. The second resulted from aridification and the spread of C4-dominated vegetation, which triggered the emergence of an ungulate fauna characterized by grazers and browsers with high and low crowned teeth, respectively.\nThe authors cite a couple of papers by Figueirido et al. in PNAS (their refs. 10 and 11) as examples of studies addressing faunal shifts in the fossil record that “have been traditionally tackled from a taxonomic perspective”. However, this applies to Figueirido et al. (2012), who analyzed the evolutionary faunas of the North American Cenozoic using a taxonomic perspective, but not to Figueirido et al. (2019), a paper where the approach was based on a taxon-free ecomorphological characterization of these evolutionary faunas. Moreover, the latter study included analyses of functional diversity of large mammals (not only ungulates) based on tooth types, as those performed in the current study by Blanco and colleagues, but also on the frequencies of locomotor types and body size classes among mammalian assemblages, an approach that provided a wider and more balanced view of the changes through time in ecomorph diversity. While acknowledging that ungulates are primary “ecosystem engineers”, it would be worth considering here if a better picture of past ecosystems could emerge from an analysis of the whole community of large mammals.\nFinally, a couple of relatively minor points that should also be addressed:\n1.- After the Mid-Miocene Climatic Optimum, the transition towards a drier and colder climate drove the development of hypsodont teeth in ungulates (and also of hypercursorial adaptations, as noted by Figueirido et al., 2019), which allowed them to process the more abrasive foodstuffs available in open habitats. However, this dietary adaptation is related to (i) the presence of external abrasives (grit) on the plants from open habitats and (ii) the consumption of soil rather than to the presence of internal abrasives (silicophytolithes) in the herbaceous plants, as demonstrated by Mendoza & Palmqvist (2008; J. Zool., 274: 134-142) and Damuth & Janis (2011; Biol. Rev., 86: 733-758).\n2.- The authors comment on only one major dispersal event between Eurasia and Africa at the Oligocene-Eocene transition. However, in their analysis of North American evolutionary faunas of large mammals Figueirido et al. (2012) also indicated that the rise of the Oligocene fauna may relate to the major immigration event from Asia in North America, which took place during the middle-late Eocene (i.e., at 37.5 Ma) following a sea level fall at this time. Also, they indicated that the rise of the North American Pliocene fauna may be related with a major immigration that occurred near the end of the early Miocene, which suggests that this evolutionary fauna began to diversify at approximately 16.5 Ma, shortly after an extensive interchange between Eurasian and North American faunas, and that this immigration event may also have played a role in the decline of the Miocene fauna.", "reviewer_2": "I found the manuscript was well-written, concise and enjoyable to read. I believe it would be a great article for Nature Communications to publish. I felt it conveyed important findings about macroevolutionary changes in ungulates. Results and discussion were kept broad and logical without over-inferring driving impacts. The story was upheld of investigating continental scale changes through time without getting lost in the details. The figures are limited to those providing key information and are visually appealing. The methods in functional diversity and network analysis are sound and possible biases in the data were addressed appropriately in the supplementary material.\nThere are a few points I would like addressed before acceptance. I added to my review document including the combination of North and South America, missing traits, the clear definition of reconfiguration vs. reassembly and the interpretation of high functional diversity. I feel confident these concerns can be appropriately addressed with logical explanations and/or additional analyses.\nThank you for your contribution to science.", "reviewer_3": "Dear authors,\nI think that this article is very important and meet some vital aspects required for a good article. The methods are clearly stated and adapted to the question. Their limitations are considered and the author’s team overcame them with classic corrections methods. The dataset build is not new, it’s made of existing data from various sources (databases, literature) but they use it at a very large scale to interrogate crucial questions. How functional diversity has changed during the last 60 million years for ungulate all around the globe? By using various classical methods (e.g. clustering, models, ordination) on network of functional type they were able to detect significant breaks and tendencies in the continuum at global and continental scales. They discovered timing of emergence of modern patterns and tried to link these breaks, emergence and tendencies with known abiotic and biotic events. I really enjoyed this paper and congratulate the authors for this quality research. After answering to reviewer remarks and questions, I would recommend a minor revision. My remarks are dedicated to make the paper easier to understand, I would try to highlight the parts where I struggled and suggest potential arguments to make a better case of the observed patterns.\nMajor remarks\nFigures are very well done and most of them are easy to interpret and understand. However, I struggled to really understand Fig. S1. The caption could be extended to better explain its construction and/or interpretation. I had difficulties to understand the statement made about it line 66 (“Our results show a decoupling of taxonomic and functional”) while looking at the figure. Please, either extend the caption, or add some text to interpret this figure in the page 1 of the supplementary data (because I’m sure you have to keep the main text as concise as possible). I think that the statement made line 68-69 is very important (even if already known in others context/studies) and could be more easily understood if the description of Fig. S1 would be extended in supplementary files.\nL. 291-293: “The aggregation of localities by continent and stage yielded a total of 78 continent-stage assemblages, which are the analytical basis for the study”. I think that’s totally needed and fair for a study at the temporal and spatial scale, but could you develop a bit the potential limitation of such approach in your method section (or in supplementary)? What could be the impact of such data agglomeration? For example, could different UFFs co-exist on one continent (e.g. SE Asia vs Central Asia) and thus create fake stability by experiencing different trajectories of changes through times?\nMinor remarks\nI think that the authors could have spend more time to compare known patterns of taxonomic diversity in mammals that took interest in clustering temporal fauna (e.g. chronofauna in both North America and Europe/Asia/Africa, see for example Figueirido, Gibert, Casanovas-Vilar) and look for breaks and tendencies but I understand that Nature format is very constrained. However, I think that the excellent work they did with functional diversity could have been even better if it would have been better compared to other dimensions of mammalian fossil biodiversity.\nL. 74-75 : “Since 10 Ma, large herbivore assemblages have experienced remarkable stability”.\nOne remark: it could be stated here that Asian UFFs are acting slightly differently than the others continental fauna.\nOne question: when looking at Fig. 1B, we can see how the average/median point of UFFs in NMDS is moving in the environmental space through time, but I struggle to understand how Fig. 1A and Fig. 1B relate.\nFor example, when looking at the orange/yellow section for Africa (i.e. the last 10 myr), we see strong displacement in the environmental space for the african UFF, but it corresponds to statis (or a single mode, turquoise-green) in Fig. 1A. In the same way, the early Cenozoic asian UFF relate to the first mode (i.e. red) and last for ~35 myr, but we can see in Fig. 1B this UFF moving a lot through the environmental fauna (i.e. black to red section of the line). Can you extend the explanation in supplementary data to understand the relationship between Fig. 1A and 1B. Fig. S19 could be used to explain this relationship with more details.\nL. 97-98 : “This was primarly driven by the emergence of functional dominance of FTs in larger body size categories in the Americas, in contract to the mid-sized herbivores that characterize the UFF of Eurasia (Fig. S17)”.\nCould it be a matter of the relationship between landmass and (paleo)latitudes? The larger southern extent of Asia could explain the survival of mid-sized herbivores and the continuous cooling of Cenozoic could act in very different direction depending on the latitudinal distribution of landmass, isn’t it? The argument developed in the rest of this section (from L.95 to L.119) could be completed by a discussion around landmass/biome (i.e. closed vs open forests) distribution. The two references associated with EOC transition (29, 40) are associated with very descriptive arguments (L.109-119) when explanatory argument (i.e. mechanism) could/should be invoked here (e.g. Latitude and landmass, Europe being restricted to high latitude and forming an archipelago making climate tracking more difficult for European fauna). On a more general statement, I think that this paper is sometimes a bit too descriptive, more potential mechanisms could be invoked to explain the observed patterns.\nL. 124-125: “The functional coalescence is reflected in the drastic reduction of the functional distance between Africa and Eurasia after 21 Ma”.\nVery interesting, could you just add the age uncertainties (i.e. the corresponding geological stage length if I correctly understood your methodology) in order to make this statement more subtle. The exact date of the main part of the interchange/faunal mixing between Africa and Eurasia remain debated and your sentence could imply that your study can precisely date this event.\nL. 298: “We excluded the species with no information about these traits”.\nDid you analyse/reanalyse/measure species traits or did you find the 13 traits being already described for the 3012 species included in your study?\nL. 306: “and assigned a weight of 1/6 to the remaining 12 traits”.\nHave you tried different weight, in other word, did you perform sensitivity analysis on this parameter ?\nL. 393-397: Thank for the clear interpretation of the Shannon index for FD. Would you say that high Shannon values in your study relate to strong ecosystem stability (high diversity + high redundancy)? The comparison of Shannon FD and Fig. 1 is not trivial.\nThank you again for your future answers and congratulation again for this very interesting study."}, "author_response": "Response to reviewers comments\nReviewer comments are highlighted in bold.\nResponses to reviewers are provided in regular text, following each\nrespective comment.\nReviewer #1:\nThis is a very interesting paper, where the authors use an approach\nbased on network analyses for studying the faunal dynamics of\nungulate continental assemblages (orders Artiodactyla,\nPerissodactyla and Proboscidea) during the Cenozoic and their\nreorganization in response to major tectonic and climatic events.\nThey find (i) a prevalence of long periods of ecological stasis; (ii) a\ndecoupling of taxonomic and functional assembly dynamics; and (iii)\ntwo major shifts (tipping points) in the functional diversity of the\nguilds of primary consumers at 26-21 Ma (i.e., at the Oligocene-\nEocene transition) and after 10 Ma (i.e., at late Micene times),\nrespectively. According to their interpretation, the first tipping point\nrelates to the formation of a land bridge between Eurasia and Africa,\nwhich led to a transition towards a prevalence of large browsers with\nmid- to high-crowned molars. The second resulted from aridification\nand the spread of C4-dominated vegetation, which triggered the\nemergence of an ungulate fauna characterized by grazers and\nbrowsers with high and low crowned teeth, respectively.\nWe sincerely appreciate the positive feedback from Reviewer #1 regarding\nour methodological approach and key findings.\nThe authors cite a couple of papers by Figueirido et al. in PNAS\n(their refs. 10 and 11) as examples of studies addressing faunal\nshifts in the fossil record that “have been traditionally tackled from\na taxonomic perspective”. However, this applies to Figueirido et al.\n(2012), who analyzed the evolutionary faunas of the North American\nCenozoic using a taxonomic perspective, but not to Figueirido et al.\n(2019), a paper where the approach was based on a taxon-free\necomorphological characterization of these evolutionary faunas.\net al. (2019) employed a taxon-free ecomorphological characterization of\nevolutionary faunas. The inclusion of this reference in that context was\nan oversight, and we have corrected it by removing Figueirido et al.\n(2019) from the citation.\nMoreover, the latter study included analyses of functional diversity\nof large mammals (not only ungulates) based on tooth types, as\nthose performed in the current study by Blanco and colleagues, but\nalso on the frequencies of locomotor types and body size classes\namong mammalian assemblages, an approach that provided a\nwider and more balanced view of the changes through time in\necomorph diversity. While acknowledging that ungulates are\nprimary “ecosystem engineers”, it would be worth considering here\nif a better picture of past ecosystems could emerge from an\nanalysis of the whole community of large mammals.\nWe appreciate the reviewer’s insightful comment regarding the potential\nbenefits of analyzing the entire large mammal community. However, we\nfocus on ungulate faunas for several reasons:\n1. Ecological significance: Ungulates have played a dominant role in\nshaping terrestrial ecosystems since the K-Pg mass extinction,\nacting as primary ecosystem engineers (Brits et al 2002, van der\nWaal et al. 2011, Smith et al.m 2016).\n2. Fossil record availability: Their widespread distribution and\nabundant fossil remains throughout the Cenozoic enable a large-\nscale macroevolutionary analysis (Zliobaite and Fortelius 2022). In\ncontrast, other mammal groups (e.g., Carnivora) tend to have a\nmore limited and sparse fossil record, making a comparable\nanalysis challenging.\n3. Temporal continuity: We focus on the orders Artiodactyla,\nPerissodactyla, and Proboscidea because they are the only large\nherbivore clades with both extant representatives and a fossil\nrecord extending back 60 million years. This continuity allows us to\nanalyze long-term trends in functional structure.\n4. Dataset scope: Thanks to these factors, we assembled a robust\ndataset comprising 22,028 occurrences of 3,046 species from\n2\n10,680 fossil sites worldwide, providing an unprecedented\ntemporal and spatial resolution for studying functional diversity in\nlarge herbivores.\nWhile incorporating additional mammalian groups could offer a broader\nperspective on ecosystem dynamics, the limitations in their fossil\nrepresentation would compromise the feasibility and comparability of\nsuch an analysis. Future progress in taxonomic and systematic\nresearch, along with new fossil discoveries from ongoing and future\nexcavations, will be crucial for expanding our understanding of the entire\nlarge herbivore guild and even mammals as a whole. As more data\nbecome available, such analyses may become feasible, potentially\noffering deeper insights into past ecosystems.\n1.- After the Mid-Miocene Climatic Optimum, the transition towards\na drier and colder climate drove the development of hypsodont\nteeth in ungulates (and also of hypercursorial adaptations, as noted\nby Figueirido et al., 2019), which allowed them to process the more\nabrasive foodstuffs available in open habitats. However, this\ndietary adaptation is related to (i) the presence of external\nabrasives (grit) on the plants from open habitats and (ii) the\nconsumption of soil rather than to the presence of internal\nabrasives (silicophytolithes) in the herbaceous plants, as\ndemonstrated by Mendoza & Palmqvist (2008; J. Zool., 274: 134-\n142) and Damuth & Janis (2011; Biol. Rev., 86: 733-758).\nWe appreciate the reviewer’s insightful comment regarding the factors\ninfluencing dental abrasion in ungulates. To address this, we have\nincorporated the relevant information on line 190, emphasizing the role\nof external abrasives (grit) and incidental soil consumption in driving the\nevolution of hypsodonty:\n“In this new context, with the need to cope with more abrasive particles\nin their diet (from plant tissues, soil and grit on the plants), several\nherbivore lineages developed new dental functional structures and traits,\nincluding increased relative tooth height (hypsodonty), larger molar\nocclusal areas, flatter occlusal topography or increasing number of\ncutting, transversal lophs, among others (45-50). ”\nheight, larger molar occlusal surfaces, flatter occlusal topography, and a\ngreater number of transverse lophs, among other adaptations in\nungulate faunas.\n2.- The authors comment on only one major dispersal event\nbetween Eurasia and Africa at the Oligocene-Eocene transition.\nHowever, in their analysis of North American evolutionary faunas of\nlarge mammals Figueirido et al. (2012) also indicated that the rise of\nthe Oligocene fauna may relate to the major immigration event from\nAsia in North America, which took place during the middle-late\nEocene (i.e., at 37.5 Ma) following a sea level fall at this time. Also,\nthey indicated that the rise of the North American Pliocene fauna\nmay be related with a major immigration that occurred near the end\nof the early Miocene, which suggests that this evolutionary fauna\nbegan to diversify at approximately 16.5 Ma, shortly after an\nextensive interchange between Eurasian and North American\nfaunas, and that this immigration event may also have played a role\nin the decline of the Miocene fauna.\nWe thank the reviewer for these comments. We have incorporated the\nsuggested information into line 94: “Functional innovations were likely\nintertwined with biogeographic events, contributing to this enrichment.\nFor example, around 37.5 Ma, a significant immigration of large\nmammals from Asia to North America occurred following a sea-level\ndrop (10).”\nAnd line 165: “Interestingly, this decline follows an immigration event\nbetween Eurasia and North America at the end of the early Miocene,\nsome 16.5 Ma (10), suggesting that faunal dispersals do not necessarily\nenhance functional diversity but may instead reduce the functional\nrichness of endemic American functional types.”\nI found the manuscript was well-written, concise and enjoyable to\nread. I believe it would be a great article for Nature Communications\nto publish. I felt it conveyed important findings about\nmacroevolutionary changes in ungulates. Results and discussion\nwere kept broad and logical without over-inferring driving impacts.\nThe story was upheld of investigating continental scale changes\nthrough time without getting lost in the details. The figures are limited\nto those providing key information and are visually appealing. The\nmethods in functional diversity and network analysis are sound and\npossible biases in the data were addressed appropriately in the\nsupplementary material. There are a few points I would like addressed\nbefore acceptance. I added to my review document including the\ncombination of North and South America, missing traits, the clear\ndefinition of reconfiguration vs. reassembly and the interpretation of\nhigh functional diversity. I feel confident these concerns can be\nappropriately addressed with logical explanations and/or additional\nanalyses.\nThank you for your contribution to science.\nWe sincerely appreciate the reviewer’s thoughtful and encouraging\ncomments on our manuscript. We are pleased that they found the\nmanuscript well-written, concise, and engaging, and that our approach to\nstudying macroevolutionary changes in ungulates was clear and\nappropriately framed. We also appreciate the reviewer’s positive\nassessment of our methodological framework and data considerations.\nWe provide detailed responses below, along with clarifications and\nchanges in the manuscript.\nIntroduction\nLines 48 – 53: I think it would be helpful for the geographic\ndistribution to be stated somewhere in the description of the data\nhere. As a reader, it was the first question that popped in my head.\nJust adding “global” would be sufficient.\nResults and Discussion\nLines 78: Can you clarify why North and South America were\nanalyzed together? Due to the independence of the two faunas, it\nseems more reasonable that they would be analyzed separately. It\nseems important to note if their separation alters the overall\nfunctional structure. Moreover, you would potentially expect a shift\nat GABI.\nWe merged North and South America in our analysis due to limited data\navailability for South America, particularly for the three focal orders\n(Artiodactyla, Perissodactyla, and Proboscidea), which include only 50\nspecies there. Several factors explain this limitation: 1) The fossil record\nof South American mammals is not as extensively studied as that of\nNorth America, which benefits from a long tradition of paleontological\nresearch. Biases in the fossil record of South America, particularly in\ntropical regions, as well as geopolitical and economic factors, contribute\nto this disparity. 2) South America hosted other large herbivore groups,\nsuch as notoungulates, which were not included in our analysis. We\nacknowledge that future excavations and ongoing studies on endemic\nSouth American herbivores may allow for their inclusion in future\nanalyses, providing a more comprehensive understanding of the\necosystems in which these groups lived. We have added text\nacknowledging this limitation in the new section “Potential limitations of\nthe spatio-temporal standardization”.\nLines 68 – 73: What defines a functional reconfiguration? And how\ndoes that vary from tipping points? You say there is a functional\nreconfiguration every ~6 Myr, but there are long periods of\necological stasis and there are two periods of reassembly. Could\nyou clarify how those different types of events are defined? I\nbecame confused when you said “with a first reconfiguration in\nUFFs around 40 Ma (line 86)\nWe define functional reconfiguration as any module change identified\nthrough network analysis. Tipping points, on the other hand, refer to\nmoments when these changes occur simultaneously across all\n6\ncontinents (i.e., a simultaneous module change worldwide). To clarify\nthis distinction, we revised the sentence in line 71:\n“Large herbivore assemblages experienced long periods of ecological\nstasis punctuated by two episodes of global and almost synchronous\nfunctional reassembly (tipping points, Fig. 1A). These tipping points\noccurred around 21 and 10 Ma, likely as a response to major abiotic\nevents. Since 10 Ma, large herbivore assemblages have experienced\nremarkable stability.”\nIn this context, we mean that the functional structure underwent a\nmodule change every ~6 million years on average, while the taxonomic\nstructure changed more rapidly, with a module change every 1.7 million\nyears on average. We revised the sentence for clarity as follows in line\n67:\n“On average, taxonomic reassembly (module change, see Fig. S1)\noccurred every 1.7±0.31 Myr, whereas functional reassembly took place\nevery 6.25±1.6 Myr. “\nRegarding UFFs (Ungulate Functional Faunas), we have clarified their\ndefinition since the introduction: ungulate continental assemblages with\nshared configurations (presence, relevance, and fidelity) of functional\ntypes (see line 67):\n“We then applied Network Analysis to identify ungulate continental\nungulate functional faunas (UFFs)(5), defined as unique associations of\nfunctional types. Notably, UFFs are not defined solely by the presence of\nFTs, but by their relevance, affinity, and fidelity in the faunas (see\nsupplementary methods).”\nUFFs are the modules in our network analysis (see Fig. 1A). The first\nfunctional reconfiguration in the UFFs (module change) occurred around\nmodule.\nWe hope that these clarifications make the terminology more precise\nand easier to follow.\nit twice because I thought you were suggesting that functional\ndiversity was recovering in those continents.\nOk, we have replaced \"recovered\" with \"identified\" to avoid any\nconfusion.\nLines 197 – 198: If you were looking at functional diversity of North\nand South America together before GABI, you wouldn’t suspect a\nchange in the overall functional pool following the interchange.\nAs the reviewer pointed out, the functional pool (i.e., functional richness,\nor the number of functional types) was not the main driver of the\nfunctional diversity increase (Fig. 2B). In fact, our results indicate a\ndecrease in functional diversity during that period. The observed\nincrease in functional diversity is primarily due to a rise in functional\nevenness (Fig. 2B). This suggests that after the GABI, there was an\nincrease in redundancy, with species becoming more evenly distributed\nacross functional types (i.e., ecological roles). While functional richness\ndecreased, the new species that evolved during this period occupied\nsimilar ecological roles to those that existed previously.\nLines 198 – 199: The functional structure has remained the same\nfollowing all the extinctions and extirpations? If you mean in the\nmost recent interval (encompassing the last 10,000 years), I would\nsay that, not until present day. The way I understand it, is all large\nungulates that have existed on the continent during any part of that\ninterval is included in determine the functional structure.\nOur results show a continuous decline in functional diversity from 4.5\nMa, primarily driven by a decrease in functional richness. However,\nfunctional evenness slightly increased during this period. While\nfunctional types are being lost, new species that evolved during this time\noccupied pre-existing functional roles, leading to increased redundancy\nand potentially contributing to the stability of the functional structure\n(modules).\nnecessary for functional types to go extinct, but these extinctions must\ninvolve types that are highly connected within a module, with minimal\nconnections to other modules (high affinity and fidelity, showing a high\nIndVal index value, see method section in the supplementary). These\ntypes are critical for defining modules due to their exclusive connections\nwithin specific groups. Furthermore, a recent study using a different\nnetwork analysis approach (bayesian network analysis of food webs)\nfound persistent stability throughout the Pleistocene in Eurasian\nmammals (Bekeraite et al. 2024), supporting our findings at a global\nscale.\nMethods\nLines 297 – 298: It is unclear to me if “no information about these\ntraits” means the species needed to be missing all of the traits or if\nthey are missing a information about a single trait. If species were\nkept that were missing 1 or multiple traits, could this bias your data\nand was a sensitivity analysis used to address any possible\nbiases?\nWe discarded species lacking data on body size or tooth shape, as\nthese traits are fundamental to defining the ecological roles that species\nplayed in their ecosystems. Only 35 species out of more than 3,000 were\nmissing these traits, representing less than 0.01% of the entire species\npool. We believe this minimal exclusion does not introduce significant\nbias to our analysis. Although no sensitivity analysis was conducted\nspecifically for this issue, given the small proportion of species affected,\nwe do not expect this omission to substantially impact the results.\nFor clarity, we have slightly rephrased the sentence as follows:\n“Since we consider body size and tooth shape fundamental traits to\ndefine the species’ role in the ecosystems, we excluded the species with\nno information about these traits (this was the case for only 34 species,\nrepresenting less than 0.01% of the entire species pool). “\nLines 395 – 397 and 404 - 406: a homogenous distribution in\nfunctional space does not necessarily represent high redundancy.\ndistribution homogeneity across functional space. So I’m not sure I\nbelieve high Shannon values always represent ecological\nredundancy, unless I am misunderstanding something that\naddresses this in the equation. Is there any reason to think that\ngeographic distribution of available occurrences could influence\nyour results during any time periods?\nTo calculate the Shannon index, we did not use the distribution in\nfunctional space. Instead, we applied the classical Shannon index\nequation, which incorporates both richness and evenness. High values\nof the Shannon index require both a high number of functional types and\na high degree of species evenly distributed among these types.\nTo account for potential biases in the geographic distribution of\noccurrences, we calculated the Shannon index for each continent-stage\nusing only the functional types (i.e., groups of species with similar\necological roles) within the same module, as identified by network\nanalysis, corresponding to the continental bin. Additionally, we limited\nour analysis to continent-stage assemblages with more than 20 species,\nas guided by the sensitivity analysis on sampling.\nTo further address sampling effects related to the differential number of\noccurrences across continent-stages, we extracted the residuals from\nlinear models where the indices were modeled as a function of sampling\n(see Supplementary Methods). The residuals from the best-fitting model\n(occurrences*continent) were then used in the subsequent analysis\n(Table S3).\nAll procedures are described in detail in the Methods section of the main\ntext and the Supplementary Methods.\nFigures\nIn Fig. 1A African data begins around 30 Ma, in Fig. S2 it begins\naround 25 Ma and in Fig. 2, African data does not start until 20 Ma.\nCan you explain why?\nBefore 30 Ma, African assemblages only included proboscideans.\nTherefore, we decided to exclude the earlier part of the African Cenozoic\n10\n(60-30 Ma) from the analysis to maintain a comparable framework with\nthe other continents.\nRegarding Fig. S2, the beta diversity difference is calculated by the\ndifference between one bin and the previous one. Since the first bin\nstarts at 30 Ma, we discarded this bin and began the calculation with the\nfollowing one (Chattian stage with a mean age of 25 Ma).\nFor Fig. 2, as explained above and detailed in the Methods and\nSupplementary Methods sections, we applied a filter requiring at least 20\nspecies per continent-stage for inclusion in the analysis. The first African\ntime bins did not meet this condition, so they were not included. The first\nbin that satisfied this condition is at 20 Ma (Burdigalian stage).\nDear authors,\nI think that this article is very important and meet some vital aspects\nrequired for a good article. The methods are clearly stated and\nadapted to the question. Their limitations are considered and the\nauthor’s team overcame them with classic corrections methods. The\ndataset build is not new, it’s made of existing data from various\nsources (databases, literature) but they use it at a very large scale to\ninterrogate crucial questions. How functional diversity has changed\nduring the last 60 million years for ungulate all around the globe? By\nusing various classical methods (e.g. clustering, models, ordination)\non network of functional type they were able to detect significant\nbreaks and tendencies in the continuum at global and continental\nscales. They discovered timing of emergence of modern patterns and\ntried to link these breaks, emergence and tendencies with known\nabiotic and biotic events. I really enjoyed this paper and congratulate\nthe authors for this quality research. After answering to reviewer\nremarks and questions, I would recommend a minor revision. My\nremarks are dedicated to make the paper easier to understand, I\nwould try to highlight the parts where I struggled and suggest\npotential arguments to make a better case of the observed patterns.\nWe sincerely appreciate the reviewer’s detailed and positive feedback on\nour manuscript. We are pleased that they found our study important,\nmethodologically sound, and relevant to understanding long-term\nchanges in ungulate functional diversity.\nWe have carefully addressed each of their points and incorporated their\nsuggestions where applicable to strengthen the presentation of our\nfindings.\nFigures are very well done and most of them are easy to interpret\nand understand. However, I struggled to really understand Fig. S1.\nThe caption could be extended to better explain its construction\nand/or interpretation. I had difficulties to understand the statement\nmade about it line 66 (“Our results show a decoupling of taxonomic\n12\nand functional”) while looking at the figure. Please, either extend\nthe caption, or add some text to interpret this figure in the page 1 of\nthe supplementary data (because I’m sure you have to keep the\nmain text as concise as possible). I think that the statement made\nline 68-69 is very important (even if already known in others\ncontext/studies) and could be more easily understood if the\ndescription of Fig. S1 would be extended in supplementary files.\nThank you for pointing this out. We have expanded the caption of Fig.\nS1 to provide a more detailed explanation of its construction and\ninterpretation. We hope this makes the figure easier to understand.\n“Fig. S1. Taxonomic structure dynamics of large herbivores. Succession\nof taxonomic modules (M1 to M35) derived from network analysis,\nplotted against time (in million years, Ma). Each module represents a\ndistinct taxonomic structure based on species and continent-stage\n(nodes in the network) associations. Different colors and shapes indicate\ncontinent-stages across the different continents (America, Europe, Asia,\nand Africa; see legend). Lines represent module transitions within each\ncontinent, illustrating large-scale patterns of taxonomic succession.”\nL. 291-293: “The aggregation of localities by continent and stage\nyielded a total of 78 continent-stage assemblages, which are the\nanalytical basis for the study”. I think that’s totally needed and fair\nfor a study at the temporal and spatial scale, but could you develop\na bit the potential limitation of such approach in your method\nsection (or in supplementary)? What could be the impact of such\ndata agglomeration? For example, could different UFFs co-exist on\none continent (e.g. SE Asia vs Central Asia) and thus create fake\nstability by experiencing different trajectories of changes through\ntimes?\nWe used the continent-stage aggregation to correct for potential biases\nin spatial and temporal sampling. Prior to this aggregation, we carefully\nchecked each occurrence to refine the temporal range, correct locality\nnames, and verify taxonomic assignments. This work was carried out in\ncollaboration with expert taxonomists (coauthors of the study).\nThis allowed us to examine macro-scale patterns related to major abiotic\nchanges, without the confounding effects of uneven geographical\nsampling. We then aggregated the data by chronostratigraphic stages,\nwhich vary in temporal range. This approach narrows the temporal range\ntoward the present to address potential bias, especially the tendency to\npool more recent data.\nThe continent-stage aggregation yielded a total of 78 continent-stage\nassemblages. We chose this approach due to the lack of precise locality\ncoordinates for most of the data, which made it unfeasible to accurately\nmap localities within specific regions or countries (smaller scales). By\nfocusing on a continental scale, we were able to utilize a large number of\nlocalities and study their patterns of change in terms of functional\nstructure.\nGiven our definition of ungulate functional faunas (ungulate continental\nassemblages with shared distributions of functional types), it is not\npossible for more than one UFF to co-exist on the same continent (e.g.,\nSoutheast Asia vs. Central Asia). Our analysis focuses on the functional\nstructure of large herbivores at a continental scale, so the stability or\nreassembly observed is at that macro scale, not at the local or regional\nlevel. Therefore, we do not believe our results represent \"fake\" stability;\nrather, they indicate real stability or reassembly, but at the continental\nlevel (UFFs) across the study period.\nIn response to the reviewer’s suggestion, we have included a section in\nthe supplementary materials discussing the potential limitations of the\nspatio-temporal standardization and how it might influence the\ninterpretation of our results.\nI think that the authors could have spend more time to compare\nknown patterns of taxonomic diversity in mammals that took\ninterest in clustering temporal fauna (e.g. chronofauna in both\nNorth America and Europe/Asia/Africa, see for example Figueirido,\nGibert, Casanovas-Vilar) and look for breaks and tendencies but I\nunderstand that Nature format is very constrained. However, I think\nthat the excellent work they did with functional diversity could have\nbeen even better if it would have been better compared to other\ndimensions of mammalian fossil biodiversity.\nAs the reviewer pointed out, due to the limited format, we could not\ndiscuss all the biotic and abiotic events that took place during the\nCenozoic. We focused our discussion on those events that had the\ngreatest impact on functional structure and diversity dynamics,\nparticularly those recorded at continental and global scales. We agree\nthat regional events are of significant interest for the evolution of certain\nmammalian families. In Blanco et al. 2021 and Nascimento et al. 2024,\nwe addressed some of these regional patterns, particularly in the Iberian\nPeninsula, due to the regional scale of those analyses. However, the\nmacro-scale approach of our current analysis led us to focus on the\nbroader patterns described in the manuscript.\nL. 74-75 : “Since 10 Ma, large herbivore assemblages have\nexperienced remarkable stability”.\nOne remark: it could be stated here that Asian UFFs are acting\nslightly differently than the others continental fauna.\nOkay, we have slightly rephrased the sentence to incorporate the\nreviewer's comment:\n“Since 10 Ma, large herbivore assemblages have experienced\nremarkable stability beginning in Africa, followed by America and Europe\naround 7 Ma, and culminating in the later stability of Asia around 4.5 Ma\n(Fig. 1A).”\nOne question: when looking at Fig. 1B, we can see how the\naverage/median point of UFFs in NMDS is moving in the\nenvironmental space through time, but I struggle to understand\nhow Fig. 1A and Fig. 1B relate.\nFor example, when looking at the orange/yellow section for Africa (i.e. the last 10 myr), we see strong displacement in the\nenvironmental space for the african UFF, but it corresponds to statis (or a single mode, turquoise-green) in Fig. 1A. In the same way, the early Cenozoic asian UFF relate to the first mode (i.e. red) and last for ~35 myr, but we can see in Fig. 1B this\nUFF moving a lot through the environmental fauna (i.e. black to red\n15\nsection of the line). Can you extend the explanation in\nsupplementary data to understand the relationship between Fig. 1A\nand 1B. Fig. S19 could be used to explain this relationship with\nmore details.\nIn Figure 1B, we plotted the functional space of all the continent-stage\nassemblages (represented as dots). The purpose of this figure is to\nillustrate the different trajectories in terms of functional beta diversity,\ncalculated here as turnover (see Methods), that the continents followed\nover time, occupying different regions of the functional space.\nThe relationship between Fig. 1A and 1B lies in the similar trajectories\nfollowed by the continents when they formed a joint module, indicating\nthat they explored similar regions of the functional space. Therefore, in\nFig. 1B, the key point is the direction of the trajectory movement and the\nexploration of these shared regions of functional space, rather than the\nmagnitude of change within the space.\nTo clarify this relationship, we have added text to the Methods section to\nexplain the functional space and its connection to Figs. 1A.\nL. 97-98 : “This was primarly driven by the emergence of functional\ndominance of FTs in larger body size categories in the Americas, in\ncontract to the mid-sized herbivores that characterize the UFF of\nEurasia (Fig. S17)”. Could it be a matter of the relationship between\nlandmass and (paleo)latitudes? The larger southern extent of Asia\ncould explain the survival of mid-sized herbivores and the\ncontinuous cooling of Cenozoic could act in very different direction\ndepending on the latitudinal distribution of landmass, isn’t it? The\nargument developed in the rest of this section (from L.95 to L.119)\ncould be completed by a discussion around landmass/biome (i.e.\nclosed vs open forests) distribution.\nWe agree with the reviewer that this is an interesting and relevant point.\nAccordingly, we have added the following text to address this comment\nL101:\n\"The geographical extent of landmasses could have significantly\ninfluenced the survival dynamics of different faunal groups. Specifically,\n16\nthe broader southern reach of the Eurasian landmass compared to North\nAmerica may have provided a more stable environment for mid-sized\nherbivores. This stability persisted despite the overarching trend of\nCenozoic cooling, which likely affected species differently depending on\ntheir latitudinal distribution.\"\nThe two references associated with EOC transition (29, 40) are\nassociated with very descriptive arguments (L.109-119) when\nexplanatory argument (i.e. mechanism) could/should be invoked\nhere (e.g. Latitude and landmass, Europe being restricted to high\nlatitude and forming an archipelago making climate tracking more\ndifficult for European fauna). On a more general statement, I think\nthat this paper is sometimes a bit too descriptive, more potential\nmechanisms could be invoked to explain the observed patterns.\nWe appreciate the reviewer’s suggestion to incorporate more\nmechanistic explanations rather than relying solely on descriptive\narguments. To address this, we have expanded our discussion of the\nEOC transition by explicitly considering the role of latitude, landmass\nconfiguration, and geographic constraints on faunal dynamics.\nSpecifically, we discuss how Europe’s restricted high-latitude position\nand archipelagic nature may have limited the ability of its fauna to track\nclimatic shifts compared to other continents. We have added the\nfollowing text to the main manuscript to strengthen this section in line\n126:\n“The biogeographic and climatic context of the EOC transition also\nplayed a crucial role in shaping the observed functional shifts (29).\nEurope’s restricted high-latitude position and its configuration as an\narchipelago during this period may have constrained faunal dispersal\nand climate tracking, making species more vulnerable to extinction as\nconditions changed. In contrast, larger and more continuous\nlandmasses, such as those in North America and Asia, provided broader\nclimatic gradients and refugia, which facilitated species persistence and\nfunctional stability (Fig. 1A). These geographic constraints, coupled with\nthe rapid cooling and heightened seasonality of the EOC transition (29),\ncontribute to explaining the pronounced shifts observed in European\nfaunal assemblages. ”\nreduction of the functional distance between Africa and Eurasia\nafter 21 Ma”. Very interesting, could you just add the age\nuncertainties (i.e. the corresponding geological stage length if I\ncorrectly understood your methodology) in order to make this\nstatement more subtle. The exact date of the main part of the\ninterchange/faunal mixing between Africa and Eurasia remain\ndebated and your sentence could imply that your study can\nprecisely date this event.\nWe thank the reviewer for pointing this out. We agree that this is a\nvaluable observation. To refine our statement, we used the break in the\nsegmented regression from the linear model (25.4249 Ma), which marks\nthe start of the reduction in functional distance between Africa and\nEurope. The reduction between Asia and Africa began later, during the\nAquitanian (21.735 Ma). Therefore, we have emphasized the beginning\nof the reduction in functional distance between Africa and Eurasia.\nWe have added the corresponding geological stage length and refined\nthe ages as follows in line 139:\n“The functional coalescence is reflected in the drastic reduction of the\nfunctional distance between Africa and Eurasia, beginning around 25\nMa, during the Chantian stage (27.82–23.03 Ma) in Europe, and later\nfollowed by the Aquitanian (23.03-20.44 Ma) in Asia (Fig. S3)”\nL. 298: “We excluded the species with no information about these\ntraits”. Did you analyse/reanalyse/measure species traits or did you\nfind the 13 traits being already described for the 3012 species\nincluded in your study?\nOur procedure for assigning traits to species was rigorous and involved\nmultiple sources and expert validation. First, we downloaded all available\ninformation from the NOW database, the PBDB, and selected\nliterature.Each species was then individually assessed by taxonomy\nspecists who are experts in their respective groups and coauthors of\nthis manuscript. Traits were also inferred from bone measurements or\nextrapolated from characteristics of closely related species when\nnecessary. Thus, these experts confirmed the validity of the trait\n18\nassignments. Certain groups, such as Proboscideans, were nearly\ncomplete, with most data sourced from a previous work by\nCantalapiedra et al. (2021). However, other groups posed challenges\ndue to the scarcity of published information (e.g., Anthracoteriids).\nOverall, the trait assignment process was a collaborative effort involving\n10 specialists over more than 2 years to complete the database.\nL. 306: “and assigned a weight of 1/6 to the remaining 12 traits”.\nHave you tried different weight, in other word, did you perform\nsensitivity analysis on this parameter ?\nWe applied different weights to the traits when calculating the distance\nmatrix using the R function daisy:\ndaisy(no_na_fdata_art, metric=\"gower\", type =\nlist(ordratio = ordratio),weights = weights)\nThis distance matrix was then used to perform a PCOA, followed by k-\nmeans clustering to define the functional types. We assigned higher\nweights to body size and tooth shape because these traits play a\nsignificant role in defining the ecological function of species. Without\nweighting, their effects would be dispersed across the other 12 dental\ntraits.\nTo balance the combined effect of body size and tooth shape with the\nremaining 12 dental traits, we used a weight ratio of 1:1/6. This ensures\nthat the influence of the 12 dental traits is roughly equivalent to the\ncombined effect of body size and tooth shape. Altering the weight values\nwould not change the key aspect of this approach, which is that the\ncombined effect of the two sets of traits (i.e., 3 for body size and tooth\nshape and 3/6 for the remaining 12 traits) remains consistent.\nThis weighting strategy follows the methodology outlined in\nCantalapiedra et al. (2021).\nL. 393-397: Thank for the clear interpretation of the Shannon index\nfor FD. Would you say that high Shannon values in your study\nrelate to strong ecosystem stability (high diversity + high\nredundancy)? The comparison of Shannon FD and Fig. 1 is not\ntrivial.\nthe Shannon index around 10 Ma (Fig. 2), which coincides with a\ntransitional phase in both functional diversity and functional structure\n(Fig. 1 and 2). After this point, we identified a break in the trend of the\ncurve through segmented regression, followed by a sustained drop in\nfunctional diversity. This shift coincided with a tipping point in the\nfunctional structure, marking the onset of a global reassembly process.\nIt is true that between 21 and 10 Ma, we observed a continuous increase\nin functional diversity (Fig. 2), which coincided with a period of functional\nstability in terms of functional structure (Fig. 1A). However, this\nsustained increase in functional diversity began earlier during the\nCenozoic, when functional reassembly occurred across all continents\n(Fig. 1A).\nMoreover, after the reassembly that began around 10 Ma, the system\nregained stability, which persisted until the present day, overcoming a\ncontinuous decline in functional diversity (Fig. 2).\nWe acknowledge that the interpretation of Shannon diversity in relation\nto ecosystem stability is complex. While high functional diversity values\nare often associated with greater ecosystem resilience and stability, our\nstudy suggests that moments of high diversity do not always correspond\nto long-term ecosystem stability. Instead, the dynamics are influenced by\nseveral factors, including the broader context of functional structure. For\ninstance, the period of high diversity observed between 21 and 10 Ma\nwas a time of functional reassembly, which was a process that unfolded\nacross continents (Fig. 1A). This highlights the notion that ecosystem\nstability is not only a result of high diversity, but also of the resilience and\nstructure of the functional groups within the system. In fact, stability can\nbe context-dependent, influenced by factors such as tectonics, climate\nchange, and other abiotic processes. Therefore, while high functional\ndiversity may suggest a resilient ecosystem, it is not always directly\nlinked to sustained stability, as the functional structure and the system’s\nability to withstand or recover from disturbances also play crucial roles.\n20"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors find major shifts in the functional diversity of ungulates over the 60 million years associated with ecological and environmental shifts. They also find a decoupling between taxonomic and functional diversity fluctuations through time. The research and methodology are sound and appropriate for the questions being asked by the authors. I continue to believe that this manuscript provides important findings about the dynamics of mammals over evolutionary timescales, providing a baseline for pre-human ecology. Moreover, I believe it will be of interest to a wide audience.\nAfter reviewing the minor revisions and responses to the reviewers, I have no further concerns or comments regarding the manuscript and would like to see this published in Nature Communications.\n(Remarks on code availability)\nThe trait and occurrence data is available for analysis and review. All necessary details are included and easy to interpret. The code is clean and easy to follow regarding reproducibility. Therefore, I feel that the provided code and data is acceptable for publication.", "reviewer_3": "Dear authors,\nThank you for taking into account all the suggestions, and modifying your text accordingly."}, "author_response": ""}]

Decision letter