Adaptive shifts in Phanaeini dung beetles of the Mexican plateau cenocron in the Mexican transition zone

The Mexican Transition Zone is a biogeographically complex area where old and new lineages of Neotropical and Nearctic affinities overlap. Its biota was assembled by successive dispersal events of cenocrons, which are sets of taxa that dispersed during a given time interval from both North and South America and then diversified in the area. The Mexican Plateau cenocron, with Neotropical affinities, is found in temperate and dry climates in the Nearctic region. We hypothesised that it underwent an adaptive shift in environmental niche. We tested this hypothesis using a phylogenetic comparative framework, measuring phylogenetic signal and fitting to single optima macroevolutionary models, and an Ornstein‐Uhlenbeck macroevolutionary model with multiple optima. We used phylogenetic and distributional information of the tribe Phanaeini to assess whether there exists a distinction in conservatism between the earliest (Mexican Plateau) and most recent (Typical Neotropical) cenocrons within the Mexican Transition Zone (MTZ) as this tribe stands as a classic example of the dispersal and diversification patterns of cenocrons originating in the Neotropics. We identified different shifts in environmental requirements that match the niche description of the Mexican Plateau cenocron, suggesting that it was established through multiple adaptive shifts in the Mexican Transition Zone.

northern range limits often found in the tropical lowlands of Mexico (Wiens & Donoghue, 2004).Many species, however, could overcome such barriers and adapted to novel environments given enough geological time, originating complex biogeographical patterns with vicariant species and high endemism.
The Mexican Transition Zone (MTZ) exemplifies this adaptive phenomenon.Located within the mountain ranges of Mexico, Guatemala, Honduras, and El Salvador (De Mendonça & Ebach, 2020;Morrone, 2020), lineages of Neotropical and Nearctic affinities converge in this area under environmental conditions that differ from their ancestral niches.Biotas of the MTZ have been exposed to several climate fluctuations that promoted allopatric and parapatric divergence (Mastretta-Yanes et al., 2015), promoting speciation, high endemism and a clear elevational zonation that reflects a gradual temperature change (Espinosa-Organista et al., 2008;Morrone, 2017).According to the original proposals of the MTZ biotic assembly, the identity of the species in this area is given by their ancestral niche (temperate vs. tropical), the dispersal time, and their response to environmental instability (tracking vs. adaptation) (Halffter & Morrone, 2017;Lobo, 2007;Morrone, 2020).This hypothesis has been used to explain the fact that some of the endemic species found in high elevations and cold environments of the MTZ have Neotropical affinities and are closely related to tropical lowland species, contradicting what could be expected from the niche conservatism hypothesis.
An example of environmental adaptive shift corresponds to the species in the Mexican Plateau, which are mainly distributed in the arid and semiarid ecosystems regardless of their tropical origin (Halffter & Morrone, 2017).According to the biotic assembly hypothesis, these species represent the oldest Neotropical dispersal event into the MTZ.They might have arrived from South America as early as the Cretaceous-Palaeocene and then diversified in the Mexican Plateau specialising to temperate and xeric conditions (Halffter, 1962(Halffter, , 1964;;Halffter & Morrone, 2017), though recent biogeographical analyses suggest a later dispersal during the Oligocene-Miocene (Gillett & Toussaint, 2020;Halffter, 2022;Moctezuma & Espinosa de los Monteros, 2024).Such an early dispersal might have given enough time for local speciation since the Mexican Plateau emerged from the sea 66 mya (Fernandes et al., 2019) and before it lifted to its current elevation 44 mya (Fernandes et al., 2019).During this time, the area faced climatic and topographical changes that formed a sky-island archipelago where populations persisted despite severe environmental changes, like the expansion of arid scrublands (Mastretta-Yanes et al., 2015).These conditions might have prevented species from tracking their ancestral niche and promoted adaptation to novel environments.
A cenocron is a set of taxa that share a geographic origin, that dispersed during a given time interval and diversified in a particular area (Fernandes et al., 2019;Morrone, 2020).The age of these dispersal events has been established previously using the taxon's distribution size, niche width, phylogenetic affinities, and differential species richness between regions (Juárez- Barrera et al., 2020), under the assumption that phylogenetic inertia controls the spatial distribution of the biota (Halffter & Morrone, 2017;Morrone, 2020).Consequently, the distinctive biota of the MTZ has been assembled by the successive dispersal events of cenocrons that arrived from both North and South America (Halffter & Morrone, 2017), the oldest being the Mexican Plateau cenocron.Subsequent cenocrons are the Mountain Mesoamerican (Oligocene-Miocene), which dispersed from the tropical lowland forests to mountain regions in Central America and is associated with pine-oak and cloud forests (Halffter & Morrone, 2017).Then, the Nearctic cenocron dispersed from North America through the mountains of the MTZ (Miocene-Pliocene, Morrone, 2020) and is found in pineoak forests and mountain grasslands.Finally, the Typical Neotropical cenocron dispersed in the Pleistocene from South America and is restricted to warm lowlands with cloud, tropical rain, and tropical dry forests.Typical Neotropical species are distributed in more continuous and stable conditions that allowed them to track their niches, thus exhibiting a higher degree of niche conservatism.
To evaluate if there is a difference in niche conservatism between the oldest (Mexican Plateau) and the most recent (Typical Neotropical) cenocrons that assembled in the MTZ, we use dung beetles, members of the tribe Phanaeini (Coleoptera: Scarabaeidae).This tribe has served as a classic example to illustrate the dispersal and diversification of the cenocrons of Neotropical origin (Halffter & Morrone, 2017).The tribe Phanaeini is endemic to the Americas (Edmonds, 1994) and has at least 216 species classified in 11 genera (Cupello et al., 2023).Recognised as one of the best studied groups of scarab beetles (Cupello et al., 2023;Gillett & Toussaint, 2020), the tribe Phanaeini stands out as one of the few scarab tribes in which phylogenetic relationships have been explored multiple times (Gillett & Toussaint, 2020;Maldaner et al., 2020;Philips et al., 2004;Price, 2007Price, , 2009)).The most complete phylogeny available to date (Gillett & Toussaint, 2020) includes 74 species and incorporates macroevolutionary and biogeographical analyses.Although several taxonomical changes have led to an unexpected increase in the number of species in recent years (Cupello et al., 2023;Lizardo et al., 2022;Moctezuma et al., 2021;Moctezuma & Halffter, 2021), Gillett & Toussaint's (2020) dated phylogeny remains a valuable resource that allows to apply phylogenetic comparative methods to test evolutionary hypotheses.As a part of the Mexican Plateau cenocron, phanaeines might have dispersed early in the Oligocene-Miocene (Gillett & Toussaint, 2020;Price, 2009), then the genus Phanaeus Macleay dispersed and diversified, where most of its species are found (Lizardo et al., 2022).Such species have an important degree of speciation unique to the MTZ (Halffter & Morrone, 2017) and exhibit a vicariant pattern due to allopatry (Edmonds, 1994) or environmental conditions (Moctezuma et al., 2021).In contrast, the species in the Typical Neotropical cenocron of this group dispersed more recently, have a continuous distribution, lack notable vicariance, and are restricted to tropical lowland forests (Halffter & Morrone, 2017).
This study aims to examine the extent of niche conservatism and adaptive shifts among Phanaeini dung beetles in the MTZ.Our primary goal is to discern differences in niche conservatism between two cenocrons within the MTZ: the oldest cenocron, the Mexican Plateau, and the most recent cenocron, the Typical Neotropical.We predict that due to their early hypothesised dispersal into the MTZ and the environmental instability resulting from topographic changes, the Mexican Plateau cenocron species will display noticeable adaptive shifts and reduced niche conservatism.In contrast, Typical Neotropical cenocron species, which hypothetically dispersed more recently and have experienced continuous and stable environments, will show a higher degree of niche conservatism and no adaptive regimes will be found.To test these hypotheses, we applied a phylogenetic comparative framework and an Ornstein-Uhlenbeck macroevolutionary model with multiple optima to identify adaptive regimes in two niche traits (precipitation and temperature).We expected to find adaptive shifts towards colder and drier environments in species found in the Mexican Plateau.These shifts can be interpreted as niche convergences due to multiple adaptive shifts and can be used to identify the Mexican Plateau cenocron species in a phylogeny.

| Data sources
We used the most complete dated phylogeny at the species level available for the tribe (Gillett & Toussaint, 2020), which is a maximum clade credibility chronogram derived from the best BEAST analysis based on two uncorrelated lognormal relaxed clocks and a Yule tree model.From this phylogeny, we selected clade IV, which corresponds to the Phanaeini excluding Gromphas Brullé, often referred to as the subtribe Phanaeina (Cupello et al., 2023).We excluded from the analysis two species: Coprophanaeus machadoi (Pereira & d'Andretta) and Megatharsis buckleyi Waterhouse.The former was initially considered a synonym of C. saphirinus (Sturm), leading to the merging of their records as a single species in the GBIF backbone.It is important to note that Coprophanaeus machadoi was revalidated by Cupello et al. (2023); however, the GBIF dataset has not been updated to reflect this nomenclatural change, resulting in the continued merging of their records and it is not possible nor advisable to separate the records by location alone.On the other hand, M. buckleyi has no georeferenced records in GBIF.Although there are references regarding the distribution of both species (Gillett et al., 2009;Maldaner et al., 2020), the occurrences presented therein have many spatial issues and it is not possible to characterize their niche.This resulted in a tree with 69 species (see supplementary material), from which we obtained occurrence records from GBIF (2019) and relevant literature sources (Edmonds & Zidek, 2010;Gillett et al., 2009;Lizardo et al., 2022;Maldaner et al., 2020).To improve data accuracy, we performed an automatised data cleaning procedure using the CoordinateCleaner R package (Zizka et al., 2019) and a manual cleaning based on specialised literature (Edmonds & Zidek, 2010;Lizardo et al., 2022).The final dataset contained 5367 records; some species hold only three records while others have more than 500 (Table 1).
We used these records to extract environmental data from WordlClim's V2 (Fick & Hijmans, 2017) at a resolution of 30 arc sec.The environmental variables analysed were mean annual temperature (bio1) and annual precipitation (bio12), from which we obtained the mean, standard deviation, and standard error for each species.This was done using the R packages terra (Hijmans et al., 2022) and tidyverse (Wickham et al., 2019).Data are available on the GitHub repository (https:// github.com/ Virid ianaL izardo/ phana eini_ niche_ shifts).

| Models for continuous character evolution
As suggested by Wiens et al. (2010) and applied by subsequent authors (Hawkins et al., 2014;López-Estrada et al., 2019;Morinière et al., 2016), we compared the fit of three models: White Noise (WN), Brownian Motion (BM), and Ornstein Uhlenbeck (OU).WN is a model without phylogenetic signal, which assumes that data come from a single normal distribution with no covariance structure among species (Pennell et al., 2014).The fit of a trait to this model implies the absence of niche conservatism.Conversely, if a trait fits a BM or OU model, it will indicate that the traits have a distribution related to the phylogenetic structure.BM assumes that the correlation structure among traits is proportional to the shared ancestry (Pennell et al., 2014).In this model, species inherit their niches from their ancestors and slowly diverge by drift, thus indicating niche conservatism (Cooper et al., 2011).OU is a modification of the BM model (Cooper et al., 2016), that incorporates a central tendency (θ) and an attraction strength proportional to the parameter α (Pennell et al., 2014).This model can also indicate niche conservatism (Wiens et al., 2010).We used the R package geiger (Pennell et al., 2014) to assess the fit of the traits to these models.Then, we compared the goodness of fit using Akaike Weights to select the best fitting evolutionary model using phytools R package (Revell, 2012).

| Phylogenetic signal
To complement the previous analyses, we evaluated phylogenetic niche conservatism by calculating Pagel's λ (which assumes Brownian motion) with adephylo (Jombart et al., 2010), which shows good behaviour when working with small trees (Münkemüller et al., 2012).

| Identifying adaptive regimes
To test if there are adaptive shifts along the phylogeny, we applied a multi-optima Ornstein-Uhlenbeck (OU) model using the R package bayou (Uyeda & Harmon, 2014), based on a Bayesian framework.This analysis identifies changes in θ which are then identified in the phylogeny.The θ represents the optimal phenotypic value followed by an independent adaptive regime.This method applies a reversible-jump Markov chain Monte Carlo (MCMC) algorithm bayou (Uyeda & Harmon, 2014) to identify such shifts by running a one-million generations chain for each trait, using weakly informative priors (half Cauchy for α and σ 2 , mean and 1.5 standard deviations of observed data for θ, and the standard error was assigned based on the observed data), from which 30% of the generations were excluded as burn-in.Then, we detected adaptive regimes using a 0.2 posterior possibility cut-off threshold.Finally, to identify the distribution of species with adaptive shifts, we generated two richness maps.One map represented species with changes in θ, while the other depicted species without such changes.This was done by counting the number of species within a hexagonal grid of 1.5° using recorded species presence data obtained from GBIF (2019) and relevant literature sources (Edmonds & Zidek, 2010;Lizardo et al., 2022).

| RESULTS
Precipitation fits an OU model (Figure 1) with α = .04.This trait has a statistically significant and high phylogenetic signal (λ = .79,p < .05).The multi-optima OU model did not find adaptive regimes in the precipitation but indicated a rate of adaptation (α) of 1.33 (SD = 10.6) and a mean evolutionary rate (σ 2 ) of 0.09 (SD = 0.24).The precipitation value at the root was 1906 mm, 95% CI [509,7510].Precipitation has 5.19 million years of adaptive half-life (the amount of time expected for a trait to evolve halfway from its ancestral value to its optimum, calculated as ln(2)/α).
Temperature fits a white noise model (Figure 2a,b).It has low and statistically non-significant values of phylogenetic signal (λ = .20,p = .3).The OU multi-optima model found a rate of adaptation (α) of .37 (SD = 0.6), an evolutionary rate (σ 2 ) of 4.51 (SD = 0.07), and an adaptive half-time of 1.88 million years.The temperature value at the root was 22.71°C, 95% CI [21.76,23.60].The model detected three adaptive regimes found in the colder end of the trait distribution (15°C to 21°C).
The regimes with change in theta values are (1) Phanaeus vindex MacLeay, P. igneus MacLeay, and P. difformis LeConte (θ = 19.29°C);(2) P. texensis Edmonds and P. quadridens (Say) (θ = 18.44°C); and (3) P. yecoraensis Edmonds (θ = 18.39°C).The branches with shifts in the thermal optima have a maximum crown age around 21 Ma, calculated as the maximum crowning age of the species with niche shift.Species with no shift in thermal optima are distributed in the Neotropics in high richness areas (Figure 2c), while species with niche shift are in the Mexican Plateau northwards to the USA East coast, with almost no overlap between species (Figure 2d).

| DISCUSSION
Our analysis of niche conservatism and adaptive shifts among Phanaeini dung beetles in the Mexican Transition Zone (MTZ) revealed differences between the Mexican Plateau and the Typical Neotropical cenocrons, and between the temperature and precipitation as niche characteristics.Multiple measurements showed the absence of niche conservatism in the thermal niche of the tribe Phanaeini which fits a non-phylogenetic model (WN) (Ferrari et al., 2012) and has a low and non-significant value of λ.The more complex multi-optima OU model, however, showed that the thermal niche does have a phylogenetic structure as a product of the convergence of three regimes in the Nearctic region, while the remaining species have a conserved niche towards warmer regions (~23°C).We interpret this as convergent evolution for the characteristic Mexican Plateau thermal conditions that have appeared three times.On the other hand, precipitation fits a phylogenetic model (OU) with a very low α value, which resembles a BM model.This trait also has a statistically significant and high phylogenetic signal (λ), suggesting that precipitation is a trait that has some degree of phylogenetic conservatism within the studied species, which is confirmed by the lack of adaptive regimes.It is possible that the seasonality of the precipitation could be more important in the dispersal dynamics and the establishment in novel environments, yet it remains to be evaluated.
F I G U R E 1 AIC weights of Phanaeini niches' annual precipitation and mean annual temperature values fitted to macroevolutionary models, the asterisk indicates a significant value (p < .05).
The pattern found here fits the niche predictions for the Mexican Plateau cenocron, where a Neotropical lineage shifted its niche to colder conditions.According to Halffter (1962Halffter ( , 1964) ) and Halffter and Morrone (2017), the ancestor of the Phanaeus species dispersed to Mesoamerica during the Late Cretaceoous-Eocene, contrastingly Gillett and Toussaint (2020) suggest it happened during the Oligocene.Then, the dispersal from Mesoamerica to the highlands is proposed to have occurred in the Miocene (Gillett & Toussaint, 2020).Notably, the differing proposed dispersal dates hint that the Mexican Plateau Cenocron's dispersal was likely influenced more by the Trans-Mexican Volcanic Belt orogeny during the Plio-Pleistocene than the Plateau's elevation.This dispersal may have occurred during the warmer climate of the Middle Miocene Transition, when the Mexican Plateau was covered by warm-temperate evergreen broadleaf and mixed forest, despite the increase in elevation.The mismatch could also be due to the updated topology of the phylogenies used to interpret the geographical patterns.The original works from Halffter (1962Halffter ( , 1964) ) relied heavily on subgenus separation, which is no longer valid  (Gillett & Toussaint, 2020;Lizardo et al., 2022;Moctezuma & Halffter, 2021;Price, 2009).
Other studies found a similar dispersal date.Halffter (2022) discovered that Boreocanthon Halffter, 1958 andCanthon imitator Brown, 1946, which are frequently cited as examples of the Mexican Plateau cenocron, originated during the Miocene.Moctezuma and Espinosa de los Monteros (2024) explored the biogeographic history of 22 Mexican Scarabaeinae species finding the same dates and suggesting that the Phanaeus clade alone represents the Mexican Plateau cenocron.All this indicates that the Miocene not only represents the period of Phanaeus' appearance but also that of Canthon, strengthening the establishment of the cenocron during these dates and suggests that the dispersal of this cenocron was not interrupted by the lifting of the Mexican Plateau.These findings complement our own, providing a more comprehensive perspective on the biogeographic history of the Mexican Plateau cenocron.
It is suggested that Phanaeus diversified and expanded northward into what today is the USA.These species in the USA are part of the Mexican Plateau cenocron and match those found in this study as niche shifters (Figure 2b), as well as their distribution (Figure 2d).Conversely, the species without niche shifts, which belong to the Typical Neotropical cenocron, are expected to occur along the coastal plains of the Gulf of Mexico to the east, and the Pacific Ocean to the south and southwest.This cenocron includes the Phanaeus species, along with some species of Coprophanaeus Olsoufieff and Sulcophanaeus Olsoufieff (Halffter & Morrone, 2017).
At the time of the predicted niche shift (21 mya), the topography and climate of Mexico were different from today.There were no mountains, but low-elevation volcanoes (Ferrari et al., 2012); then, from the late Miocene to the early Pliocene, plateaus covering large regions were created and the elevation range increased hundreds of metres over 1-2 myr.Afterwards, from the Pliocene to the Pleistocene, the Trans-Mexican Volcanic Belt (TMVB) lifted >2000 m above the plateau level (Mastretta-Yanes et al., 2015).At the same time, the climate changed from warmer conditions to a glaciation during what is called the Middle Miocene Climate Transition (Frigola et al., 2018).During these environmental fluctuations, the main response of the species was dispersal according to the niche conservatism hypothesis (Donoghue, 2008;Lobo, 2007).In this situation, the TMVB rose, making it impossible to disperse back.The lineages had to adapt to the new climate to avoid extinction, and they also had the chance to diversify in association with the Cenozoic mammal megafauna in the MTZ (Edmonds, 1994;Favila, 2012;Halffter & Morrone, 2017;Joaqui et al., 2019;Miller, 1983).In summary, we propose that it was not only the orographic change that prevented the dispersal, but the climatic fluctuations and the adaptative shifts caused by both factors.
Our results show that cenocrons are evolutionary biotic units affected by macroevolutionary forces.Nevertheless, a word of caution is necessary since the phylogenetic hypothesis that we used for the analysis does not include all the species and many taxonomic changes have occurred since its publication.Therefore, it is not possible to associate the species with the occurrence data with total accuracy.This has occurred mainly in Phanaeus, where some species have been split into several different ones (Lizardo et al., 2022;Moctezuma et al., 2021;Moctezuma & Halffter, 2021).Furthermore, the full variety of environmental conditions might not be completely sampled, and more fieldwork is still needed.Although the present work summarizes the available data on Phanaeini and their distribution, it is necessary to keep working on the phylogenetic reconstruction of the group and to fill the Linnean and Wallacean shortfalls.Despite these shortcomings, the analyses implemented herein were able to identify the Mexican Plateau cenocron by using the phylogenetic comparative method.We suggest validating the method proposed herein by applying similar methods to larger phylogenies and novel approaches to characterize the niche.
We found distinctions between the evolution of the niche of species of the Mexican Plateau and Typical Neotropical cenocrons, as well as differences in the temperature and precipitation traits.Our results align with predictions for the Mexican Plateau cenocron, where Neotropical lineages shifted their niches to adapt to colder conditions and converged into the Mexican Plateau region.Finally, we propose that cenocrons exhibit adaptive regimes to certain climatic variables and macroevolutionary models can unify classical biogeographic theories with evolutionary concepts.Here, we show that temperature is indeed a central driver in the assembly of the MTZ as previously suggested, and that niche conservatism is not the only process involved in the establishment of diversity gradients, but its interaction with parallel niche shifts that result in multiple lineages converging to the same novel climatic conditions.

F
I G U R E 2 Phanaeini's mean annual temperature value fitted to an OU multi-optimum model and its adaptive regimes.The posterior probability distribution of each regime is shown in (a), where the dotted lines represent the theta (θ) value of each regime.(b) shows the ancestral trait reconstruction with the adaptive regimes shown in the same colour as in (a).In (a, b), grey indicates species with no change in theta and the colours indicate the selective regimes detected with the analysis that have a different theta value.In the lower panel, the distribution of species richness is compared between those species having conserved niches (c) and those species with shifted adaptive regimes (d).