Body mass evolution in extinct mammals is an important factor in understanding their diversification and niche occupation dynamics over geological time. To know how Macraucheniidae, South American native ungulates, diversified over 30 million years, we estimated body mass for 15 genera within this family and investigated which evolutionary process best explains body mass evolution in this lineage. For this purpose, we applied a set of body mass estimation equations and conducted comparative phylogenetic analyses to assess the presence of phylogenetic signals and identify the best-fitting model of trait evolution. Our results show that all Macraucheniidae genera can be classified as megafauna (i.e., mammals over 45 kg). The evolutionary model that best explains body mass diversification in this family is a directional, gradual increase over time. The subfamilies Cramaucheniinae and Macraucheniinae did not differ significantly in body mass evolutionary rates; however, the negative relationship between body mass and temperature, together with the most parsimonious evolutionary scenario, supports temperature as the primary ecological driver of body mass evolution in macraucheniids.      

Institutional abbreviations: MACN – Museo Argentino de Ciencias Naturais ‘Bernardino Rivadavia’, Buenos Aires, Argentina; MPEF – Museo Paleontológico Egídio Feruglio, Trelew, Argentina; MCL – PUC/MG – Museu de Ciências Naturais da Pontifícia Universidade Católica, Belo Horizonte, Brazil; MLP-PV – División Paleontología de Vertebrados, Museo de La Plata, La Plata, Argentina.

We sampled 86 specimens from 15 genera, representing all genera of the Macraucheniidae family (Supplementary Data 1). The dataset consists of linear measurements of skeletal elements. Part of these measurements were taken by the author (LSL) at MACN, MCL, MLP-PV, and MPEF, while the remaining measurements were compiled from the literature: Mercerat (1891); Ameghino (1904); Scott (1910); Kraglievich and Parodi (1931); Soria (1981); Cifelli and Soria (1983b); Soria and Hoffstetter (1985); Fariña et al. (1998); Schmidt and Ferrero (2014); Cheme-Arriaga et al. (2016); McGrath et al. (2018). It is important to acknowledge the inherent incompleteness of the fossil record, which results in heterogeneous anatomical representation among genera. For instance, the skulls and most postcranial elements of Coniopternium and Cullinia remain unknown, whereas genera such as Theosodon and Macrauchenia are represented by nearly complete skeletons.  Consequently, the dataset analyzed here reflects an uneven distribution of predictive information across the family.

We selected univariable predictive equations to include a greater number of predictive models and to ensure that all genera of the family were incorporated in this contribution, given the uneven degree of anatomical knowledge among the representatives of the group described above. Head–body length was not included among the predictors used in this study, although it is a potentially informative proxy for well-preserved taxa and a useful complement to other body-mass estimators (Nelson et al. 2023); however, we did not have the opportunity to obtain such measurements for the present dataset. A total of 49 predictive equations from Damuth (1990), Janis (1990), and Scott (1990) were selected. These equations were grouped into three anatomical categories: cranial (10 equations), dental (24 equations), and postcranial (15 equations) (see their descriptions in Supplementary Data 2). Although these predictive equations were calibrated on specimen-level datasets, they were applied here to homologous anatomical measurements rather than taxonomic units. Accordingly, all equations were applied to individual specimens, and genus-level body mass values were obtained by aggregating multiple independent estimates using robust central tendency metrics. This approach reduces the impact of outliers and uneven sampling, which are common in fossil datasets, while ensuring that genera are treated as statistical summaries of specimen-based predictions rather than as direct inputs to the regressions. Consequently, the models are applied within a comparable anatomical framework and remain appropriate for taxa with heterogeneous and incomplete anatomical representation. Although direct validation of prediction ranges is not possible for fossil taxa lacking extant analogues, the use of multiple independent estimators helps mitigate potential departures from the original calibration space. Moreover, because no living relatives are available for direct comparison, the selected predictive equations were applied using the broad ungulate calibration, except for equation F3, which required the equid calibration because the third trochanter is absent in Artiodactyla (O’Leary and Gatesy 2008).

Owing to heterogeneity in the number of equations and specimens across genera, as well as the broad ranges observed in some taxa, the median was selected for the comparative phylogenetic analyses, as it is less sensitive to outliers. This approach has also been adopted in other studies (e.g., Rothier et al. 2025). In addition, the arithmetic mean and a weighted mean were calculated to provide alternative central estimates of genus-specific body mass. The weighted mean was defined to assign greater weight to equations supported by a higher number of observations and lower variability, while reducing the impact of equations based on a few data points and showing high dispersion. To achieve this, the square root of the number of observations (combining taxon and equation) was used, thereby reducing the dominance of overrepresented equations. Additionally, the equation range plus 0.01 was incorporated to penalize equations with high dispersion, while still allowing the inclusion of equations represented by a single estimate. The final formula for the weighted mean is: Final weight = sqrt(n_obs) / (range_eq + 0.01), where sqrt(n_obs) is the square root of the number of observations, and range_eq is the range of the equation. Finally, a box plot was generated for data visualization. All analyses were performed in R version 4.5.2 (R Core Team**** 2025) using the following packages: dplyr version 1.2.1 (Wickham et al. 2025), ggplot2 version 4.0.3 (Wickham 2016), and viridis version 0.6.5 (Garnier et al. 2024).

Phylogenetic hypotheses and temporal occurrence intervals (FAD–LAD) for macraucheniid genera were obtained from Lobo et al. (2024) using the ape R package version 5.8.1 (Paradis and Schliep 2019). The phylogeny adopted here corresponds to the hypothesis of Lobo et al. (2024), based on an implied-weights analysis with a concavity constant k = 6. Terminal taxa were time-calibrated using uniform distributions between their first (FAD) and last (LAD) appearance dates, and node ages were estimated using the paleotree R package version 3.4.7 (Bapst 2012) to account for stratigraphic uncertainty. Paleotemperature data were obtained from Hansen et al. (2013), a database providing global temperature estimates (°C) and their corresponding reliability weights for Cenozoic intervals. For each genus, a weighted mean temperature was calculated using Hansen et al. (2013) weight values, restricted to data within the taxon’s temporal range (between its FAD and LAD). This approach approximates the average climatic conditions experienced by each taxon while reducing the impact of anomalous temperature spikes, as isolated extreme values contribute little when assigned low weights. Moreover, the broad temporal windows of older taxa also further smoothed short-term fluctuations.

Phylogenetic signal in body mass was evaluated by estimating Pagel’s λ using the phytools R package version 2.5.2 (Revell 2023). To identify the best-fitting evolutionary model for body mass evolution in Macraucheniidae, we compared five models: pure Brownian Motion (BM; tested via λ), Early Burst (EB), Kappa, Brownian Motion with directional trend (BMT), and Ornstein-Uhlenbeck (OU) using the geiger R package version 2.0.11 (Pennell et al. 2014).  Model selection was based on AIC/AICw values. Ancestral body mass states were subsequently reconstructed with anc.trend (BMT model, phytools) to produce a phenogram illustrating node reconstructions alongside terminal taxa across the geological time.

Evolutionary rates of body mass were estimated using RRphylo package 3.0.2 (Castiglione et al. 2018) using phylogenetic ridge regression, and the significance of accelerations and decelerations in rates was assessed with search.shift. Branch-specific rates were extracted using standard (RRphylo$rates) and shift-adjusted regimes (RRphylo$multiple.rates), when significant shifts occurred. To evaluate differences in evolutionary dynamics between Cramaucheniinae and Macraucheniinae, a permutation test (999 iterations) was performed to assess whether they exhibit distinct body mass regimes.

Given that the directional evolution model provided the best fit according to AIC/AICw criteria, we tested whether changes in body mass were correlated with environmental variation, specifically climate change over the past 30 million years. Ordinary Least Squares (OLS) and Phylogenetic Generalized Least Squares (PGLS) regressions were performed, with temperature (°C) as the independent variable and log-transformed body mass as the dependent variable. Two PGLS models were fitted by varying the λ parameter: (1) PGLS with an estimated λ and (2) PGLS with λ = 1, representing a pure Brownian motion model. Model selection was performed using the sample-size-corrected Akaike Information Criterion (AICc), and model comparison assessed whether the mass-temperature relationship is better explained within a phylogenetic framework (PGLS) or a non-phylogenetic framework (OLS), using the caper R package version 1.0.4 (Orme et al.** 2023).