Diverging facets of grassland ant diversity along a Mediterranean elevational gradient

1. We address associations of taxonomic diversity (TD), functional diversity (FD), and phylogenetic diversity (PD) of ant assemblages with gradients of elevation to assess whether energetic limitations or deterministic or stochastic niche‐building processes are more relevant to the assembly of communities.


Introduction
Understanding how communities are assembled and structured in space is one of the primary goals in community ecology and conservation. From a theoretical perspective, assembly rules inform about the processes governing species interactions, community composition, and ecosystem functioning (Hooper et al., 2005). From a more application point of view, characterising diversity patterns helps understanding Correspondence: Javier Seoane, Universidad Autónoma de Madrid. Darwin, 2. Madrid 28049, Spain. Email: javier.seoane@uam.es potential threats and suggesting courses of action under global change (Walther, 2010). However, diversity patterns can differ across scales of analysis, so that it is necessary to analyse patterns at different scales to have a full understanding about the drivers of diversity (Whittaker et al., 2001). Until the last decade, studies of diversity patterns generally considered taxonomic diversity (TD), with a strong focus on species richness (Pavoine & Bonsall, 2011). However, the functional and phylogenetic facets (related to species' traits and phylogeny, respectively) of diversity are increasingly studied because they can reveal complementary information about the organisation of community assemblages and the processes structuring local communities (Götzenberger et al., 2012).
Spatial characterisation techniques generally include estimating gamma diversity, which reflects the total diversity of the species pool in a whole region. Gamma diversity can be then partitioned into its alpha and beta components, which allows examining and comparing biodiversity within and among regions (Mori et al., 2018). The alpha component refers to local diversity-a snapshot taken in a particular study unit-and is led primarily (although by no means uniquely) by immediate factors and conditions (Pavoine & Bonsall, 2011). On the other hand, the beta component shows the variation among sites, and it is best explained through landscape features and biogeographical context (Chase & Myers, 2011). Some difficulties to understand these components arise because patterns of species diversity at the local level and between assemblages can differ between habitats (Jankowski et al., 2009), while the size of study units may affect the estimations of alpha and beta diversity of assemblages (Magurran, 2004).
Functional diversity (FD) reflects differences in morphological, physiological, and ecological traits among species in a community . This facet of diversity is particularly interesting because it informs about the use of resources and the compatibility of species from the regional species pool to local habitat characteristics (environmental filters) (Martello et al., 2018). Although most FD analyses consider a single average value for each species within an assemblage, it is increasingly recognised that intraspecific variability may affect the performance and ecological success of populations and species (Forsman & Wennersten, 2016), alter the detection of patterns of ecological variation , and impact the functioning of the assemblages (Szefer et al., 2017). On the other hand, phylogenetic diversity (PD) expresses the accumulated evolutionary history, and its spatial variation suggests differences in the dispersal and sensitivity to environmental conditions among taxa (Webb et al., 2002). The study of the relationships among the three facets of biodiversity (taxonomic, functional, and phylogenetic) might hint at whether the dominant processes in community assembly are deterministic through environmental filtering and competition, or stochastic by dispersal and ecological drift (Pavoine & Bonsall, 2011). Although it is commonly expected that these facets of biodiversity are correlated, particularly if traits are evolutionarily conserved, their relationship is often complex and depends on different ecological, evolutionary, and methodological issues (Cadotte et al., 2019). Such non-obvious correlations have been shown both on broad geographic scales (Arnan et al., 2017) and in elevational gradients at smaller scales (Liu et al., 2018).
Elevation is the dominant environmental gradient in mountain habitats, and studying its role on species diversity can shed light on the possible effects of climate change on assemblages . Moreover, elevational gradients provide suitable grounds for observational and experimental studies to understand spatial variations in the diversity and functioning of ecosystems (Sanders & Rahbek, 2012). Mountains are directly related to decreasing temperature and various other abiotic and biotic environmental factors, which lead to changing ecological processes at small geographical scales (Brehm et al., 2003). Specifically, Mediterranean mountains present a progressively cooler and wetter environment at higher elevations, while valleys suffer more severe summer drought conditions (Schöb et al., 2013). Thus, the elevational gradient in Mediterranean mountains offers a suitable scenario to explore the spatial variation of the facets of diversity and assess the roles of temperature and primary productivity, which have been reported to explain species richness in arid and semi-arid ecosystems for a variety of taxa (Flores et al., 2018;Prieto-Torres et al., 2019). Among these, insects, and specifically ants, are considered good candidates to look into the mechanisms that shape diversity patterns, because they may informatively track changing ecological conditions at small spatio-temporal scales (Martello et al., 2018). Ants are globally ubiquitous, diverse, and abundant, and are credited to show strong interspecific interactions that contribute in important ways to the functioning of ecosystems (Arnan et al., 2017).
In this study, we analysed associations of elevation with TD, FD, and PD of ant assemblages in a Mediterranean mountain region. First, we examined how alpha and beta diversities vary with elevation, direct factors (primary productivity and soil temperature), and geographic distance among sampling plots (to account for distance-decay patterns). Second, we explored the relationships among the different facets of diversity to infer what are the dominant ecological processes structuring assemblages. We expected that responses of TD, FD, and PD were generally correlated if these facets of diversity were congruent in our study system. More specifically, we hypothesised (a) that energetic constraints based on temperature and productivity limit the total amount of species, the set of taxonomic groups, and the amount of viable functional space, so that we predict aligned responses of alpha TD, FD, and PD to such constraints (Seoane et al., 2017); (b) that energetic constraints determine the way species change their niches through their evolutionary history (Cadotte et al., 2008), so that we predict aligned responses of functional and phylogenetic dissimilarities between species within communities to energetic constraints; (c) that niche differentiation is determinant for species to coexist in assemblages (Liu et al., 2018), so that we predict (alpha) functional and phylogenetic overdispersion (higher variance of traits and phylogeny than random expectations); and (d) that energetic limitations are more relevant than stochastic processes (dispersal and drift) (Götzenberger et al., 2012), so that we predict larger effects of environmental distances over geographical distances on dissimilarities among communities (beta diversity). Finally, we hypothesised that if the above ecological mechanisms are not at work, no coherent patterns should be detected. Although we are mostly interested here in our study system to address those general hypotheses, we also underline that this is, to the best of our knowledge, the first assessment of the different components and facets of diversity in ant communities and its relationship to elevational gradients.

Study area and sampling design
We studied the spatial variation of ant diversity along an elevational gradient in the Guadarrama Mountains of Central Spain. This area has a typical continental Mediterranean climate, with severe summer droughts ameliorated with elevation. Mean annual temperatures range from 15 ∘ C in the foothills (500 m.a.s.l.) to 4 ∘ C at summits (peaking at 2428 m.a.s.l.), and mean annual rainfall ranges from 550 to 1500 mm (Ninyerola et al., 2005). Forests of holm oak (Quercus ilex subsp. rotundifolia), Pyrenean oak (Quercus pyrenaica), and Scots pine (Pinus sylvestris) change with elevation and intersperse with natural and semi-natural grasslands, which are abundant due to extensive livestock grazing and the activity of wild herbivores.
The sampling was conducted in summer, the period of higher ant activity in the mountains of central Spain. In July 2016, we sampled ant assemblage composition in 18 localities, all of them south-facing grassland sites, along the range (six elevation belts replicated three times, approximately at 650, 1000, 1300, 1700, 2000, and 2250 m) that covered an elevational gradient of 1712 m (range 640-2352 m; Fig. 1, Appendix S1). At each locality, we established three 20 × 15 m plots separated by 50 m, avoiding heterogeneous zones, strong slopes, or scrub patches (54 plots in total). We placed sets of 12 pitfall traps per plot in a grid of 3 × 4 traps distanced by 5 m for 7 consecutive days (648 traps in total). Traps were small plastic tubes (2.5 cm diameter × 5 cm depth) buried at ground level and contained 70% ethanol and 30% ethylene glycol solution as preservative. We then sorted ant workers (including all non-reproductive castes) in the laboratory and identified them to the species level with a binocular microscope. Reproductive individuals were excluded, since they could be in the nuptial flight or dispersal phases, and therefore we cannot infer a successful colonisation of the species in the local habitat from their occurrence. We followed the nomenclature considered as accepted at antcat.org and used the general keys available at ant.org, complemented, when necessary, with specific keys for certain genera. Pitfall trapping is an effective method for sampling grassland ant communities, and it has not been reported to cause significant damage to the ant colonies. The small diameter of the traps minimised the capture of larger sized invertebrates. One species in the area is listed as vulnerable by The International Union for Conservation of Nature (Formica dusmeti); however, it occurs primarily in forest habitats, so, as expected, very few workers were captured.

Direct factors
We estimated both soil temperature and productivity in our study plots. We buried data loggers (HOBO Pendant, 0.4 ∘ C precision) at ground level to obtain mean soil temperatures. The loggers recorded soil temperatures every half an hour during a complete year (March 2016 to March 2017). As a proxy for primary productivity in plots, we measured the Normalised Difference Vegetation Index (NDVI). We used bottom of atmosphere reflectance values for bands four (red) and eight (near infrared) of the MultiSpectral Instrument sensor onboard Sentinel 2A satellite at 10 m spatial resolution (provided by the Spanish National Geographic Institute; http://www.ign.es/iberpix2/ visor/). Raw NDVI values were provided per each 10-day period and averaged for a whole year. By considering yearly measures, Fig. 1. Map of survey sites where grassland ant communities were sampled (n = 18 sites, with three sampling plots in each). Shapes represent each of the three spatial replicates, while colours show the six elevational belts we considered along the Guadarrama range (from light grey at around 650 m to dark grey at >2200 m). The inset shows the location of the study area (rectangle in red) within the Iberian Peninsula.
we assumed that abiotic limiting factors act throughout the year, since ant colonies remain in the field and this group of insects cannot use avoidance strategies such as migrations or the adoption of resistance states.

Ant functional traits
We selected a set of morphological traits most relevant to feeding and foraging strategies determining species' ecological roles (Arnan et al., 2014). Thus, we measured head width, eye width, scape length, and tibia length of hind leg for three randomly selected individuals of each species present in each sampling plot. Six rare species, accounting for a marginal proportion in the assemblages, had fewer individuals. However, this shortage of information has been reported not to affect significantly previous examinations of FD of ants (Martello et al., 2018). Due to the substantial differences in size among species, traits were log-transformed before any calculations. All the considered traits are inherently related to body size, and hence strongly correlated between them (Pearson's correlation between pairs of traits was >0.8 in all cases). Because of this, for eye width, scape length, and tibia length, we performed regressions on head width and used the residuals of these regressions as response variables in the following models instead of the original traits. Thus, if, for example, an individual has a large residual value in the model for eye width, then the eyes are wider than expected for its head width. Following Martello et al.'s (2018) study, we estimated the mean trait values for each species and used these average values to perform a Principal Components Analysis (PCA). We then used trait values measured on individuals to predict their position in this PCA. This strategy gives equal weight to all species for the estimation of the PCA, regardless of their prevalence in the sampled region, but allows us to consider measures at the individual level. We used Horn's parallel analysis to determine the number of principal components retained (R package paran: Dinno, 2018).
We used the scores of the individuals in the retained orthogonal axes to estimate the trait probability density (TPD; Carmona et al., 2016) functions of each species using the TPD package (Carmona, 2019). The TPD function of species (TPDs) reflects the abundance of the different trait combinations among their individuals. TPDs functions can be combined to estimate the TPD function of a community (TPDc), which reflects the relative abundance of traits in an assemblage and can be used to estimate several aspects of FD . We estimated the TPDc of each plot considering the species that were present in the pitfall traps.

Ant phylogeny
We built a complete phylogeny for the 43 ant species of our study (Appendix S2). This was based on the species-level tree in Arnan et al. (2017), which in turn started from the super tree of Moreau and Bell (2013). To get our tree, we manually deleted unrecorded species and added species (n = 11) to the basal tree with Mesquite 3.0 by combining the species-level phylogeny with our own expert opinion based on the taxonomy and morphology of the study species. We then built an ultrametric tree with Grafen's Rho transformation because reliable estimates of branch length are currently unavailable.

Alpha and beta taxonomic, functional, and phylogenetic ant diversities
We estimated alpha diversity in plots for the three facets of diversity. For TD, we considered species richness (the number of species recorded) in each plot. For the functional and phylogenetic facets, we estimated two indicators of alpha diversity: first, functional richness (the total amount of functional space occupied by the species present in each plot; Carmona et al., 2016) and Faith's PD (the total evolutionary history contained in each plot ;Faith, 1992). These indices can be considered as the functional and phylogenetic analogues of species richness. Indicators of richness can convey information about the total diversity of functional and evolutionary strategies that can be found in a community, and they can be used to detect environmental filtering (e.g. Carmona et al., 2012). In addition, we estimated the mean pairwise dissimilarity (MPD) between species (using the 'melodic' R function: de Bello et al., 2016). MPD represents the expected dissimilarity (functional or phylogenetic) between two different species chosen at random from the assemblage. The functional dissimilarity between pairs of species was estimated as 1 − overlap in their respective TPDs functions. Phylogenetic dissimilarities between pairs of species were estimated as the cophenetic distance in the phylogenetic tree, which was standardised to have a maximum value of 1 by dividing all dissimilarities by the maximum dissimilarity between species in the dataset. Indicators of the average dissimilarity between species are helpful to detect changes in niche differentiation patterns along the gradient .
In addition, following Hevia et al. (2016), we performed a partition of diversity at different hierarchical levels: within plots, within localities, within elevational belts, and between elevational belts. For this, we started estimating the different indicators of richness (taxonomic, functional, and phylogenetic) at the plot level, as described earlier. Then, we sequentially pooled the plots from the next hierarchical level (first the three plots in a locality, then the three localities in an altitudinal belt, and finally the six altitudinal belts) to estimate the indicators of richness at that level. This analysis can only be performed with richness indicators, since they increase monotonically with sample size (e.g. species richness cannot decrease when a new plot is added), which is not true for MPD.
We then estimated beta diversities between ant assemblages. Taxonomic beta diversity was calculated as Sørensen dissimilarity between pairs of plots ('beta.pair' in package betapart: Baselga et al., 2018). Functional beta diversity was estimated as the dissimilarity between the TPDc functions of each pair of plots (TPD package: Carmona, 2019). Finally, phylogenetic beta diversity was estimated with the PhyloSor index based on the Sørensen dissimilarity index (we used the package betapart, following Bryant et al., 2008).
FD and PD indexes could be related to species richness (e.g. functional richness and Faith's PD can only increase as new species are added). To correct for this effect, we compared the observed diversity values with values obtained using a null model generated by a matrix swap of the species × plot presence-absence matrix (Oksanen et al., 2019). We performed 1000 randomisations using the quasiswap algorithm ('permatswap' function from the package vegan: Oksanen et al., 2019). Then, we estimated the standardised effect size for each diversity index (SES index = [observed index − mean of simulated indexes]/SD of simulated index).

Statistical analysis
We analysed the relationships among alpha diversity components, first with elevation (to assess elevational patterns) and then with temperature and productivity (to understand those patterns) using generalised additive mixed models that included locality as a random effect (GAMMs; 'gamm4' function from gamm4 package). We also used generalised additive models (GAMs) to examine the relationships of taxonomic, functional, and phylogenetic beta diversities with the differences among assemblages in elevation, productivity, and temperature, and the geographic distances among them (with 'mgcv' function from 'mgcv' package). GAMs and GAMMs model responses as smooth functions of predictors. They may fit linear responses but are particularly useful to describe nonlinear relationships without the need to impose in advance a particular shape (as, for instance, a quadratic polynomial). Unanticipated or little-known non-linear responses can then be accommodated. Poisson errors were used for ant species richness (TD), and Gaussian errors were used for the rest of the response variables (FD and PD related). Significance of terms was assessed by nested comparisons with log-likelihood ratio tests. In addition, for the beta diversity models, we performed partial GAMs to estimate the proportion of each facet of diversity explained by pure and joint effects of the predictors. To correct for using the same plot several times in these models (one in each pairwise comparison to calculate beta diversities), we calculated P values by resampling the response variable (n = 1000 resamples) and comparison with null models (Gotelli & Graves, 1996). The significance of each predictor was estimated comparing the explained deviance in the original model with the distribution of explained deviances resulting from fitting the model to randomly permuted matrices of beta diversity indexes (N = 999 times). Pure and joint effects were estimated according to a variation partitioning scheme based on simple algebra with the R 2 of the models (for a similar approach, see Seoane & Carrascal, 2008).
Finally, to characterise the degree of niche conservatism in ant species' functional traits along the elevational gradient, we tested the presence of a phylogenetic signal in each trait. We used Pagel's lambda ( ) test (Pagel, 1999), which is one of the most widely used in ecology (with 'phylosig' function in phytools package in R). This index assumes a Brownian motion evolutionary process. Its values range from 0 to 1: values close to 0 indicate no phylogenetic signal (the trait has evolved independently of phylogeny), and values close to 1 indicate trait evolution according to Brownian motion.
Sample completeness was assessed by building rarefaction curves for species richness using sampling-unit-based incidence data, according to methods in Chao et al. (2014) as implemented in Hsieh et al. (2016). The curves allow estimation of sample coverage, which gives the proportion of the total number of individuals in a community that belong to the species already represented in the sample (Chao & Jost, 2012). We also estimated asymptotic species richness to compare with observed species richness (Chao, 1984).

Results
We found a total of 43 ant species from 18 genera and three subfamilies (Myrmicinae, Formicinae, and Dolichoderinae) along our elevational gradient (see Appendices S2 and S3 for phylogenetic tree used and the list of species recorded in each plot, respectively). The three subfamilies, as well as several intermediate and lower level clades (e.g. genera Tetramorium, Cataglyphis, Tapinoma, Temnothorax) were distributed along most of the elevational gradient (Appendix S3). All ant species were native to the Iberian peninsula. Ant occurrence data are available at figshare with DOI 10.6084/m9.figshare.14596014.
Observed and estimated (asymptotic) species richness were highly correlated (r = 0.93) and sampling completeness was very high (>90% in 52 out of 54 plots; see Appendices S7-S10). This suggests that our sampling effort was satisfactory, and we opted for using observed species richness in the following analyses.
Horn's parallel analysis retained the first two components of the PCA (which accounted for 81% of the total variance in species mean trait values). Eye width, scape length, and tibia length (relative to head width) were positively related to the first component (which accounted for 56% of trait variance), whereas the second component (accounting for 25% of variance) was strongly related to head width (Appendix S4). Ant trait data are available at figshare with DOI 10.6084/m9.figshare.14596065.
Soil temperature during sampling periods ranged from 9.0 to 62.1 ∘ C, averaging 20.8 ∘ C in the highest locality and 36.0 ∘ C in the lowest one (Appendix S1). Elevation was linearly related to soil temperature (linear regression of mean temperature per locality on elevation R 2 = 0.94) and quadratically to productivity (R 2 = 0.58; highest productivities were observed at 1000-1400 m). Given the strong collinearity between temperature and elevation, we decided not to explore temperature further.
The relationship between species richness and elevation was quadratic (R 2 = 0.65; P < 0.001; Fig. 2a) showing a peak at ∼1200 m, while the relationship with productivity was strictly increasing (R 2 = 0.45; P = 0.002; Fig. 3a). The unexplained variation of both models could be attributed mainly to differences among localities (elevation: locality = 2.41; 72% of residual variation; productivity: locality = 2.32; 70% of residual variation), and secondarily to differences among plots ( residual = 1.52; 28%; productivity: residual = 1.52; 30%), all due to unmodeled factors. Functional richness was negatively related with elevation ( Fig. 2b, left; R 2 = 0.36, P < 0.001) and positively (albeit not reaching significance) with productivity (Fig. 3b, left; R 2 = 0.12, P = 0.085), while phylogenetic richness did not show any significant response to either factor. When functional and phylogenetic richness values were not controlled by species richness, the relationships of these metrics with elevation (decreasing) and productivity (increasing) were stronger and significant (Appendix S5), which suggests that phylogenetic, but not all of the functional, variation along the gradients was driven by species richness in the plots. Standardised functional and phylogenetic alpha MPDs had positive linear relationships with elevation (Fig. 2b,c, right) and negative relationships with productivity (Fig. 3b,c, right). However, the slopes of these relationships were not significant in all cases. The unexplained variation in functional MPD could be attributed largely to differences (other than those due to elevation) among plots (elevation and productivity models: locality = 0.03 and ε = 0.70, the last one accounting for 95% of residual variation). The unexplained variation in phylogenetic MPD was more evenly attributed to the differences among localities and plots (elevation: locality = 0.81; 60%; productivity: locality = 0.88; 61%). Assemblages did not exhibit patterns of clustering or overdispersion in neither functional nor phylogenetic alpha MPDs (SES values were mostly within 2 units from mean values except for a few exceptionally low diverse plots, Figs 2b,c and 3b,c).  Bryant et al., 2008). The standardised effect sizes of functional and phylogenetic richness show, respectively, the functional and phylogenetic space occupied by the assemblages, while the corresponding main pairwise dissimilarities show changes in species' niche and phylogeny differentiation.  Bryant et al., 2008). The standardised effect sizes of functional and phylogenetic richness show, respectively, the functional and phylogenetic space occupied by the assemblages, while the corresponding main pairwise dissimilarities show changes in species' niche and phylogeny differentiation.

Fig. 4.
Partitioning of species (taxonomic diversity, TD), functional diversity, and phylogenetic diversity richness within assemblages at the different scales of the sampling hierarchy (plot, locality, elevational belt, and study region). Diversity within sampling plots is the average amount of diversity that is contained at the sampling unit scale (e.g. the average number of species in plots in the case of TD); diversity within sampling localities is the average amount of diversity that is attributed to differences between the plots in a locality (e.g. how many species are added from the previous level when we aggregate the three plots in a locality); diversity within altitudinal belts is the average amount of diversity that is attributed to differences between the localities in an altitudinal level (e.g. how many species are added from the previous level when we aggregate the three localities in the same altitudinal level); finally, total diversity represents the amount of diversity contained in the whole dataset (e.g. how many new species added from the previous level when we aggregate the six altitudinal belts).
Partitioning of TD (Fig. 4) showed that most of the species richness of ant assemblages was attributed to differences between elevations (19 species), in contrast to differences between plots (7 species) or localities (10 species). However, the partitioning patterns of functional and phylogenetic richness (Fig. 4) differed from species richness. Most of the variability in functional and, to a lesser extent, phylogenetic richness occurred within the sampling plot scale, and then, between elevations.
Taxonomic beta diversity was related to differences in elevation and productivity, and with geographic distances between plots (GAM: R 2 = 0.61; P < 0.001 for all terms; Fig. 5). The dissimilarity between the compositions of ant assemblages increased with the difference in elevation between sites, for the most part of the elevational gradient. Similarly, larger differences in productivity and larger geographical distances between sites corresponded to larger differences in taxonomic composition, although for these variables the relationship with beta diversity was more quadratic in shape and peaked at mid-values of the productivity and distance gradient (Fig. 5a). Beta (standardised) FD was also related to elevation, productivity, and geographic distances, although the explanatory capacity of these predictors was weaker than for TD (GAM: R 2 = 0.27; P < 0.001 for all terms; Fig. 5b). Beta FD increased with differences in elevation with a diminishing slope along the gradient, and also with geographical distances, in this case reaching a plateau at the end of the gradient. The relationship of beta FD with between-plot differences in productivity was flatter, although it showed an intermediate peak (Fig. 5b). Beta PD was less explained by the variables used in our models (GAM: R 2 = 0.16, Fig. 5c). The relationship with elevational and geographical distances was positive with a diminishing slope along the gradient that eventually plateaued (P < 0.001 for both terms, Fig. 5c). No relationship was shown with differences in productivity (P = 0.150, Fig. 5c). Pure and joint effects of elevation were largely the most relevant determinants of beta diversities (Fig. 5, rightmost column). However, pure and joint effects of geographical distances were also relevant to explain FD and PD.
Ant functional traits, except scape length, showed significant phylogenetic signals along the elevational gradient (Table 1).

Discussion
We found some revealing correlations among the responses of the different facets (and metrics) of ant diversity to environmental and geographical gradients. The patterns of variation of TD, FD, and PD were remarkably similar when these facets were described with raw metrics that do not tease apart the effect of species richness; they all showed strong and coincident responses to the gradients. However, when the facets were described factoring out the effect of species richness, their responses to the gradients and the similarities among them weakened much (cf. Fig. 2 and Appendix S5). TD and FD changed roughly in a similar manner with elevation and  Bryant et al., 2008). We calculated P values by resampling the response variable (n = 1000 resamples) and comparison with null models (Gotelli & Graves, 1996). productivity, both within communities, as predicted by our hypothesis (a), and across communities, with altitude-that is, energetic limitations-playing a more important role than geographical distance-that is, dispersal and drift-following hypothesis (d). However, PD did not show clear responses to any of the studied factors, hence deviating from hypothesis (b) that predicted changes in PD along the gradient. Overall, these partial redundancies among facets of diversity suggest that the spatial variation of diversity in our study system is mainly related to the sheer number of species. However, we cannot discard that both environmental filtering (evidenced by decreasing functional richness with altitude) and limiting similarity (suggested by the relatively constant dissimilarity between coexisting species among altitudes, which provides support for hypothesis (c)) could be acting in concordance, reducing the number of species that are finally present in communities. In the rest of this section, we discuss these results and their potential implications. The variation of taxonomic and functional richness with elevation and productivity matches global patterns of distribution of diversity. TD for a variety of taxa has been repeatedly found to increase with sampling plot area (Nogués-Bravo et al., 2008), heterogeneity (Stein et al., 2014), and productivity (Hawkins et al., 2003). We controlled sampling area and heterogeneity by design, sampling equally sized plots with a narrow range of conditions (constraining slope and aspect) within a single, simple habitat, which rules out their contribution to the patterns we obtained. On the other hand, water and temperature interact to compound productivity, and there is some evidence that species richness is primarily related to water availability in arid ecosystems and to energy regimes in cooler environments (Whittaker et al., 2007). The underlying mechanisms are not completely known, but prominent explanations revolve about more productive sites offering more avenues for niche packing and complementarity, or allowing the growth of more abundant populations that are thus less vulnerable to local extinction (Hurlbert & Jetz, 2010, but see Storch et al., 2018. Our study tracks smaller gradients of the same energy variables in a mountain system within a Mediterranean climate, and results suggest that ant species richness at local scales may be more limited by cool temperatures at high elevations and by hot temperatures and low water availability at valleys. If productivity increases species richness by favouring niche diversification, then TD and FD should roughly correlate along the gradients, as we found. Contrastingly, PD was not related to our contemporary descriptors of available energy. Note, however, that phylogenetic and functional distances between pairs of species showed a very poor correlation (Spearman rank correlation = 0.13), with some closely related species having much larger functional distances (Appendix S6). This suggests that changes in functional richness along the elevational gradient were to a large extent due to replacement of functionally different species from the same clades (de Bello et al., 2017). The pattern goes against the view that phylogenies may effectively substitute functional descriptors to understand the assembly of communities and their effects on ecosystems, and suggests that different facets of diversity contribute with complementary information (Cadotte et al., 2008;Devictor et al., 2010). However, we acknowledge that our study is based on a moderately sized regional species pool (around 100 species in all the habitats in the region, according to Martínez, 1987) occurring on a small area (roughly 3000 km 2 ), which may be too constraining for the assemblages to show a clear spatial structure of their evolutionary history.
The average functional and phylogenetic distances between pairs of species within a community (the mean pairwise dissimilarities in Figs 2 and 3) did not change along the gradients of elevation and productivity. Under the assumption that phylogenetically close or functional similar species would likely share resources due to their physiological and ecological affinities, we could take these between-species dissimilarities as an indirect measure of competition pressure (Wiens et al., 2010). If so, we have to conclude that current interspecific competition is not modulated by the contrasting environmental gradients we studied. On the same line, we did not find overdispersion of functional and phylogenetic diversities within a community, which again casts doubt about our initial assumption that niche differences among ant species were determinant for coexistence (Nipperess & Beattie, 2004). These results are unexpected as ants have been frequently reported to be fierce competitors and beg the questions of what could be promoting coexistence among similar species (Ellwood et al., 2016). Several mechanisms that could act simultaneously have been shown to facilitate species coexistence in competition-laden tightly packed ant assemblages (Andersen, 2008). We suspect that trade-offs based on behavioural differences in the use of resources and temporal shifts of activity, such as those between discovery and dominance abilities (Bertelsmeier et al., 2015;but see: Parr & Gibb, 2012) or between thermal tolerance and competitive ability (Cerdá & Retana, 1997) allow niche partition and coexistence in our studied assemblages. An alternative view is that interspecific competition among ant species that commonly have broad diets might indeed be weak and it has been an overrated process to explain the assembly of ant communities (Cerdá et al., 2013;Gotelli & Ellison, 2002).
Interestingly, that the average phenotypic dissimilarity between coexisting species remained roughly constant along the gradients suggests a role for limiting similarity as a mechanism precluding the coexistence of species with very similar traits (Wilson, 2007). Thus, interspecific competition might be a determining factor for the structure of communities just beyond the threshold imposed by such limiting similarity (Ellwood et al., 2016;Nipperess & Beattie, 2004). Moreover, limiting similarity is consistent with the fact that most of the functional variability within elevational belts (i.e. including the localities at similar elevations in the different gradients) was found at the plot scale. This result suggests that the viable portion of the functional space at a given elevation is quickly filled at fine scale, as other studies have found (with ants: Hevia et al., 2016;and saproxylic beetles: Micó et al., 2020). Functional redundancy is therefore low and independent of elevation. This implies that the functional effect of losing one species should be similar in all localities, despite losing one species in the poorer assemblages (those at higher elevations) implies a larger biodiversity loss in relative terms. Ant species loss is regarded as a plausible scenario in Mediterranean mountains in the face of climate change (Flores et al., 2018), and it is likely to affect ecosystem functioning to a large extent due to the key role played by ants on processes and interactions such as seed predation, seed dispersal, soil formation, invertebrate predation, or a number of arthropod-ant mutualisms (Del Toro et al., 2012;Silvestre et al., 2019). Finally, limiting similarity and environmental filtering acting simultaneously could be a mechanism contributing to the loss of species along elevation that we found. This is because if the amount of functional space that is viable decreases along elevation, but species similarities remain constant, then a lower amount of species could enter in communities as elevation increases.
Then, functional richness and the average functional dissimilarity between species change differently along the gradients, although these differences are more conclusive regarding elevation than productivity (compare plots in Figs 2b and 3b). The total amount of functional space occupied by the ant assemblages diminishes with elevation but, at the same time, the differences between their constituent species do not change, which is much coherent with an environmental filtering of ants' functional responses. Increases in elevation continuously erode the functional space occupied by the assemblage, forbidding some combinations of traits, but the species that occupy the remaining functional space are as (dis)similar as species at low altitudes are. Thus, species loss with elevation is not random, and the assembly of ant communities must be conditioned by the selective pressure of direct factors that correlate spatially with elevation (see a similar result in Reymond et al., 2013).
Finally, environmental and geographical gradients explained the variation in metrics of beta diversity describing differences among assemblages, although to a different extent and, to our view, somewhat surprisingly. Elevation and productivity gradients accounted for much of the taxonomic variation among assemblages but only explained a moderate to low variation of functional traits and phylogenetic relationships, while geographical distances were more explanatory of PD. We assumed that the gradients in elevation and productivity describe energetic limitations (temperature, water availability, feeding resources) that may act as selective agents, while geographical distances would be more related to stochastic processes of community assembly (dispersal and drift). If so, our results would suggest a main role of spatially variable selection on species composition (Vellend, 2016) so that different species are at advantage under the different environmental conditions imposed indirectly by elevation (the species sorting paradigm of (Leibold et al., 2004). However, such selection is not clearly shown in a similar sorting of functional responses and phylogenetic lineages. Instead, geographical distances appear more relevant to explain functional and phylogenetic similarities among species, which suggests a more important role of stochastic processes determining these facets of diversity. This major role of distance and the different spatial patterns of the facets or diversity are intriguing, particularly when considering that we surveyed a moderately sized area within a single biogeographical region (see Liu et al., 2018, who found weaker effects of geographical distance on a wider study area). Besides, species composition changes more rapidly along the environmental gradients than species' traits and relatedness (note that the slopes of the diversity-elevation and diversity-productivity relationships are higher for taxonomic beta diversity). This faster species' turnover implies some redundancy: along the gradients, species are substituted by others similar in functional traits and close in phylogeny.
In conclusion, our study recalls that diversity is a complex concept, and that elevation is a strong determinant of the diversity patterns of ant assemblages. We have shown that, in the studied Mediterranean gradient, elevation has a strong effect on TD and FD, but not on PD, possibly due to functional divergence being faster than phylogenetic differentiation between evolutionarily close species. This contributes to our knowledge of the structuring of assemblages and to design more effective conservation strategies. For example, the fact that functional and phylogenetic dissimilarity between species did not change across the gradient implies that levels of functional redundancy are relatively low even in assemblages with many species. Since all species bear some unique phenotype or evolutionary history, the functioning of these systems is probably rather sensitive to species losses that could impact on several ant-mediated processes and ecosystem services (Arnan et al., 2017;Azcárate et al., 2005). However, localities where ant assemblages rely on a few species, such as high mountain assemblages, seem to be particularly vulnerable to species losses due to the warning projections for climate change and land use modifications in the Mediterranean mountains (Nogués-Bravo et al., 2008). work within the Madrid's Government research group network REMEDINAL3-CM (S-2013/MAE-2719). Carlos P. Carmona was supported by the Estonian Research Council (PSG293) and the European Union through the European Regional Development Fund (Centre of Excellence EcolChange). Mariola Silvestre was supported by an FPI grant from MINECO. Special thanks to Héctor Miranda, Laura Morgado, and Miguel Blázquez for laboratory assistance and Violeta Hevia, Joaquin Calatayud, and Anders Forsman for your tips. Alberto González helped us to make the map and Xavier Espadaler to identify some species. We appreciate the use of data from Instituto Geográfico Nacional and Nuria Plaza for provided access with the satellite imagery.
The authors declare that there is no conflict of interest.

Author contributions
JS, FMA, and MS conceived the ideas and designed the experiments. MS, JS, and FMA performed the experiments. CC, MS and JS analysed the data. MS and JS led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.