Spatial Scale Dependence of Ecological Factors That Regulate Functional and Phylogenetic Assembly in a Mediterranean High Mountain Grassland

Introduction

Coexistence theory considers that community assembly is a process regulated by multiple mechanisms, both stochastic and deterministic, which act at different spatial scales and that are subject to different ecological factors (Palmer, 1994; Chesson, 2000). Understanding how these assembly mechanisms and their drivers operate simultaneously at different spatial scales and along environmental gradients is critical for forecasting community responses to environmental changes, including the major issue of ongoing global warming. Approaches based on functional traits, which are valuable indicators of resource use strategies and fitness, can provide critical information about how ecological factors affect assembly because only species with the appropriate set of traits can be part of the local (or realized) plant community (Cornelissen et al., 2003; Violle et al., 2007). In addition, phylogeny seems to be involved in assembly processes because it summarizes the evolutionary history of the local species pool (Losos, 1996; Webb et al., 2002; Cavender-Bares et al., 2009; Kraft and Ackerly, 2010; Mouquet et al., 2012; Gerhold et al., 2015), and it is often used as a proxy of current functional information as it is assumed that evolutionary diversification has generated predictable trait diversification (Flynn et al., 2011). Thus, comparisons of the functional (and/or phylogenetic) diversity of species assemblages and those expected from null assemblages have been performed routinely during the last two decades to identify the assembly processes that structure plant communities at different spatial scales (Díaz et al., 1998; Cavender-Bares et al., 2006; Kraft and Ackerly, 2010; Chalmandrier et al., 2017; López-Angulo et al., 2021).

Among the so-called deterministic processes, environmental filtering and limiting similarity are two major ecological forces that drive community assembly. It is known that these mechanisms act at different spatial scales on the regional species pool (Levin, 1992; Swenson et al., 2006; Emerson and Gillespie, 2008; Götzenberger et al., 2012) but their importance in specific plant communities has traditionally been studied at a single spatial scale under similar environmental conditions (but see Kraft and Ackerly, 2010; Messier et al., 2010; de Bello et al., 2013b; Carboni et al., 2014; Chalmandrier et al., 2017) and/or in systems with high species turnover. Environmental filtering, which is considered to act at relatively large scales, is expected to lead to a community where species are more functionally similar or closely related (i.e., functional and phylogenetic clustering; Weiher and Keddy, 1995; Leibold, 1998; Webb, 2000) because filters will select species that can thrive under the local conditions (but see Gerhold et al., 2015). By contrast, limiting similarity (sensu MacArthur and Levins, 1967), which is related to the existence of niche differences where plant-to-plant competitive interactions prevail (i.e., smaller scales), predicts that competitiveness among species with similar requirements leads to niche segregation where the coexisting species are functionally and/or phylogenetically different (Mason and Wilson, 2006; Kraft et al., 2008), thereby implying trait and phylogenetic overdispersion for exploiting the diversity of niches at small scales. However, alternative combinations have also been described. For example, strong environmental heterogeneity may potentially lead to a functional overdispersion pattern at relatively small spatial scales because different niches are available (Mayfield and Levine, 2010; de Bello et al., 2013b), which would determine intense turnover. In addition, trait clustering may be linked to biotic competitive processes at small spatial scales when the dominance of closely related species excludes weaker competitors (Chesson, 2000; Mayfield and Levine, 2010).

Evaluating patterns of site-to-site functional variation, i.e., a type of β-diversity referred to as functional heterogeneity in the following, among spatial scales may be necessary for understanding the importance of mechanisms involved at both local and regional scales in community assembly (Seabloom et al., 2005; Siefert et al., 2013; Yang et al., 2015). Thus, high functional heterogeneity would be expected when comparing the functional variation between sites or assemblages (at large spatial scales) simply because geographic and environmental distances often imply strong species turnover (i.e., different community types) or shifts in their functional traits (Myers et al., 2013). However, functional heterogeneity could be stabilized at the regional scale when a strong environmental filter acts and selects species with a particular combination of associated functional traits (Körner, 2003). In addition, functional heterogeneity could be increased further at local scales (within sites) by small-scale environmental variations such as those induced by soil heterogeneity (Freestone and Inouye, 2006) or species interactions (Pescador et al., 2014).

Mountains exhibit steeper environmental gradients over short geographic distances and at various scales (Körner, 2003). Thus, mountains are excellent systems for investigating how assembly processes operate at different spatial scales and the relevance of ecological factors for structuring plant communities (McCain and Colwell, 2011). Previous trait- and phylogenetic-based studies in high mountain communities suggest that environmental filtering increases along the elevation by affecting functional and phylogenetic patterns (de Bello et al., 2013a; Pistón et al., 2016; Xu et al., 2017; Le Bagousse-Pinguet et al., 2018). In addition, these patterns can be shaped by local factors in alpine communities, such as the potential solar radiation, soil pH, or nutrient heterogeneity (Cadotte et al., 2011; Bernard-Verdier et al., 2012; Luo et al., 2016). However, functional and phylogenetic patterns have generally been assessed across large environmental gradients, which imply high species turnover with major changes in the habitat type, and thus it is difficult to understand assembly mechanisms at different scales.

In the present study, we investigated functional and phylogenetic patterns as well as their ecological predictors by considering a nested spatial design in a single plant community type located along a 500-m elevation gradient in a Mediterranean high mountain area. This community is characterized by strong environmental filtering acting along the whole study gradient, which excludes species that are vulnerable to cold (Pescador et al., 2016), and a severe summer drought, especially at lower elevations (Giménez-Benavides et al., 2007; Olano et al., 2013). These conditions limit the species pool to specialists that are well adapted to such conditions, thereby resulting in the same type of plant community with low species turnover (Pescador et al., 2015). To the best of our knowledge, a low turnover scenario has not been considered in previous studies for assessing functional and phylogenetic patterns, but it probably represents the best system for investigating how a habitat-specific species pool is affected by ongoing global warming before a change in turnover. Thus, we hypothesized that an environmental adjustment would prevail in this community with environmental filtering acting at large scales and across the elevation gradient, which would define the realized assemblages through functional and phylogenetic clustering, especially at higher elevations. By contrast, micro-environmental heterogeneity and species interactions would prevail at smaller scales, thereby leading to functional and phylogenetic divergence, which should be correlated with the soil characteristics at these scales. As a consequence, the functional heterogeneity should be relatively similar across the spatial scales studied because the primary environmental filters should stabilize the large functional heterogeneity between sites (regional scale), whereas micro-environmental heterogeneity and species interactions could force the opposite. Thus, in the present study, we addressed the following three questions: (i) Are the functional and phylogenetic clustering–overdispersion patterns within assemblages similar across spatial scales? (ii) Does the functional heterogeneity between assemblages vary among spatial scales? (iii) How do functional and phylogenetic patterns vary along the elevation gradient and with other key ecological factors (e.g., potential solar radiation, soil pH or nutrient heterogeneity, and species richness) at each spatial scale?

Materials and Methods

Study Site and Field Sampling

This study was conducted above the treeline at Sierra de Guadarrama National Park located northwest of Madrid in Spain (40°46′39′′ to 40°51′8″ N; 3°49′44″ to 4°4′59″ W; 1,940 to 2,419 m a.s.l.). The climate is Mediterranean, with a mean annual temperature and precipitation of 6.4°C and 1,350 mm, respectively, and a drought from May to October (less than 10% of the annual precipitation), which is more intense at lower elevations (Giménez-Benavides et al., 2007; Olano et al., 2013).

The treeline is located between 1,900 and 2,000 m a.s.l., and it is dominated by Scots pine (Pinus sylvestris L.) interspersed in a shrubby matrix of Cytisus oromediterraneus Rivas Mart. et al. and Juniperus communis L. subsp. alpina (Suter) Čelak. Above this elevation, the dominant shrubby matrix extends up to 2,100–2,200 m a.s.l., with dry and cryophilic grassland dominated by Festuca curvifolia Lag. ex Lange in less steep locations. In the highest areas (from 2,200 to 2,419 m a.s.l.), grassland represents the dominant vegetation and it is structured in a mosaic of patches dominated by F. curvifolia and some creeping and cushion chamaephytes clumped in a bare ground matrix (Pescador et al., 2014).

Thirty-nine sites covering the whole elevation range of F. curvifolia grasslands were selected in the summer of 2011 (Supplementary Figure 1A). At each site, we established a sampling plot of 20 m × 20 m and surveyed the vegetation at three different spatial scales: cell, subplot, and plot (Supplementary Figure 1B). First, a square measuring 2.4 m × 2.4 m was haphazardly located in the center of the plot and divided into 64 cells of 30 cm × 30 cm, which represented the cell scale (smallest scale). Second, five subplots of 2.4 m × 2.4 m (including the grid of cells) were established in the center and four corners of the plot as the intermediate scale. The coverage was visually estimated for every plant species, soil, and rock at the three spatial scales. At the plot scale, all plant species and their percentage cover were estimated as the average cover per species in the five subplots after also considering the presence and cover of the plant species not found in the cells or subplots but present in the plot (for more details, see López-Angulo et al., 2019).

Each plot was also characterized based on a set of environmental variables. In particular, we measured the elevation and orientation using a GPS system (Garmin Colorado-300) and the slope with a clinometer (Silva Clinomaster CM-360-%, LA). Orientation and slope values were employed to estimate the potential solar radiation according to Gandullo’s method (Gandullo, 1974; López-Angulo et al., 2018). Finally, 15 soil samples (diameter = 5 cm, depth = 10 cm) were randomly collected from each plot (Supplementary Figure 1B) under vegetated patches of the dominant grass (usually F. curvifolia; five samples), shrub canopies (C. oromediterraneus and J. communis subsp. alpina; five samples), and on bare ground (five samples). The soil samples were sieved (2-mm mesh) and air dried for 1 month, before the pH and soil organic carbon (SOC, %) were determined at the NutriLab-URJC laboratory. The pH was determined using a pH meter GLP 21 (Crison, Barcelona, Spain) and SOC was determined by colorimetry after oxidation with K₂Cr₂O₂ and H₂SO₄.

Plant Functional Traits

During the summer of 2011 and 2012, five plant functional traits were measured in a total of 56 species (Supplementary Table 1). These species represented more than 83% of the species in the community study and approximately 99% of the species cover. At least 10 individuals per species were selected and five functional traits were measured according to standardized protocols (Cornelissen et al., 2003): plant height (Hmax; distance between the ground and top of photosynthetic tissues), foliar thickness (FT), specific leaf area (SLA; ratio of leaf fresh area relative to dry weight), leaf dry matter content (LDMC; ratio of leaf dry weight relative to fresh weight), and seed mass (for more details, see Pescador et al., 2015).

DNA Isolation, Sequencing, and Phylogenetic Analysis

We collected fresh young leaves from three individuals per species and they were stored under dry conditions in silica gel for 1 month. Up to 20 mg of each dry leaf sample was disrupted and homogenized with a TissueLyser LT (Qiagen, Valencia, CA, United States) using glass beads. Total genomic DNA was extracted with a DNeasy Plant Mini-Kit (Qiagen, Valencia, CA, United States) according to the manufacturer’s instructions. Standard polymerase chain reactions (PCRs) were conducted to amplify two barcoding loci (i.e., rbcL and matK). In the PCR amplification mixture, 2 μl of DNA was added to 23 μl of reaction mixture comprising 2.5 μl of Taq buffer 2 mM with MgCl₂, 1 μl of dNTP Mix (0.4 mM), 1.25 μl of reverse and forward primers, and 1.25 U Taq DNA Polymerase (Biotools, Madrid, Spain). PCR was performed with an S1000 Thermal Cycler (Bio-Rad, United States) using the specific PCR conditions given in Supplementary Table 2. The PCR products were cleaned up using ExoSap-IT^® (USB Corporation, Cleveland, OH, United States) and submitted to Macrogen Inc. (Seoul, South Korea) for sequencing in the forward and reverse directions. In total, 288 sequences [three individuals per species × (56 species successfully sequenced with rbcL + 40 species successfully sequenced with matK)] were edited with Sequencher 4.1.4 (Genes Code Corporation), and the consensus sequences in each species were determined from the forward and reverse sequences. The rbcL and matK markers were aligned independently using Mafft online-version 7 (Katoh and Standley, 2013) following the G-INS-I strategy recommended for < 200 sequences with global homology. MacClade 4.06 (Maddison and Maddison, 2003) was used to concatenate primer configurations (rbcL + matK) and to create a combined matrix with missing data when the marker could not be sequenced.

The rbcL + matK combined matrix was used to reconstruct the phylogeny topology and the divergence times for our community using Bayesian inference and the Markov chain Monte Carlo (MCMC) approach in the BEAST 1.6.1 program (Drummond and Rambaut, 2007). We ran three MCMC chains with 50 million generations and sampled one tree and their parameters every 10,000 generations (partition). We estimate the rooted, time-measured phylogeny inferred using an uncorrelated relaxed molecular clock, which assumes that the substitution rates associated with each branch are drawn independently from a single lognormal distribution (Drummond et al., 2006). The substitution model used for each partition was GTR + Γ. The substitution model and estimated frequencies were unlinked among partitions. To calibrate the tree temporally, we constrained the divergence times for the minimum crown age at seven nodes based on the fossil record (Supplementary Table 3). Each constrained node was drawn following a lognormal or exponential prior with different parameters (for more details, see Supplementary Table 3). The Tracer 1.6 program (Rambaut et al., 2014) was used to evaluate the effective sample size for parameters, and we verified that three runs converged on the posterior distribution and reached stationarity. We combined the independent MCMC parameters and sample trees using the LogCombiner 1.6.1 program. Finally, all trees were summarized in a single maximum clade credibility tree (Supplementary Figure 2) after removing 2,500 trees as the burn-in and calculating the clade posterior probabilities and 95% honest posterior density intervals for node-specific parameters using TreeAnnotator 1.6.1.

Community Structure Metrics

We calculated the mean pairwise distance (MPD_obs; see de Bello et al., 2016) between the coexisting species observed in each cell, subplot, or plot to assess the functional and phylogenetic clustering–overdispersion patterns at different spatial scales. The functional and phylogenetic distance between each pair of coexisting species were calculated based on Gower’s dissimilarity matrix (Pavoine et al., 2009) using the information for the plant functional traits and the cophenetic correlation for the phylogenetic approach. Each MPD_obs value was expressed relative to the MPD simulated (MPD_exp) based on a null model to calculate a standardized effect size (SES; Gotelli and McCabe, 2002) for each scale, trait, all traits together, and phylogeny. Choosing an appropriate null model is a critical step when testing the effects of different mechanisms (Perronne et al., 2017), so we designed two “realistic” null models where the species names in each simulated assemblage were the names of the species in the plot to which the cell/subplot belonged (i.e., local pool, NM1), or the species names in all plots (i.e., regional pool, NM2). Thus, at the cell and subplot scales, the species that generated the null assemblage could come from the local or regional pool, whereas only the regional pool was considered at the plot scale. In addition, irrespective of the null model used, the number of species and total cover in each cell, subplot, or plot were kept fixed to control for differences in resource availability among scales and species (Gotelli, 2000; Le Bagousse-Pinguet et al., 2018). The likelihood that one species appeared in null assemblages was directly related to its abundance at the plot or regional scale, and the cover for each species in the null distribution was proportional to its abundance at the plot or regional scale. Based on these constraints, we generated 1,000 null assemblages for each cell (64 cells/plot × 39 plots), subplot (5 subplots/plot × 39 plots), and plot (39 plots), where the MPD_exp values were estimated using multi-trait, individual trait, and phylogenetic information, and each null model (i.e., NM1 and NM2). Comparing each MPD_obs to the corresponding 1,000 random MPD_exp values allowed us to calculate the SES of the MPD as follows:

where MPD_obs is the MPD obtained from the data observed at the pertinent working scale and using the functional or phylogenetic distance, MPDexp is the mean value for the 1,000 null MPD distributions, and sd(MPD_exp) is the standard deviation value for the 1,000 null MPD distributions. It should be noted that this index is equivalent to the inverse of the net relatedness index used in many phylogenetic approaches (Webb et al., 2002), and positive values indicate that the co-occurring species are more functional or phylogenetically distant than in the null assemblages (overdispersion pattern), whereas negative values denote that species are more closely related (clustering pattern). In total, we calculated SES–MPD values for 64 × 39 cells, 5 × 39 subplots, and 39 plots by considering each plant functional trait individually (i.e., Hmax, FT, SLA, LDMC, and seed mass), all traits together, and phylogenetic distances. At the cell and subplot scales, the null assemblage was built by considering the local (NM1) or regional (NM2) species pool. In each situation, the average SES–MPD value per plot was considered for each of the 64 cells and five subplots.

We estimated the coefficient of variation of the community weighted mean (CV-CWM; ratio of the standard deviation to the mean for average traits weighted by species abundances) for the realized assemblages at each spatial scale (i.e., cell, subplot, and plot) to characterize the functional heterogeneity at three spatial scales. We considered the information for 64 cells per plot to estimate the functional heterogeneity within each subplot (i.e., CV-CWM_w) because the summed area of the 64 cells was equivalent to the subplot area. We calculated the between-subplot heterogeneity by considering the five subplots per plot (i.e., CV-CWM_b). Finally, a unique CV was obtained for the 39 plots at the regional scale (i.e., CV-CWM_r). The latter variation component can be considered as the total variation between plots, and thus of the plant community in study area. All of these functional heterogeneity components were estimated by considering all traits together (in a multi-trait approach) and each trait independently (i.e., Hmax, FT, SLA, LDMC, and seed mass).

Statistical Analyses

We performed t-tests using the SES–MPD values at each spatial scale to determine whether the functional and phylogenetic diversities differed significantly from expectations. In particular, SES–MPD values significantly higher than 0 indicated overdispersion among assemblage species, whereas SES–MPD values lower than zero denoted clustering patterns. We performed paired t-test to evaluate the differences between functional heterogeneity within (CV-CWM_w) and between subplots (CV-CWM_b). We established linear models to assess the response of ecological factors on each SES–MPD following a multi-model inference (MMI) over all possible combinations of the factors and with AIC_c. A total of 35 independent linear models were drawn as a result of considering SES–MPD values for all traits together, each individual trait, and phylogeny at each of the three spatial scales and two null models. We used non-correlated ecological factors (i.e., elevation, potential solar radiation, pH, SOC, and species richness). We also considered the functional heterogeneity as an additional predictor factor. Within each independent model, all possible combinations of ecological factors were considered, and the resulting “partial models” were ranked according to the AIC_c value with respect to the “partial model” with the minimum AIC_c (ΔAIC_c = AIC_i – AIC_min). “Partial models” with ΔAIC_c < 2 were considered to be indistinguishable (Burnham and Anderson, 2002) and an average explanatory model was defined according to this threshold. The average estimate, significance, and relative importance of each predictor (w+) were estimated for each model. The relative importance of a predictor was evaluated by summing the w_i values for the partial models that contained the predictor of interest. Finally, the phylogenetic signal was calculated for each functional trait with Blomberg’s K index to capture the effect of trait evolution and to assess the existence of phylogenetic niche conservatism in our system. Basically, Blomberg’s K represents the ratio of the observed phylogenetically correct mean-square error divided by the mean-square error using a variance matrix derived from the phylogeny standardized by the expectation under a Brownian motion model (Blomberg et al., 2003). The significance of K (p-value) was calculated by comparing it to the null distribution with 9,999 replicates. Thus, K < 1 indicates that closely related species resemble each other less than expected under the Brownian motion model of trait evolution (no phylogenetic signal), whereas K > 1 suggests that closely related species are more similar than predicted by the model (stronger phylogenetic signal).

All statistical analyses were conducted with R 3.6.0 (R Core Team, 2019) and the following functions: dbFD and gowdis functions implemented in the FD package (Laliberté and Legendre, 2010; Laliberté and Shipley, 2011), dredge and importance functions in the MuMIn package (Kamil, 2013), melodic function (de Bello et al., 2016), and phylosignal function in the picante package (Kembel et al., 2010).

Results

Overall, we identified significant deviations in the functional and phylogenetic clustering–overdispersion patterns depending on the spatial scale (i.e., cell, subplot, or plot scale) and null model (i.e., species consideration from the plot pool, NM1, or from the regional pool, NM2; Figure 1 and Supplementary Table 4). At the plot scale, the trait clustering pattern (SES–MPD values < 0) was most relevant for all single-trait and multi-trait approaches. At the subplot and cell spatial scales, significant trait overdispersion (SES–MPD values > 0) had considerable importance for most traits irrespective of whether we considered the plot (NM1) or regional pool species (NM2). By contrast, Hmax had a significant clustering pattern at both spatial scales and when considering both null models (Figure 1B). The community had a weak phylogenetic clustering pattern at the plot scale (species at plot scales were more closely related than expected from the null model), and significant clustering was found at the subplot and cell scales when both null models were used (Figure 1G).

Figure 1. Density plots of standardized effect size of mean pairwise distance (SES–MPD) for (A) multi-trait, (B–F) each individual plant functional trait and (G) phylogeny. Three different operating spatial scales (cell, subplot, or plot scale) and two null models (simulated assemblages using plot pool, NM1, and regional pool, NM2) are included in each plot. Negative values indicate that co-occurring species within the assemblages are more related than expected by chance (i.e., clustering), whereas positive values denote distances greater than expected by chance (i.e., overdispersion). Black-dashed arrows represent the SES threshold of 0. Each panel shows the results of one-sample t-tests conducted to determine whether SES–MPD values for each approach, null model, and scale were positive (+) or negative (–), and significantly different from 0 as: ^∗∗∗p < 0.001; ^∗∗p < 0.01; ^∗p < 0.05. Full scores for the one-sample t-tests are presented in Supplementary Table 4. Hmax, plant height; SLA, specific leaf area; LDMC, leaf dry matter content; FT, foliar thickness.

The functional heterogeneity across spatial scales exhibited similar patterns for the multi-trait approach and based on each plant functional trait (Figure 2 and Supplementary Table 5). On average, the functional heterogeneity was usually greater at the regional scale (i.e., CV-CWM_r) than those within and between subplots (i.e., CV-CWM_w and CV-CWM_b, respectively), except for traits related to the leaf economic spectrum (i.e., FT, SLA, and LDMC), where the values of CV-CWM_w were similar or even higher than those of CV-CWM_r. Paired t-tests between the CVs at the plot level (i.e., CV-CWM_w and CV-CWM_b) indicated lower variability between subplots than within subplots in most cases, except for Hmax (Figure 2B) and FT (Figure 2C), where both components were not significantly different (p-values according to paired t-test = 0.50 and 0.09, respectively; Supplementary Table 5).

Figure 2. Relationships among coefficient of variation (CV) of community weighted mean (CWM) at the three spatial scales assessed [within subplots (CV-CWM_w; y-axis in each panel), between subplots (CV-CWM_b; x-axis in each panel), and between plots (CV-CWM_r; solid lines in each panel)] and considering (A) multi-trait approach and (B–F) each individual plant functional trait. Dashed lines represent a 1:1 relationship. In each panel, p-values obtained from paired t-tests for within subplots (w), between subplots (b), and between plots (r) are denoted as: ***p < 0.001; **p < 0.01; *p < 0.05. Full results for these paired t-tests are presented in Supplementary Table 5. Hmax, plant height; FT, foliar thickness; SLA, specific leaf area; LDMC, leaf dry matter content.

The effects of ecological factors on the SES–MPD values depended on the spatial scale and the null model considered (Figure 3). In particular, clustering patterns (negative relationships) were detected for multi-trait, SLA, and seed mass as the elevation increased at the plot scale and under NM2. By contrast, a phylogenetic overdispersion pattern (positive relationship) was more obvious at the smallest spatial scale and under NM2 as the elevation increased. The potential solar radiation effect led to clustering patterns for multi-trait and seed mass for co-existing species at the subplot scale. Among the soil factors, an increase in SOC resulted in the overdispersion of SLA at different scales but a clustering pattern of FT at the plot and subplot scales and under NM2. Higher soil pH levels resulted in overdispersion patterns for Hmax and seed mass at the subplot and cell scale, respectively. By contrast, clustering patterns were found for FT, SLA, and LDMC when the soil pH levels increased at the subplot scale and under NM1. Increased species richness was related to clustering patterns for multi-trait, Hmax, FT, and phylogeny at the three spatial scales. Finally, overdispersion patterns for Hmax and FT at the cell scale were directly dependent on the functional heterogeneity characterized as the CV of the corresponding CWM.

Figure 3. Linear averaged models showing the responses of non-correlated ecological factors on each standardized effect size of mean pairwise distance (SES–MPD) for multi-trait, individual functional traits, and phylogenetic approaches. SES–MPD values were considered at three spatial scales (i.e., cell, subplot, or plot scale) and under two null models (i.e., simulated assemblages using plot pool, NM1, and regional pool, NM2). In each case, all possible predictor combinations were tested and “partial models” were ranked according to Akaike’s information criterion (AIC_c). Averaged models were drawn from the best partial models based on ΔAIC_c < 2 (ΔAIC_c = AIC_i – AIC_min). Predictors included in each averaged model are denoted by a dot (.) and positive or negative significant relationships (according to the average-slope value) are highlighted in light gray or dark gray, respectively. R² (mean ± standard deviation) and AIC_c values are shown in each case for the best partial models. Full results for the averaged linear models are shown in Supplementary Table 6. MT, multi-trait; Hmax, plant height; FT, foliar thickness; SLA, specific leaf area; LDMC, leaf dry matter content; SM, seed mass; Phy, phylogeny; Elev, elevation (m a.s.l.); PSR, potential solar radiation; SOC, soil organic carbon (%); pH, soil pH; Sp, number of species per plot; FH, functional heterogeneity (CV-CWM).

Discussion

Our results indicate the existence of three assembly mechanisms operating at different spatial scales in the same plant community type on a Mediterranean high mountain. First, we established the scale dependence of the mechanisms involved in the assembly process ranging from functional clustering at regional spatial scales to functional overdispersion at the smallest scale (Mason et al., 2011; de Bello et al., 2013b; Carboni et al., 2014), and phylogenetic clustering irrespective of the spatial scale. Second, we found that the functional heterogeneity at the smallest scale represented a source of variation, which was as important as the heterogeneity detected along the whole elevation gradient, especially for the traits related to the leaf economic spectrum. Third, we showed that different ecological factors (i.e., elevation, potential solar radiation, SOC, pH, species richness, and functional heterogeneity) regulated the functional and phylogenetic patterns among co-occurring species and ultimately determined the community assembly. Overall, the results obtained in our study confirmed the importance of considering approaches based on functional traits and phylogenetic information across spatial scales in plant coexistence assessments (Kraft and Ackerly, 2010; Perronne et al., 2014).

Diversity Patterns Across Spatial Scales

By assessing the SES–MPD according to two null models (simulated community assemblages using plot pool, NM1, and regional pool, NM2) and at three spatial scales (cell, subplot, or plot scale), we detected contrasting clustering–overdispersion patterns in the functional and phylogenetic diversity (Freschet et al., 2011; Bernard-Verdier et al., 2012; de Bello et al., 2013b). We found that an increase in the spatial granularity (i.e., plot scale of 20 m × 20 m) implied a clustering pattern for most plant functional traits, thereby supporting the idea that the effects of environmental filters are more conspicuous at larger scales (Weiher and Keddy, 1995). By contrast, when the spatial granularity was reduced (i.e., subplot and cell scales), functional overdispersion was the most frequent pattern for all functional traits, except for Hmax. Similarly, some previous studies found an overdispersion pattern at small spatial scales (Stubbs and Wilson, 2004; de Bello et al., 2013b), and our results support the limiting similarity hypothesis for most functional traits at these scales. Hmax was the only functional trait that exhibited a clustering pattern at these scales (Figure 1B), and thus the heights of species in a subplot or cell tended to be similar. In high mountain areas, cold is a strong abiotic filter and greater height may be considered a limitation for the survival of species (Körner, 2003). Consequently, low stature appears to be the most suitable strategy for this type of environment. Furthermore, Mayfield and Levine (2010) suggested the existence of similarity patterns in Hmax linked to competitive ability. The phylogenetic relationship assessment detected clustering patterns at three study scales, thereby indicating that the species were more closely related than the null expectation. This pattern is consistent with the prevalence of environmental filtering indicated for most of our functional traits and the idea that phylogenetic relatedness suggests functional proximity (Webb, 2000; Cavender-Bares et al., 2006; Kraft and Ackerly, 2010). In fact, environmental filtering could be identified with phylogeny as a proxy of functional differences among species at larger scales, but this pattern might have been caused by competitive exclusion leading to clustering after the elimination of the weakest competitors (Mayfield and Levine, 2010). It should be noted that this pattern was more obvious with NM1, which implies that the species pool selected passed a specific environmental filter, and thus the patterns detected could probably be attributed to the incidence of biotic factors (species competition for the same resources) rather than environmental heterogeneity (Stubbs and Wilson, 2004).

Functional Heterogeneity Across Spatial Scales

Overall, functional heterogeneity at the smallest scale (CV-CWM_w) for the traits related to the leaf economic spectrum showed a variation similar to or even higher than that observed at the regional scale (CV-CWM_r). These findings suggest the existence of (a) a primary abiotic filter at the regional scale determining low functional turnover throughout the distribution area and minor differences between plots (i.e., low CV-CWM_r); (b) the existence of local environmental heterogeneity acting within each plot to yield functional heterogeneity (i.e., plots with relatively similar CV-CWM_b values rather than CV-CWM_r values); and (c) different and significant determinist processes at the smallest scale, which would also have been conditioning the species assembly (i.e., similar or higher CV-CWM_w values to those between plots). As noted above, cold and summer drought are primary filters in Mediterranean high mountains that promote the existence of a similar and well-adapted species pool among sites (i.e., low composition and functional turnover along elevation; Helmus and Ives, 2012). This relatively limited pool of species (no more than 60 species along a 500-m elevation gradient in Sierra de Guadarrama National Park) may be subject to the effects of other sources of environmental heterogeneity (potential solar radiation, light competition, soil characteristics, and shrub encroachment) or simply to stochasticity. These sources could constrain the species assemblage in each plot and lead to functional heterogeneity between plots (i.e., CV-CWM_r), which is also limited by the primary filters. Nevertheless, the ecological variables that lead to other sources of heterogeneity are relatively similar between the subplots in each plot (same elevation, potential solar radiation, or soil characteristics), thereby resulting in redundant species assemblages between subplots with similar functional heterogeneity (i.e., CV-CWM_b). By contrast, at small spatial scales, the species in a plot are subject to additional assembly mechanisms based on micro-environmental gradients (Hutchings et al., 2003; Burton et al., 2011) or species interactions (Pescador et al., 2014), which may increase the functional heterogeneity at the cell scale. These additional assembly mechanisms lead to a scenario where the deterministic processes that operate at small scales are as important for generating functional variability as those observed at larger scales. In general, if the functional heterogeneity at the smallest scale is in the range of that at the regional scale, then any model that attempts to explore the species distribution and community composition based on the environmental conditions and species functional traits (see Shipley et al., 2006) should consider at least these sources of spatial variability.

Relevance of Ecological Factors on the Diversity Patterns

Linear averaged models showed that assembly processes were influenced by ecological factors that depended on the scale, functional traits, and null model considered (Mason et al., 2011; Spasojevic and Suding, 2012; de Bello et al., 2013a; Price et al., 2014; Reitalu et al., 2014). Species in plots located at the highest elevations exhibited a clustering pattern according to the multi-trait approach as well as for the SLA and seed mass traits. However, when the sampling granularity decreased to the cell level, phylogenetic overdispersion increased with the elevation but with no effects on functional responses. These results suggest that species at the plot level use a similar functional strategy to cope with the harsh conditions on the summits (Weiher and Keddy, 1995; Stubbs and Wilson, 2004; Schwilk and Ackerly, 2005; Chalmandrier et al., 2017). However, a phylogenetic overdispersion pattern is needed at the smallest scale to reduce competition with direct neighbors for the same local resources (Cavender-Bares et al., 2004, 2006; Pescador et al., 2014). Phylogenetic patterns have been investigated previously along elevation gradients with contrasting findings (Bryant et al., 2008; Culmsee and Leuschner, 2013; González-Caro et al., 2014; Li et al., 2014; Qian et al., 2014). In particular, similarly to our results, Bryant et al. (2008) found that plant communities tended to exhibit phylogenetic overdispersion at higher elevations, probably due to the prevalence and intensification of facilitation with elevation, whereas competition was more important at lower elevations. Furthermore, high levels of stress can lead to more open and discontinuous vegetation cover providing more gaps for colonization, and thus higher phylogenetic diversity (Reitalu et al., 2014). Based on the conceptual framework of functional and phylogenetic patterns (Webb et al., 2002; Cavender-Bares et al., 2004; Kraft and Ackerly, 2010), our findings imply that important traits for habitat specialization in harsh conditions probably exhibit convergent evolution (i.e., a low phylogenetic signal; Supplementary Table 7).

The absence of responses in terms of the functional and phylogenetic diversity to potential solar radiation at the plot and cell scales suggest that difference in the total annual radiation between north- and south-facing slopes was less important as a filter (Méndez-Toribio et al., 2017). However, multi-trait and seed clustering was detected at the subplot scale when the potential solar radiation values were higher. Thus, at this scale, co-occurring species must employ similar physiological and regenerative strategies to cope with water stress in plots with high potential solar radiation, especially in the Mediterranean-type climate. Irrespective of the spatial scale considered, the multiple effects of soil factors on functional patterns suggest high heterogeneity in the soil characteristics at small scales (Jackson and Caldwell, 1993; Hutchings et al., 2003; Gazol et al., 2013), especially given the strong effects of soil pH (Chytrý et al., 2003) and SOC in grasslands. Remarkably, we found clustering patterns for the traits related to the leaf economic spectrum and overdispersion patterns for Hmax and seed mass at higher pH values, and these patterns have not been described previously to the best of our knowledge. Furthermore, the prevalence of clustering as the species richness increases was found previously using functional (Carboni et al., 2014; Niu et al., 2014) and phylogenetic approaches (Reitalu et al., 2014). Thus, if the environment favors the co-occurrence of a great number of species, they will tend to be functionally and phylogenetically more redundant and similar (Cadotte et al., 2011), possibly because the regional species pool is unable to provide new candidates with different functional strategies to cope with harsh conditions (Yachi and Loreau, 2007; Pakeman, 2011). As a consequence, this predictor together with elevation had hardly any effect on the ecological clustering–overdispersion patterns when NM1 was used (i.e., the species pool was represented by the species present in the plot considered), probably because the species in null assemblages were from a species pool subject to the same elevation and species richness as the realized assemblage. Finally, the overdispersion patterns for the Hmax and FT traits were modulated by an increase in the functional heterogeneity at the smallest scale (CV-CWM_w), and thus site-to-site functional variation should be considered as a relevant factor when assessing community assemblies.

Conclusion

Our multi-scale approach applied to a single plant community type with a low species turnover highlights the importance of evaluating functional and phylogenetic patterns at different spatial scales to understand the current community assembly. We showed that complex functional and phylogenetic patterns emerge depending on the spatial scale considered. Environmental filtering was found at larger spatial scales where species coexist under particular conditions, whereas functional overdispersion was present at small scales, thereby suggesting secondary niche differentiation after the effects of primary filters such as cold and summer drought. Phylogenetically based analysis indicated a pattern consistent with the functional analysis at larger scales, but this congruence was absent at the subplot and cell scales where the phylogenetic pattern remained clustered. Our results also demonstrated the importance of functional heterogeneity at the smallest scale in our plant community, thereby confirming the presence of primary environmental filters and suggesting the existence of additional assembly mechanisms based on micro-environmental gradients and species interactions at the smallest scale. Finally, our results show that functional and phylogenetic diversity patterns differ in terms of their responses to ecological factors across spatial scales, which highlights the need to develop adequate tools for predicting the response of vegetation to global warming (Reu et al., 2011) or the distributions of species based on the relationships between environmental conditions and species functional traits (Shipley et al., 2006). Given the shifts in the functional and phylogenetic responses along the elevation and among spatial scales, both sources of diversity should be considered by employing a multi-scale approach to help understand the future of alpine biodiversity under ongoing warming.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, MW927222–MW927321. Data sets on vegetation co-occurence patterns, environmental variables, and functional traits for each plant species can be found in López-Angulo et al. (2019).

Author Contributions

DP, AE, and FV conceived and designed the study. DP and JL-A conducted the fieldwork. DP and FB performed the statistical analysis. DP wrote the manuscript. All authors contributed to the discussion and interpretation of data, revised the manuscript, and approved the submitted version.

Funding

This research was supported by the Spanish Ministry of Economy and Competitiveness (PHENOTYPES project, PGC2018-099115-B-100), Community of Madrid (REMEDINAL TE-CM project, S2018/EMT-4338), and Grant Agency of the Czech Republic (P505/12/1296).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank Carlos Díaz Palomo, Felipe García García, and Jan Lepš for their technical and methodological assistance in this work. We would also like to thank Ana Millanes for help with phylogenetic reconstruction, Lars Markesteijn and Duncan E. Jackson for English language editing assistance, and staff at the Parque Nacional de la Sierra de Guadarrama for providing permission to work in the field area.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.622148/full#supplementary-material

References