Phylogenetic α- and β-Diversities Jointly Reveal Leaf-Litter Ant Community Assembly Mechanisms Along a Tropical Elevational Gradient

Despite the long-standing interest in the organization of ant communities across elevational gradients, few studies have incorporated the evolutionary information to understand the historical processes that underlay such patterns. Through the evaluation of phylogenetic α and β-diversity, we analyzed the structure of leaf-litter ant communities along the Cofre de Perote mountain in Mexico and inferred its putative driving forces. Lowland and some highland sites showed phylogenetic clustering, whereas intermediate elevations and the highest site presented phylogenetic overdispersion. We infer that strong environmental constrains found at the bottom and the top elevations are favoring closely-related species to prevail at those elevations. Conversely, more benign conditions at intermediate elevations suggest interspecic interactions being more important in these environments. Total phylogenetic dissimilarity was driven by the turnover component, indicating that the turnover of ant species along the mountain is actually shifts of lineages adapted to particular locations resembling their ancestral niche. The greater phylogenetic dissimilarity between communities was related to greater temperature distances probably due to narrow thermal tolerances inherit to several ant lineages that evolved in more stable conditions. Our results suggest that the interplay between environmental ltering, interspecic competition and habitat specialization plays an important role in the assembly of leaf-litter ant communities along elevational gradients. together, these results highlight the role of species sorting processes 67 , where the phylogenetic composition is mainly driven by deterministic process (i.e. habitat ltering) in response to local environmental conditions rather than stochastic processes (i.e. dispersal limitation).The simultaneous examination of the phylogenetic α- and β- diversity (and its components of turnover and nestedness) enhances our understanding of the relative importance of deterministic and stochastic processes in structuring patterns of ant diversity. Here, we showed that environmental ltering, interspecic competition and habitat specialization jointly the leaf-litter ant communities along the Cofre de Perote. These results highlight the importance of deterministic (niche-based) processes over stochastic processes. Further, our results provide insights about phylogenetic niche conservatism since some ant lineages have retained the necessary harsher (the colder at the summit). Additionally, large each site total


Introduction
Elucidating the mechanisms underlying the natural variation in species richness and composition across local communities has challenged ecologists for more than a century 1 . Since the development of community phylogenetics, a vast number of studies have used the phylogenetic approach to disentangle the relative importance of both deterministic and stochastic processes involved in the coexistence of species within a community 2,3 . The cornerstone of the phylogenetic approach is the assumption of phylogenetic niche conservatism, i.e., the tendency of species to conserve their niches over evolutionary history, with closely related species being ecologically more similar to each other than to distantly related species 4 . Based on this assumption, phylogenetically clustered communities (i.e., co-occurrence of closely related species) are typically thought to result from environmental ltering, where taxa are ltered by constraints imposed by the environment 5 . However, modern coexistence theory considers that phylogenetic clustering could also be driven by competitive exclusion where entire clades have higher competitive abilities leading to an exclusion of other lineages 6 . Conversely, phylogenetically overdispersed communities (i.e., co-occurrence of distantly related species) are thought to be structured by competition, which tends to select species with different ecological traits and thus low niche overlap 7 .
Finally, a lack of phylogenetic structure suggests a predominance of neutral processes 8 or the counteraction of both habitat ltering and interspeci c competition 9 . Since these mechanisms can act simultaneously within communities 10 , the critical question is no longer which mechanism structure communities but rather which one plays a dominant role in community assembly along environmental gradients 11 .
PBD explicitly adds phylogenetic information to evaluate how evolutionary relationships of component lineages change between communities or biomes across space 12 . Furthermore, PBD links local processes (e.g. biotic interactions or environmental ltering) to more regional evolutionary processes (e.g. trait evolution, speciation and dispersal), hence providing further information about how current or historical environmental factors in uence the variation in species compositions of communities across space 12,15 .
Environmental ltering and dispersal limitation have been suggested to largely determine patterns of community composition 16,17 . Habitat ltering is expected to limit community members to habitats that resembles the ancestral niche where its lineage originated, i.e., habitat specialization 18 , resulting in closely-related species occupying similar portions of regional-scale climatic gradients. On the other hand, the dispersal limitation processes predicts that community dissimilarity between sites will be correlated to geographical distances separating those sites, regardless of the environmental gradients 16 . Since environmental ltering and dispersal limitation are not mutually exclusive, a greater variation in community dissimilarity predicted by geography re ects a greater importance of dispersal limitation, whereas greater variation explained by environmental distances indicates that habitat ltering is the structuring force 19 .
Recent methodological advances have improved our understanding of the origin and maintenance of geographic patterns of PBD through its decomposition into two antithetic components that account for the replacement of lineages between sites (the turnover-resultant component) and differences in the phylogenetic diversity between nested assemblages (the nestedness-resultant component; 20,21 ). Although the turnover and nestedness components both contribute to total dissimilarity, their relative importance depends on the processes structuring communities. For instance, if environmental conditions vary across space, and species are adapted to particular conditions, the turnover component of PBD is more likely to shape community composition under environmental ltering. Conversely, nestedness tends to be a more common component due to a number of processes as selective colonization, selective extinction, nestedness of habitats or passive sampling 21 . Analyzing the relationship between PBD (and its components) with climatic and geographic variables can provide a better understanding of community assembly along environmental gradients 22 .
Several studies have placed the phylogenetic α-and β-diversity into ecological and biogeographic contexts, such as the latitudinal pattern of phylogenetic diversity 23 , the phylogenetic modi cation of native communities under species invasions 24 and the community structure along elevational gradients [25][26][27] . Particularly, ecological gradients occurring in mountains have successfully served to test community assembly processes given the strong dependence of species elevational distributions to the environmental conditions found at particular elevations 28, 29 . Particularly for thermophilic taxa, such as ants, deterministic processes related with both temperature [30][31][32] and precipitation 26, 33 have been shown to in uence community membership at elevational gradients. In tropical mountains, as elevation increases, temperature tends to decrease whereas precipitation increases 34 and, as a result, the abiotic conditions (cold and wet) at high-elevation sites become physiologically stressful for a majority of ant species in comparison with relative benign environments (warm and humid) at low elevations. This spatial structure created by the interplay of temperature and precipitation has been used to explain the changing phylogenetic structure from overdispersed ant communities at low elevations towards clustered communities at higher elevations [25][26][27] . Despite these ndings, the relative importance of deterministic and stochastic processes driving ant community organization under the phylogenetic perspective is still incipient, particularly when tropical regions are considered.
In this study, we measured the phylogenetic α-and β-diversity to infer whether deterministic and stochastic processes are driving the organization of leaf-litter ant communities occurring along a tropical elevational gradient in Mexico. Leaf-litter ants have numerous attributes that make them an ideal system to explore assembly mechanisms. For instance, they exhibit high levels of local co-occurrence; with up to 35 species co-existing in only one square meter 35 . Further, ant distributions are highly constrained by local and regional climate 36,37 . They display a wide variety of both individual 38 and colony-level 39 thermoregulatory strategies to cope with cold and hot conditions. Lastly, molecular phylogenetic analyses of the major ant lineages are providing a stable framework to understand the evolutionary relationships of the group [40][41][42] . Hinged on the assumption that ant species have retained their traits along their evolutionary history (i.e., phylogenetic niche conservatism 23,43 , we expect that (1) ant communities inhabiting in more benign and stable habitats (warm and humid habitats) found at low elevations will show phylogenetic overdispersion since negative interspeci c interactions might be more intense in these environments. Conversely, cold and wet conditions at high elevations would lead to phylogenetically clustered communities considering that only closely related species of a subset of lineages possess the physiological traits that allow them to persist in the harsh conditions present at those elevations. In terms of PBD, it is expected that (2) lineages would have high habitat specialization to the conditions where they originated and thus drastic environmental changes along the elevational gradient would lead to high lineage replacement. This will result in a greater importance of turnover on total phylogenetic dissimilarity along the elevational gradient. Finally, we can expect that (3) pairwise dissimilarity values for PBD and their components would be mainly explained by environmental ltering (i.e., climatic distances), following their ancestral climatic a nities rather than dispersal limitation.

α-diversity
Overall, the tendency of SES.PD, SES.MPD and SES.MNTD was negative values at low and some high elevational sites, whereas positive values at mid-elevations and the highest site. Considering the SES.PD, this means that three intermediate-elevational sites (1000, 1500 and 3000) present high evolutionary history, whereas four sites (0, 500, 2000, 2500) have low phylogenetic diversity (Fig. 1a). This same trend was observed for the SES.MNTD (Fig. 1b), indicating that sites containing higher evolutionary history are related to more distantly species (overdispersion), whereas sites with lower evolutionary history are inhabited by more closely-related species (clustering). For the SES.MPD, overdispersion was observed at three mid-elevation sites (1000, 1500, 2000) and the highest site (3000), whilst the rest of the sites (0, 500 and 2500) presented a clustering pattern (Fig. 1c). All those trends were consistent across the 1000 phylogenetic trees and remained the same when the MCC tree was used (Fig. 1). All these α-diversity trends were supported by a variable number of trees being statistically different from a null expectation at α = 0.05, but most common at 2500 m for SES.PD and SES.MPD (Supplementary Table S4).
When we assessed the contribution of climatic variables on explaining the alpha metrics patterns, the additive effect of both PC1 Temperature and PC1 Precipitation explained a larger amount of variation in comparison with their individual effects ( Table 1; Supplementary Fig. S2). Moreover, the additive model was considered the best statistically supported model in explaining SES.PD, SES.MPD and SES.MNTD patterns. However, in some instances, the individual effect of PC1 Temperature and PC1 Precipitation resulted as equally probable models ( Table 1, Supplementary Fig. S4), suggesting that the main effect of both variables is biologically important, yet its combined effect re nes model performance. Further, the general trend for both linear and multiple regressions was a positive relationship (β > 0) between PC1 Temperature with all phylogenetic α-metrics, whereas PC1 Precipitation related in a negative way (β < 0; Table 1).

β-diversity
We found that total phylogenetic dissimilarity derived from multiple-site calculations exhibited considerable high values (PBD multi.sor =0.73±0.006). Decomposition demonstrated that the turnover component (PBD multi.sim =0.52±0.02) had a greater contribution to total dissimilarity compared with the relatively low values of the nestedness-resultant component (PBD multi.nes =0.21±0.01). These patterns remained when the MCC tree was used ( Supplementary Fig. S5). When assessing the patterns under the adjacent approach, we observed that the total dissimilarity (PBD adj.sor ) increased with elevation. For instance, lower values (<0.5) were observed at lowland sites (< 2000), whereas higher values (>0.5) at the highlands (>2000; Fig. 2a). In most cases, the β ratio calculation displayed values higher than 0.5 indicating that total dissimilarity is determined dominantly by turnover. A deviation of this general trend was observed at 1500-2000 sites where the nestedness (PBD adj.nes ) dominated over the turnover component ( Fig. 2b). Pairwise patterns of PBD showed that an increase of elevational distance results in greater total dissimilarity (PBD pair.sor ) and is produced by a mixture of turnover (PBD pair.sim ) and nestedness (PBD pair.nes ) values operating between different elevational ranges ( Supplementary Fig. S6).
According to GDM analysis, the full model explained a considerably high variation of PBD pairwise patterns. For instance, 79.2% of deviance was explained for total dissimilarity (PBD pair.sor ), followed by 49.9% for nestedness (PBD pair.nes ) and in minor instance the turnover (PBD pair.sim ) with 26.4%. All models presented statistical signi cance for total dissimilarity, whereas turnover and nestedness were poorly supported with only three and seventy-two (out of 1,000) statistically signi cant models respectively (Table 2). These results are consistent when the MCC tree was used in the analysis (Supplementary Table  S4). In all PBD components, the variance explained by the PC1 Temperature was much larger than explained by geographic distance or PC1 Precipitation , indicating that thermal a nity has played a much greater role than a nity to precipitation and dispersal limitation in shaping the phylogenetic composition of ant assemblages across our study area (Table 2). Table 1 Summary statistics (mean ±SD) of the coe cient of determination (R 2 ), bayesian information criterion (BIC) and beta coe cient (β) extracted from the set of linear models (1000 per alpha metric) adjusted independently against the standardized effect size of phylogenetic diversity (SES.PD), mean pairwise distance (SES.MPD) and mean nearest taxon distance (SES.MNTD) against the principal component of temperature (PC1 Temperature ), precipitation (PC1 Precipitation ) and the additive effect of both factors (PC1 Temp+ PC1 Prec ). Equally probable models were considered if the difference in BIC (ΔBIC) between the focal model and the model with the lowest BIC were < 2. Value obtained from the MCC tree is expressed in the bold number within parenthesis.

Discussion
Over the last two decades since Webb et al., 7 developed the baseline of community phylogenetics, several studies in this eld have been conducted to enhance our understanding of how deterministic and stochastic processes operate to assemble communities across space. In this study, we analyzed the phylogenetic α-and β-diversity to unravel the assembly of leaf-litter ant communities along a tropical mountain in Mexico. Overall, lowland and some highland sites contained communities with low evolutionary history and composed by closely-related species (phylogenetic clustering), whereas communities at intermediate elevations and at the highest site presented high evolutionary history and were composed by distantly-related species (phylogenetic overdispersion). All these patterns were highly supported not only across all the α-phylogenetic metrics (PD, MPD, MNTD) but also across 1000 hypothetical trees and the maximum clade credibility tree, thus patterns remained regardless of phylogenetic uncertainty. We found a positive relationship between all metrics and temperature, whereas a negative relationship with precipitation. Both the main and additive effect of precipitation with temperature were equally supported, yet the additive effect explained higher variation of phylogenetic αdiversity patterns for all metrics. Phylogenetic total dissimilarity was mainly driven by turnover, thus highland lineages are not subsets of lowland assemblages but instead they are communities with different lineage compositions. Finally, the temperature differences between sites emerged as the most important driver of total dissimilarity, turnover and nestedness components. In sum, our ndings provide additional evidence that evolutionary processes mediated by climate-related (deterministic) mechanisms are strongly involved in the assembly of ant communities along elevational gradients.
Ant communities at the lowland and highland sites (except for the 3000 site) displayed phylogenetic clustering. Such a pattern has been typically associated to the gradual loss of ant lineages due to abiotic constraints (i.e. environmental ltering) where only a subset of species possesses the necessary adaptations to persist under stressful conditions [25][26][27] . In our study area, the drastic reduction of temperature coupled with higher precipitation with increasing elevation may act as a strong lter on total diversity and the members of the regional species pool, such that the resultant high-elevation community would be phylogenetically clustered with taxa possessing wet and cold-resistant traits for dealing with this stress. Here, we found that >50% of species inhabiting phylogenetically clustered communities at high elevations belong mainly to two cold-specialist genus: Stenamma and Temnothorax (Supplementary Table S1). Strategies such as the hibernation of the ant colony during the winter season in Stenamma species 44 Table S1). The well-known effects of these ant species in disrupting and displacing native ants [49][50][51] suggest that patterns of phylogenetic clustering observed at lowland communities may not be driven by contemporary climate alone, but is also a result of the invasion of these ant species.
We found intermediate elevations constituted by phylogenetically overdispersed communities. Typically, phylogenetic overdispersion is interpreted as evidence of interspeci c competition since a long history of competitive interactions should cause evolutionary divergence in species niches 52 . Interspeci c competition among ant species is intense and often involved in the con guration of ant communities 53 , yet evidence suggest that the importance of competition may be higher in relatively benign, stable environments where abiotically stressful factors are absent 27,54 . Intermediate elevations at Cofre de Perote re ect these conditions considering that at these elevations temperature is not too low to freeze available water nor too high to evaporate it ( Supplementary Fig. S1). Besides, productivity, an ecological proxy of the amount of niches and resource heterogeneity in an ecosystem 55 , is expected to peak at this point since productivity is limited by drought at lowlands and cold temperatures at highlands. Thus, if competition is the driving mechanism at these intermediate elevations, we should observed communities containing a series of species with different evolutionary histories 56 . Indeed, clades are well represented in these communities with seventeen tribes (out of eighteen) containing a mixture of species from both tropical and temperate origins distributed at low and high elevations respectively.
Whilst competitive interaction is congruent as a structuring force in phylogenetically overdispersed communities at more benign habitats found at intermediate elevations, it is unlikely that this hypothesis stands for the isolated overdispersed community at the highest stressful elevation (i.e., 3000 m). Particularly, we observed that the ant community at this elevation was composed by eight species, each one belonging to different genera dispersed across the phylogeny. Some studies have posited that geographic isolation for historical climatic variations has played a key role in the distribution of species at high-elevation habitats 57 . Thus, the presence of species with contrasting evolutionary histories may suggest that communities at 3000 m could be acting as refugia, maintaining relict lineages that migrated from distant regions with temperate climate or were more widespread in the past but became geographically isolated as a consequence of habitat contractions in the last glaciations 58 .
The complementary use of the multiple-site, adjacent, and pairwise approaches, coupled with the decomposition of total dissimilarity into the turnover and nestedness components, signi cantly contributed to unveil the underlying mechanisms in uencing dissimilarity variation. For instance, all approaches showed a high dominance of the turnover component for total dissimilarity. This result indicates that the turnover of species is actually a turnover of entire lineages or clades 12 Table S1). Taken together, our results bring evidence that habitat specialization is not only key driver of compositional dissimilarity of species (i.e., 22,25,32 ), but also an important processes scaling up to entire lineages (i.e. phylogenetic niche conservatism), in such a way that evolutionary history strongly constraints the elevational distribution of ant species 25 .
Despite adjacent turnover component having a preponderant impact on total dissimilarity, we observed a breakpoint between 1500 and 2000 where the phylogenetic nestedness prevailed over turnover. The possible reason for this dominant nestedness in PBD is that selective extinction, operating through environmental ltering, is playing an important role in shaping the patterns of phylogenetic dissimilarity at this particular elevation 21 . Because environmental ltering favors certain traits over others, a high number of ant lineages are lost as a result of the temperature decrease and humidity increase occurring from 1500 to 2000 site 64 . Indeed, we noticed that the distribution of certain ant species was restricted to low and intermediate elevations (Supplementary Table S2). Particularly, several well-known neotropical species belonging to Brachymyrmex, Europhalothrix and Strumigenys genera completely disappear at 2000 site (or higher elevation sites; Supplementary Table S1). These ant genera are mainly restricted to more tropical environments found in lowland areas, with one or few species within the same clade capable to cross to high elevations 61 . Therefore, we speculate that this point of the mountain can serve to broadly separate two well-distinguished ant fauna: low-montane fauna, distributed from sea level to below the 2000 m of elevation, and high-montane fauna habiting at 2000 elevation and greater.
We show here that phylogenetic total dissimilarity was best explained by temperature distance among pairwise sites regardless of the precipitation or geographical distance among them. More speci cally, the further two elevational sites are in terms of their temperature, higher total dissimilarity between them is expected in comparison with two sites sharing similar temperatures. This result agrees with those found by Liu et al., 25 who documented that PBD in ant composition dissimilarity of the Hengduan mountain was mainly driven by climatic distance (in which a set of temperature-related variables were included). This high importance of temperature distance shaping PBD may be explained by the climate variability hypothesis proposed by Janzen 65 . The climate variability hypothesis proposes that species exposed to variable climates (such that occurring at higher elevations) evolve broad thermal tolerances, allowing those species to traverse climatic gradients found across elevations, resulting in a wider geographic distribution than thermal specialists from stable climates (generally found at lower elevations). Therefore, we should expect lineages which originated in more tropical regions (e.g., neotropical) to be characterized by narrower thermal tolerances resulting in more restricted distribution in comparison with lineages originated in more variable climates (e.g., nearctic). This consistent effect of thermal adaptations constraining not only species distributions 32,66 but those of entire lineages support that climatic niches are conserved over evolutionary history of the ant clade ( 25,43 ). Taken together, these results highlight the role of species sorting processes 67 , where the phylogenetic composition is mainly driven by deterministic process (i.e. habitat ltering) in response to local environmental conditions rather than stochastic processes (i.e. dispersal limitation).The simultaneous examination of the phylogenetic α-and β-diversity (and its components of turnover and nestedness) enhances our understanding of the relative importance of deterministic and stochastic processes in structuring patterns of ant diversity.
Here, we showed that environmental ltering, interspeci c competition and habitat specialization jointly structure the leaf-litter ant communities along the Cofre de Perote. These results highlight the importance of deterministic (niche-based) processes over stochastic processes. Further, our results provide insights about phylogenetic niche conservatism since some ant lineages have retained the necessary traits to colonize harsher environments (the colder habitats at the summit). Additionally, the large evolutionary history accumulated in the lineages inhabiting each elevational site along with the remarkable rates of phylogenetic turnover contributing to total phylogenetic dissimilarity con rms the importance of mountains not only as centers of species diversity but crucial reservoirs of unique evolutionary history 68,69 . Altogether, this work builds up on the theory that not only contemporary but historical factors also in uence the structure of leaf-litter ant assemblages along environmental gradients and this can be detected by integrating α-and β-phylogenetic diversities.

Study area and ant sampling
This study was conducted along the eastern slope of the Cofre de Perote mountain, in Veracruz, Mexico.
This region is located at the junction of the Trans-Mexican volcanic belt and the Sierra Madre Oriental. We selected eight study sites spanning an elevational gradient of 3500 meters of altitude (Fig. 3). Regardless of the geographical distance, all sites were systematically separated with an elevational difference of 500 meters on average between each other. We placed our study sites at the following elevations above sea level We obtained 320 m 2 leaf-litter samples in 8 study sites (see 32 for a complete description of sampling design). In each 1-m 2 quadrat, we collected the leaf litter inside and sifted it through a coarse mesh screen of 1-cm grid size to remove the largest fragments and concentrate the ne litter. The concentrated ne litter from each sample was suspended in independent mini-Winkler sacks for 3 days in the laboratory. Falling arthropods were collected into a container with 95% ethanol. Ant workers were removed from each container and identi ed at the species level or assigned to a morphospecies number.

Ethics approvals
All applicable international, national and institutional guidelines for the collection of ant specimens and leaf-litter material were followed. All procedures performed in studies involving animals were in accordance with the ethical standards of the institution at which the studies were conducted. The material was collected in accordance with the permits issued by SEMARNAT (license number: FAUT-0312).

Phylogenetic tree constructions
Ideally, one would use a complete, species-level phylogeny of all ant species present in your study area to calculate phylogenetic diversity, yet our current understanding of ant relationships is still limited. As an alternative, we built a genus-level phylogeny based on the tree by Moreau & Bell 70 , but using the phylogenetic relationships and divergence times within Myrmicinae from Ward 42 . This phylogeny was then pruned to keep only a single species per genus to generate a genus-level phylogeny. To maximize taxonomic coverage, we replaced genera that were missing from those studies by closely-related lineages that were not present in our dataset using other phylogenetic studies [71][72][73] . We then used the list of species (Supplementary Table S1) in our dataset to simulate a species-level phylogeny in which the relationships within genera were obtained from a Yule (pure-birth) process using the genus.to.species.tree function in the "phytools" package 74 . A total of 1000 simulated trees were obtained to account for phylogenetic uncertainty [see 75 and 76 for similar approach]. Additionally, we constructed a maximum clade credibility tree (hereafter MCC tree) which was used to summarize the uncertainty of the 1000 simulated trees. The MCC tree was constructed from the sample of the 1000 trees with the maxCladeCred function incorporated in the "ape" package 77 . Both the 1000 hypothetical trees and the MCC tree were used in downstream analyses (Supplementary Fig. S1).

Alpha and beta phylogenetic diversity
Phylogenetic alpha diversity patterns of leaf-litter ant assemblages at each site was estimated using three metrics: Faith's PD 78 , mean pairwise distance (MPD) and mean neighbor taxon distance (MNTD).
PD measures the total phylogenetic branch length that joins the basal node to the tips of all the species in the sample 78 , and is used as a measure of absolute phylogenetic diversity 79 . MPD and MNTD are two phylogenetic divergence metrics 7,80 , with the former being an estimate of the average phylogenetic relatedness between all possible pairs of taxa in a local community, whereas the latter is an estimate of the mean phylogenetic relatedness between each taxon in a local community and its nearest relatives.
To investigate which processes may be in uencing the patterns of phylogenetic alpha diversity, we calculated the standardized effect sizes (SES) for each α-diversity metric (i.e., SES.PD, SES.MPD, SES.MNTD). SES values were calculated by taking the difference between the observed value of PD, MPD and MNTD and their corresponding mean random values derived from null communities, then dividing these differences by the standard deviation across randomizations 81 . The null communities were generated by randomizing the community data matrix using the "independentswap" algorithm with 1000 iterations. Positive and negative SES.PD values indicate species in a community accumulated more or less evolutionary history than expected by null communities, respectively 82 . Whereas for MPD and MNTD, positive SES values suggest phylogenetic overdispersion, whereas clustering is inferred by negative values 7 . Statistical signi cance is inferred if SES values are greater than 1.96 or less than -1.96. All these analyses were conducted using the ses.pd, ses.mpd and ses.mntd functions incorporated in the "picante" package 83 of the R software 84 .
We investigated the patterns of β-diversity through three complementary approaches: (i) the multiple-site approach (PBD multi ), used to summarize in one value the overall dissimilarity in the mountain, (ii) the adjacent approach (PBDadj), used to investigate the unidirectional β-diversity variation focusing only on adjacent sites towards the summit, and (iii) the pairwise approach (PBDpair), used to investigate how βdiversity patterns were related to the environmental and geographical distance between all pairwise sites.
For all approaches, we calculated the total dissimilarity through the PhyloSor distance (multiple-site: PBD multi.sor , adjacent: PBDadj.sor; pairwise: PBDpair.sor) and further decomposed it into the turnover (multiple-site: PBD multi.sim, adjacent: PBDadj.sim; pairwise: PBDpair.sim) and the nestedness (multiplesite: PBD multi.nes , adjacent: PBDadj.nes; pairwise: PBDpair.nes) components. Under the phylogenetic framework, total dissimilarity captures the proportion of shared and exclusive branch lengths among assemblages, turnover measures 'true' lineage turnover and nestedness considers the differences in Faith's PD between assemblages 21 . Multiple-site calculations were obtained using the multi.phylo function, whereas pairwise and adjacent values with the phylo.beta.pair function incorporated in the "betapart" package 85 of the R software 84 .
To assess the relative contribution of the spatial turnover component to the total dissimilarity between adjacent sites, we calculated the ratio of turnover over total dissimilarity (hereafter β ratio ) following Dobrovolski et al. 86 : PBDadj.sim/ PBDadj.sor. Thus, β ratio > 0.5 indicates that total dissimilarity is determined dominantly by the turnover, and β ratio < 0.5 indicates nestedness is the dominant component 86,87 . We did not conduct such analyses for the pairwise approach since the raw-unconverted data is necessary for GDM analyses (see statistical analyses section).

Climatic predictors
To evaluate whether local climate may account for the observed phylogenetic diversity patterns, we extracted the 19 climatic variables from the bioclimatic raster available for Mexico at 3 arc-second resolution (~90 m; 88 ) coincident with the coordinates of our forty sampling points (Supplementary Fig.  S1). For practical purposes, each site was characterized by averaging these forty values. We rst divided those 19 variables into temperature-and precipitation-related subsets. Then, we used separate principal components analyses (PCA) to generate a synthetic uncorrelated climatic variable that represents the original variables contained in each climatic subset. Before PCA analyses all variables were standardized to remove the unit and were centered (mean=0, SD=1). Since the rst principal component accounted for a high variation contained in each subset of temperature (PC1 Temperature :85.3%) and precipitation (PC1 Precipitation :67.5%), we conducted the consecutive analyses using only these vectors.
The examination of variable loading in each principal component revealed that almost all variables included in the analysis (75% of the total) highly contributed (i.e., large weights) to each rst component (Supplementary Table S2). Therefore, any interpretation using PC1 Temperature and PC1 Precipitation should largely re ect the broad variation in terms of temperature and precipitation occurring along the Cofre de Perote mountain.

Statistical analyses
To evaluate which climatic variables (i.e., temperature and precipitation) better explained the phylogenetic alpha diversity, we implemented simple linear models where the 1000 SES values of each alpha metric (SES.PD, SES.MPD, SES.MNTD) were regressed individually with the rst principal component of each climatic subset (PC1 Temperature and PC1 Precipitation ). Simultaneously, we constructed a multiple regression model with the same response variables but now modeled against the additive effect of both climatic factors (PC1 Temperature + PC1 Precipitation ). Normality assumption was checked in the residuals of all the adjusted models using the Shapiro test at α=0.05. To avoid spurious interpretation, a second run of regression models was conducted including only models which met the normality assumption (models syntaxis and number of trees included in nal analyses are condensed in Supplementary Table S3).
We selected the best simple or multiple regression model explaining the phylogenetic α-diversity patterns.
For this purpose, all linear and multiple regression models were evaluated and the model with the lowest Bayesian Information Criterion (BIC) was chosen as the best model 89 . We selected BIC over the Akaike information criterion (AIC) since BIC is based on the assumption that a true model exists among the set of candidate models 90 . We considered this scenario true since temperature and precipitation (and their interaction) have been documented as the most important predictors of ant diversity (e.g., 30,33,61 ). We considered a model equally probable to the best t model if the difference in BIC (ΔBIC) between the focal model and the model with the lowest BIC were < 2. Further, we extracted the coe cients of determination (R 2 ) and the slope coe cient (β) to evaluate the proportion of variance explained by each model and the relationship between each phylogenetic alpha metric with the climatic predictors respectively. Regression models and normality tests were conducted through the lm and shapiro.test functions respectively, whereas model performance was conducted using the bictab function incorporated in the "AICmodavg" package 91 . All functions are incorporated in the R project software 84 .
To assess whether environmental ltering (climate distances) or dispersal limitation (geographical distances) better explained pairwise PBD patterns, we used Generalized Dissimilarity Modelling (GDM; 92 ). GDM uses a nonlinear matrix regression technique for analyzing spatial patterns in compositional dissimilarity, providing tted I-splines to describe the relationships between a dissimilarity matrix (response) and both climatic and geographical predictors, coupled with the partial deviance explained by each predictor 93 . Moreover, GDM standardize variables so they can be directly compared with one another and is highly robust to multicollinearity among predictors 92 . To conduct GDM, we rst converted the observed pairwise dissimilarity matrices (PBD pair.sor , PBD pair.sim , PBD pair.nes ) into a GDM site-pair table using the formatsitepair function setting the type 3 in the "bioFormat" argument. The gdm function was used to t the model which included the climatic variables (PC1 Temperature and PC1 Precipitation ) and the geographical coordinates corresponding to the centroid of the total sampling points located at each site. Finally, the function gdm.varImp was used to extract the total deviance explained by each model, the signi cance of the full model and the importance of each predictor.
Predictor importance is quanti ed as the percent change in deviance explained by the full model and the deviance explained by a model t with that variable permuted 94 . We used 1000 permutations to estimate predictor importance and full model signi cance. Since this complete procedure was ran across the 1000 matrices of each PDB component, we thus calculated the ratio between the number of signi cant values (p<0.05) out of the 1000 phylogenetic trees. GDM analyses were conducted using the functions incorporated in the "gdm" package 94 of the R-project software 84 .

Data availability
The taxonomic matrix along with the 1000 simulated trees and the Maximum Clade Credibility tree constructed in each elevational site are available on the Zenodo digital repository (doi: 10.5281/zenodo.5646220).