Frequency and intensity of facilitation reveal opposing patterns along a stress gradient

Abstract Disentangling the different processes structuring ecological communities is a long‐standing challenge. In species‐rich ecosystems, most emphasis has so far been given to environmental filtering and competition processes, while facilitative interactions between species remain insufficiently studied. Here, we propose an analysis framework that not only allows for identifying pairs of facilitating and facilitated species, but also estimates the strength of facilitation and its variation along environmental gradients. Our framework combines the analysis of both co‐occurrence and co‐abundance patterns using a moving window approach along environmental gradients to control for potentially confounding effects of environmental filtering in the co‐abundance analysis. We first validate our new approach against community assembly simulations, and exemplify its potential on a large 1,134 plant community plots dataset. Our results generally show that facilitation intensity was strongest under cold stress, whereas the proportion of facilitating and facilitated species was higher under drought stress. Moreover, the functional distance between individual facilitated species and their facilitating species significantly changed along the temperature–moisture gradient, and seemed to influence facilitation intensity, although no general positive or general negative trend was discernible among species. The main advantages of our robust framework are as follows: It enables detecting facilitating and facilitated species in species‐rich systems, and it allows identifying the directionality and intensity of facilitation in species pairs as well as its variation across long environmental gradients. It thus opens numerous opportunities for incorporating functional (and phylogenetic) information in the analysis of facilitation patterns. Our case study indicated high complexity in facilitative interactions across the stress gradient and revealed new evidence that facilitation, similarly to competition, can operate between functionally similar and dissimilar species. Extending the analyses to other taxa and ecosystems will foster our understanding how complex interspecific interactions promote biodiversity.


| INTRODUCTION
As the rise of biogeography, researchers have sought to understand how plant-plant interactions change along environmental gradients, and what consequences this has for the composition of plant communities (e.g., von Humboldt & Bonpland, 1807). Two types of interactions are dominant in shaping community composition: competition and facilitation (Brooker & Callaghan, 1998). Competitive or facilitative interactions can be defined as interactions in which the presence of one species alters the environment (or occupies space) in a way that reduces or enhances growth, survival, and reproduction of a second species (Bronstein, 2009;Craine, Fargione, & Sugita, 2005;McIntire & Fajardo, 2014). The relative importance of these two processes has been shown to vary along environmental gradients, with competition generally dominating in communities of low-abiotic stress, while facilitation increases in importance with abiotic stress (framed in the stress gradient hypothesis; Bertness & Callaway, 1994;Callaway & Walker, 1997;Choler, Michalet, & Callaway, 2001;Callaway et al., 2002;Michalet, Schöb, Lortie, Brooker, & Callaway, 2014).
Previous work on facilitative interactions has repeatedly demonstrated that facilitation can act as a major force structuring plant communities, and helped identifying its putative underlying mechanisms (McIntire & Fajardo, 2014). Nonetheless, our understanding of this process remains limited. On the one hand, even though facilitation is usually thought to be more important under stressful conditions (Callaway, 2007), it may not necessarily be restricted to stressful conditions only (Holmgren & Scheffer, 2010;McIntire & Fajardo, 2014).
Indeed, it may happen that only few keystone species provide important facilitative services to many facilitated species under stressful conditions, while under less stressful conditions facilitation may be of lower intensity but provided by a larger number of species. Such a situation may explain why signals of facilitation are often lost under environmental conditions that are favorable to plant growth. Therefore, to be able to capture the full extent of facilitative interactions, we need to develop a community-level understanding of how facilitation varies, both in intensity and frequency, along large environmental stress gradients. Key components of such a community-level assessment should include both, the identification of each facilitating-facilitated species pair in the communities, and the estimation of the degree to which these facilitative interactions contribute to the increase in fitness of facilitated species (hereafter called facilitation intensity; Welden & Slauson, 1986).
On the other hand, the mechanisms determining the nature and magnitude of facilitation remain poorly understood (McIntire & Fajardo, 2014;Schöb, Butterfield, & Pugnaire, 2012). Facilitation mechanisms can be symmetric or asymmetric, and involve direct or indirect drivers. Asymmetric facilitation indicates that one species (the benefactor or facilitating species) will disproportionately favor another species (the beneficiary or facilitated species) more than it can mutually profit from this species. For instance, tall plants may protect shorter plants from ultraviolet radiations (asymmetric facilitation or commensalism, as shorter plants do not protect taller plants from radiations), while species with similar flower color may attract the same pollinators (symmetric facilitation or mutualism; Brooker & Callaghan, 1998;Chu et al., 2009;Lin, Berger, Grimm, & Ji, 2012). Within a functional framework, if species facilitate each other via the same mechanism (e.g., pollinator attraction via similar flower color, or soil stabilization via root reinforcement) then the intensity of symmetric facilitation should increase with species functional similarities, whereas that of asymmetric facilitation should increase with functional dissimilarities (Butterfield & Briggs, 2011;Cavieres & Badano, 2009;Gross et al., 2009 are expected to relate to their facilitation intensity. Therefore, a key challenge today is to quantify the relationship between facilitation (a) symmetry, intensity, and species functional (dis)similarities.
These knowledge gaps-about the relationship between facilitation (a)symmetry, intensity, and species functional (dis)similarities-are not due to a lack of experiments or observational studies, but for the large part rather due to a lack of methodological approaches allowing for investigations of large environmental gradients and of species-rich communities, where multispecies interactions are not known a priori and where indirect interactions may be frequent (such as intransitive competition; Gallien, 2017;Gallien, Zimmermann, Levine, & Adler, 2017). Indeed, most studies on facilitation mechanisms to date have relied on: (1) the comparison of communities in paired plots containing or not the facilitating species (e.g., species growing inside versus outside of a cushion plant; Butterfield et al., 2013); (2) experiments testing the effect of removing the facilitating species (e.g., Albrecht et al., 2015;Callaway et al., 2002;Cipriotti & Aguiar, 2015;Michalet et al., 2015); or (3) monitoring long-term changes in community composition during primary succession (e.g., Martorell & Freckleton, 2014).
These approaches are all valuable, yet they strongly rely on a priori knowledge about the facilitating species and/or extensive monitoring efforts. There is, thus, a strong need for screening methods based on comparably simple data, which allow for analyzing multispecies interaction links without experiments.
Here, we propose and apply a simple but robust framework for exploring facilitation patterns without a priori information on the local species and the processes that drive species co-occurrences, and without need for experimental manipulation. This screening procedure allows for identifying pairwise facilitative interactions in species-rich communities and for tracking their variation along large environmental gradients. We use the output of this approach to specifically investigate the relationships between facilitation intensity and species functional (dis)similarities along a long stress gradient using a large community dataset. This helps us progressing toward a better understanding of the facilitation process in plant communities and toward designing more complex and targeted experiments.
We first describe our proposed approach and evaluate its performance using a community assembly simulation model (VirtualCom; Münkemüller & Gallien, 2015). As these simulations show that our approach works well and facilitation is accurately detected, we are confident to apply it on a large dataset of 1,134 plant community plots from the Zermatt region (Switzerland) and tackle key questions related to facilitative interactions. Specifically, we asked the following: (1) How do facilitation frequency and intensity change along environmental gradients? (2) Does the functional distance between the facilitating and facilitated species change along environmental gradients? (3) Is facilitation intensity influenced by the functional similarity between the involved species? Finally, we discuss future avenues and potential research questions that can be answered using our approach.

| The facilitation screening procedure
We propose a screening framework that elaborates on the widely studied co-occurrence patterns (e.g., Boulangeat, Gravel, & Thuiller, 2012;Diamond, 1975;Jackson, Somers, & Harvey, 1989;López, Valdivia, Rivera, & Rios, 2013;Ulrich & Gotelli, 2007), and thus, only requires community relevés with recordings of the relative cover of coexisting species (i.e., at a relatively small grain size at which species interact). By combining co-occurrence analyses with analyses of co-abundance, we aim at detecting facilitation for species pairs within a specific environment. Our method estimates for each pair of cooccurring species A and B, whether species A facilitates species B and by which intensity.
To avoid confusion with environmental filtering signals, our method groups community relevés into ecologically very narrow bins of similar environmental conditions (Figure 1, steps 2). Within each bin, we then identify facilitating-facilitated species pairs by testing all possible species pairs (Figure 1, steps 3). For each species pair A and B, species A is considered as facilitating species B if it fulfills the following F I G U R E 1 The six major steps of the proposed analysis framework. Once the community relevés have been sampled (step 1), the main environmental gradient(s) among them shall be identified (for instance with a principal component analysis), and then the communities are grouped into "bins" of similar environmental conditions (step 2). Next, for each bin, all possible species pairs are tested for facilitative interactions (see also Table 1), and the facilitation intensity is estimated as the difference in abundance of the facilitated species when the facilitating species is present versus absent (step 3). The performance of the applied methodology is then evaluated on artificial communities simulated with different assembly rules using the VirtualCom simulation model (step 4). After this preliminary test, we calculate a number of facilitation metrics within each bin (such as the number of facilitation links or the average facilitation intensity received by the facilitated species), and analyze how they change along the stress gradients (for instance with regression models; step 5). Finally, using functional trait information one can test for each pair of species the relationship between the facilitation intensity, the functional distance, and the environmental gradient (step 6) per bin

A facilitates B if:
(1) Abundance of A is higher than abundance of B Abundance A as threshold (2) B co-occurs with A

Randomization test (3) B is absent when A is absent
Randomization test (4) abundance of B is higher when A present ANOVA; facilitation intensity is given by increase in abundance of B Step 3: Identify facilitation pairs and intensity for each bin Step 5: Calculate summary statistics per bins Step 6 Step 2: Group communities in bins Step 1: Sampling

Sampled sites
Functional dissimilarity Facilittation intensity

?
Step 4: Method evaluation with VirtualCom (1) The relative cover of A is higher than the relative cover of B, meaning that we assume a facilitator to have higher plant cover than the facilitated species, (2) B co-occurs with A more often than expected at random, and (3) B is more often absent when A is absent than expected at random. These two latter requirements were tested for significance via randomization tests where species B occurrences were permutated among communities independently of the presence of species A, with 499 randomizations per species (using a .025 significance threshold). Finally, (4) the relative cover of B is significantly higher in community relevés where A is present than when A is absent (Table 1). This was assessed using ANOVA tests. When significant, the amount of increase in relative cover of species B (when A is present vs. absent) was used as an estimate of facilitation intensity received by B (Figure 1, step 3). In other words, this framework allows for testing whether a facilitated species benefits from a facilitator more than can be expected by chance, both regarding its presence and its abundance: The presence of the facilitating species A increases both the likelihood of occurrence and the abundance of the species B, while A is not necessarily affected by B (Table 1 and Figure 1 step 3). Note that constraining the relative cover of the facilitator species to be higher than the one of the facilitated species generally brings a stronger focus on asymmetric facilitation patterns (e.g., facilitation via shading), but this constrain could be loosened to integrate symmetric facilitation (e.g., facilitation via pollinator attraction).
Our method relies on two fundamental assumptions: (1) The environmental heterogeneity among the considered communities is negligibly small, and (2) the within site microhabitat heterogeneity is negligible. Indeed, if environmental heterogeneity is too high, finescale environmental filtering processes may lead to differences in cooccurrence and co-abundance patterns similar to those expected from facilitation (i.e., if the niche of the species B is nested within the one of A and A's abundance is generally higher than that of B). We note that such assumption about environmental homogeneity is similarly made (although not always explicitly) in most analyses of community functional similarity patterns (e.g., when inferring environmental filtering and competition processes; Münkemüller et al., 2014;Willis et al., 2010). Additionally, our estimation of facilitation intensity relies on the assumption that an increase in relative cover of species is associated with an increase in its fitness. Although this assumption is likely to be verified in most situations, some systems might present exceptions that would preclude the utilization of our methodology.

| Method validation with processed-based community assembly simulations
As proof of concept that our approach is capable of reliably detecting facilitation and that it does not confound facilitation with other coexistence mechanisms (e.g., environmental filtering, competitive interactions, and neutral mechanisms), we used a virtual ecologist approach (Gallien, Carboni, & Münkemüller, 2014;Zurell et al., 2010), and compared four different simulation scenarios: facilitation, environmental filtering, competitive interactions, and neutral coexistence. To do so, co-occurrence patterns were generated using the recently published community assembly model VirtualCom (Münkemüller & Gallien, 2015). VirtualCom has originally been developed to simulate commu- The abundance of A is higher than that of B?
A and B co-occur more than by chance?
Green ticks indicate significant positive responses to the identification rules, while red crosses indicate significant negative responses (ticks and crosses are represented together when both responses are possible).
repetitions × 50 community per repetition = 17,500 communities overall). For each repetition, we assessed the false-positive (proportion of pairs identified as facilitating while they were not) and falsenegative (proportion of "true" facilitating pairs not identified by our method) error rates.

| Method application with the Zermatt dataset
We then used this new screening method to detect and quantify facilitation in empirical data, using phytosociological relevés of ca.
2 m × 2 m in the Zermatt mountain region in Switzerland, composed of 1,242 plots sampled in natural and seminatural vegetation during the 1990s by several persons and summarized in Steiner (2002).
The sampling covers an elevation gradient ranging from 1,536 m to 3,390 m a.s.l. (Appendix S1: Figure S1). When cleaning the dataset, we identified 108 sites containing species typical of very wet habitats indicating local water sources independent of climatic humidity gradients. We removed them from the dataset to avoid potential confounding effects of mixing different habitat types and microhabitat heterogeneity (which left us with 1,134 sites). Overall, the dataset contained a total of 574 species. Within each community plot, species relative cover was recorded using the Braun-Blanquet cover scheme (see Appendix S1 for more details). In order to avoid statistical errors due to low sample size, we chose to work with those species that were present in at least 20 community plots, which left us with 262 species for further analyses (representing 87% of the vegetation cover on average).  (6) the topographic wetness index (following Beven & Kirkby, 1979). These variables are considered to have direct physiological effects on species distributions and were used in many previous studies successfully (e.g., Randin et al., 2006;Zimmermann & Kienast, 1999). All variables were available at a 25 m spatial resolution, which is of fine enough grain to match the 2 × 2 m resolution of the community plots. The uncertainty in the temperature and precipitation data is summarized in Zimmermann and Kienast (1999), and is small enough to not confound the results along this steep and climatically very long gradient. We then chose the first PCA axis as representative of the stress gradient among sites for all further analyses because it revealed a warm/dry-to-cold/wet gradient (representing 57% of the intersite environmental differences, Figure 2a).

| Sampling along environmental gradients
Note that we used indicator species to remove sites that indicated local water sources independent of the climatic humidity gradient as

| Species-level functional traits
To investigate differences in functional similarity between the facilitating and facilitated species along the studied environmental gradient and with changing facilitation intensity, we used six species-specific, functional traits. These traits relate to the species' microhabitat preferences and life history strategies (available in Flora Indicativa; Landolt, 2010), and thus to facilitation. The traits related to microhabitat preferences included species preferences for light availability, soil moisture level, humus level, and soil aeration; traits related to species life history strategy included species average leaf life span and CRS life strategies as defined by Grime (2001). For 10 species, these traits were not available, and therefore, all functional analyses were run on 252 instead of 262 species.

| Statistical analyses
To answer our initial three questions, we followed three consecutive analytical steps. First, we tested for general trends in the frequency and intensity of facilitation along the environmental gradient (the PCA 1st axis). Second, we characterized the level of environmental stress (hereafter called environmental filtering) along the environmental gradient using functional diversity indices (Figure 1, step 5). This allowed us to then quantify the changes in functional dissimilarity among the facilitating (and among facilitated) species along the stress gradient.
Third, for each facilitated species, we investigated whether: (1) the facilitation intensity received changed along the gradient, (2) the functional distance to the facilitators changed along the gradient, and (3) whether the facilitation intensity received by a facilitating species could be related to the functional distance to its facilitators (Figure 1, step 6).

| Functional trends along the environmental gradients
In order to estimate the intensity of environmental stress along our warm/dry-to-cold/wet gradient (PCA axis 1), we calculated the mean functional distance (MFD) between all pairs of species within each community from our set of six traits (using the Gower distance that can handle both continuous and categorical variables; Gower, 1971).
Thereby, we expected that the stronger the environmental stress the more functionally similar is the coexisting species (compared to the full set of species in the dataset), as they should have similar traits to cope with the environmentally stressful conditions (Weiher & Keddy, 1999). We used the MFD standardized effect size (hereafter called MFD SES ) to estimate the strength of this environmental filter in each bin. MFD SES was obtained from null models by randomizing the functional distances among species, and thus by controlling for the community richness (999 repetitions). MFD SES varies between 0 (perfectly similar species) and 1 (completely dissimilar species; details in Appendix S1). We then tested whether the MFD SES scores changed along the environmental gradient using GLMs with linear and/or quadratic relationship followed by stepwise AIC variable selections.
Analogous to the procedure outlined above, we calculated func-

| Linking facilitation intensity with functional information
At the species level, we further investigated (1)

| A new approach to detect facilitative interactions from community data
Our proposed approach to detect pairs of facilitating and facilitated species revealed great performance in tests using simulated data.
First, our approach did not detect any facilitating interactions when there were none; that is when we simulated community assembly with scenarios of: environmental filtering (rate of false positive = 0), competition (rate of false positive = 0), and neutral coexistence (rate of false positive = 0). Second, under facilitation scenarios, we could identify the correct facilitating-facilitated species pairs, although the rates of false positives (i.e., species were wrongly identified as facilitating or facilitated while they were not) and false negatives (i.e., facilitating or facilitated species were not detected) were not null (Appendix S1: Figure S2). However, the false-negative error rates were generally low (mean error rates < 0.05 for all scenarios) and decreased when the simulated number of facilitated species increased. The false-positive error rates were very low (error rates < 0.01 for all scenarios), although they slightly increased when the number of facilitated species increased (Appendix S1: Figure S2). Overall, these error rates indicate that our test is generally conservative, especially when many species are facilitated in the studied system, meaning that some true facilitation pairs may be overlooked, but the probability of falsely identifying a species as facilitating is less than 1%. Thus, our approach is able to correctly identify facilitating and facilitated species pairs, given that the underlying assumptions are met (i.e., the environmental heterogeneity within the bins and within the communities is negligible).

| Facilitation increases with environmental severity
Community species richness showed a unimodal response along the environmental gradient and was significantly higher at intermediate

| Functional patterns of facilitation
The mean functional distance (MFD SES ) between all species in a bin showed a significant Gaussian response along the environmental gradient (p-val < .001, R 2 = .40; Figure 3a). Facilitating species tended to be more similar among each other than expected by chance at the warm/ dry edge, but this functional distance became more random at the cold/ wet edge of the gradient (p-val < .05, R 2 = .27, Figure 3b). Facilitated species showed the same pattern (p-val < .05, R 2 = .11, Figure 3c), although less pronounced than among the facilitating species.

| Species-specific trends in facilitation intensity
When considering each facilitated species independently, 33 of 46 facilitated species (72%) showed significant trends in facilitation intensity received along the environmental gradient (19 positive and 14 negative trends; average R 2 = .51; Figure 4a). For each facilitated species, we also tested whether the functional distance to its facilitating species changed along the gradient. Thereby, 22 of the 46 facilitated species (48%) showed significant trends along the environmental gradient (12 positive and 10 negative trends) with an average R 2 = .37 ( Figure 4b). Finally, considering each facilitated species independently, we found that 26 of 46 facilitated species (57%) showed significant relationships between their facilitation intensity received and their functional distance to their facilitators (15 positive and 11 negative relationships) with an average R 2 = .25 (Figure 4c).

| DISCUSSION
By introducing and validating a new analytical protocol for assessing facilitative interactions based on species co-occurrence and co-abundance patterns in community data, we are able to identify complex trends of facilitative interactions in species-rich communities and along extended environmental gradients. First, in the case study of the Zermatt region, facilitation intensity was generally strongest at high elevation where species were exposed to cold (but not drought) stress, although these communities contained also fewer facilitating species than dryer/warmer communities. Second, the functional distance between facilitating and facilitated species changed along the stress gradient and seemed to influence the facilitation intensity but with no general trend across species. Below, we discuss our results and evaluate the strengths and weaknesses of this new method.

| General trends
In agreement with theoretical expectations, plant community richness and functional diversity were highest at intermediate elevation in the Zermatt region (Michalet et al., 2006). This indicates stronger environmental filtering and, thus, stronger abiotic stress in communities at both the warm/dry and the cold/wet edge of our steep gradient (ranging from an average moisture index of 8 mm to 136 mm of remaining, nonevaporated precipitation per month; Lavergne, Mouquet, Thuiller, & Ronce, 2010;Webb, Ackerly, McPeek, & Donoghue, 2002).
Yet, the proportion of facilitating species was lowest at the cold/wet edge of the gradient (high elevations). This may be explained by the fact that facilitation under dry-warm conditions-via desiccation protection through shading for instance-is frequent and has a limited cost for the species able to grow in very dry-warm sites (e.g., Barbier, Couteron, Lefever, Deblauwe, & Lejeune, 2008;Maestre, Bautista, & Cortina, 2003; but see Maestre, Callaway, Valladares, & Lortie, 2009), whereas facilitation under cold condition-by sharing sparse nutrients or forming strong shelters for instance-comes at a greater cost, and thus, only few facilitating species may be able to provide it (e.g., cushion plants; Butterfield et al., 2013; but see Maestre et al., 2009).
Although the proportion of facilitating species was sparse at high elevations, it showed highest impact on the abundance of facilitated species (i.e., highest facilitation intensity). Our results are thus in line with previous findings, which stated that the intensity of facilitation is higher in cold environments (Callaway et al., 2002). But our results also demonstrate that facilitation is frequent in other types of stressful conditions (e.g.,

| Facilitation intensity and functional distances are linked
To better understand the facilitation interactions along our stress gradients, we explored the functional relationship between facilitating and facilitated species. We found three major results. First, facilitating species significantly resembled each other at the warm/dry edge of our gradient, while they tended to be functionally different at

| Limitation of the methodology and perspectives
The new methodology proposed here is simple and has strengths and weaknesses. On the one hand, it allows for identifying broad patterns of facilitation without experimental manipulations of the system (i.e., avoiding the introduction of unnatural levels of abiotic stress to the system; Körner, 2003) and It should be noted that we only grouped communities along one environmental gradient (i.e., moving window approach along the 1st PCA axis only). In this specific study system, the PCA axis used is in fact highly correlated with many other environmental variables, such as temperature, evapotranspiration, and moisture level (see Figure 2a). This is because our study system shows a very strong and dominating elevation gradient within a very small region (ca. 160 km 2 ), which did not provide sufficient independence in moisture and temperature to study these gradients separately. However, in highly heterogeneous systems where the main environmental drivers are less or not correlated, it would certainly be necessary to instead group communities along two or more environmental axes. In such a case, an even larger database of community plots might be required to have sufficient material for statistical analyses available.
Another important point not yet analyzed in our framework is that, at the community scale, co-occurring species can be at the same time facilitating (e.g., by modifying the local conditions) and competing with each other (e.g., by consuming the local resources). However, for predicting community dynamics, for example, under global change, quantifying the relative importance of competition and facilitation within communities is of utmost importance (McIntire & Fajardo, 2014). An analogous framework to the one presented here could be employed to Overall, our method provides information useful for the refinement of coexistence theory from a functional perspective. Indeed, we have shown that along the studied environmental gradient some species tend to be facilitated by functionally dissimilar species and others by functionally similar species. This is in contrast to prevailing predictions in community ecology that functionally dissimilar species rather co-occur due to competitive interactions, while similar species are expected to co-occur due to environmental filtering (Weiher & Keddy, 1999). Therefore, our results call for caution when using only the functional distance between species as an indicator of the underlying coexistence mechanisms, as facilitation processes alone may favor the co-occurrence of either similar or dissimilar species.
To conclude, we have introduced a simple and tractable method to identify and quantify facilitative interactions. Applying this method over a long moisture/temperature gradient in a species-rich system revealed new evidence that facilitation, similarly to competition, can operate between functionally similar and dissimilar species, and that these differences can change along environmental gradients. Applying this approach to other systems (e.g., savanna, tropics, and forest) and biotic levels (e.g., birds, amphibians, and arthropods) will offer vast opportunities to identify the main stress gradients for different taxonomic groups and regions, and help better understand the facilitation mechanisms prevailing in different environments.