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: 30-50 m, 610-670 m, 900-1010 m, 1470-1650 m, 2020-2230 m, 2470-2600 m, 3070-3160 m and 3480-3540 m, however, for simplicity, we will refer to each site as discrete units (i.e. 0, 600, 1000, 1500, 2100, 2500, 3100, 3500 m).
We obtained 320 m2 leaf-litter samples in 8 study sites (see 32 for a complete description of sampling design). In each 1-m2 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 fine litter. The concentrated fine 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 identified 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 & Bell70, but using the phylogenetic relationships and divergence times within Myrmicinae from Ward42. 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 studies71–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” package74. 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” package77. 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 PD78, 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 sample78, and is used as a measure of absolute phylogenetic diversity79. MPD and MNTD are two phylogenetic divergence metrics7,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 influencing 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 randomizations81. 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, respectively82. Whereas for MPD and MNTD, positive SES values suggest phylogenetic overdispersion, whereas clustering is inferred by negative values7. Statistical significance 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” package83 of the R software84.
We investigated the patterns of β-diversity through three complementary approaches: (i) the multiple-site approach (PBDmulti), 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: PBDmulti.sor, adjacent: PBDadj.sor; pairwise: PBDpair.sor) and further decomposed it into the turnover (multiple-site: PBDmulti.sim, adjacent: PBDadj.sim; pairwise: PBDpair.sim) and the nestedness (multiple-site: PBDmulti.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 assemblages21. 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” package85 of the R software84.
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 component86,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 first 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 first principal component accounted for a high variation contained in each subset of temperature (PC1Temperature:85.3%) and precipitation (PC1Precipitation: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 first component (Supplementary Table S2). Therefore, any interpretation using PC1Temperature and PC1Precipitation should largely reflect 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 first principal component of each climatic subset (PC1Temperature and PC1Precipitation). Simultaneously, we constructed a multiple regression model with the same response variables but now modeled against the additive effect of both climatic factors (PC1Temperature + PC1Precipitation). 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 final 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 model89. 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 models90. 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 fit model if the difference in BIC (ΔBIC) between the focal model and the model with the lowest BIC were < 2. Further, we extracted the coefficients of determination (R2) and the slope coefficient (β) 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” package91. All functions are incorporated in the R project software84.
To assess whether environmental filtering (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 fitted 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 predictors92. To conduct GDM, we first converted the observed pairwise dissimilarity matrices (PBDpair.sor, PBDpair.sim, PBDpair.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 fit the model which included the climatic variables (PC1Temperature and PC1Precipitation) 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 significance of the full model and the importance of each predictor. Predictor importance is quantified as the percent change in deviance explained by the full model and the deviance explained by a model fit with that variable permuted94. We used 1000 permutations to estimate predictor importance and full model significance. Since this complete procedure was ran across the 1000 matrices of each PDB component, we thus calculated the ratio between the number of significant values (p<0.05) out of the 1000 phylogenetic trees. GDM analyses were conducted using the functions incorporated in the “gdm” package94 of the R-project software84.