Conditional love? Co‐occurrence patterns of drought‐sensitive species in European grasslands are consistent with the stress‐gradient hypothesis

Abstract Aim The stress‐gradient hypothesis (SGH) postulates that species interactions shift from negative to positive with increasing abiotic stress. Interactions between species are increasingly being recognized as important drivers of species distributions, but it is still unclear whether stress‐induced changes in interactions affect continental‐to‐global scale species distributions. Here, we tested whether associations of vascular plant species in dry grasslands in Europe follow the SGH along a climatic water deficit (CWD) gradient across the continent. Location Dry grasslands in Europe. Time period Present. Major taxa studied Vascular plants. Methods We built a context‐dependent joint species distribution model (JSDM) to estimate the residual associations (i.e., associations that are not explained by the abiotic environment) of 161 plant species as a function of the CWD based on community data from 8,660 vegetation plots. We evaluated changes in residual associations between species for pairs and on the community level, and we compared responses for groups of species with different drought tolerances. Results We found contrasting shifts in associations for drought‐sensitive and drought‐tolerant species. For drought‐sensitive species, 21% of the pairwise associations became more positive with increasing CWD, whereas 17% became more negative. In contrast, only 17% of the pairwise associations involving drought‐tolerant species became more positive, whereas 27% became more negative in areas with a high CWD. Additionally, the incidence of positive associations increased with drought for drought‐sensitive species and decreased for drought‐tolerant species. Main conclusions We found that associations of drought‐sensitive plant species became more positive with drought, in line with the SGH. In contrast, associations of drought‐tolerant species became more negative. Additionally, changes in associations of single species pairs were highly variable. Our results indicate that stress‐modulated species associations might influence the distribution of species over large geographical extents, thus leading to unexpected responses under climate change through shifts in species associations.


| INTRODUC TI ON
). Many studies have provided evidence for the SGH along direct and indirect physical stress gradients, such as salinity and elevation (Bertness & Ewanchuck, 2002;Callaway et al., 2002), and for resource gradients, such as water availability (Dohn et al., 2013;He et al., 2013;Liancourt et al., 2005;López et al., 2016;Ziffer-Berger et al., 2014), and biotic stress gradients, such as herbivory (Graff & Aguiar, 2011). Of these, the gradient in water availability is particularly relevant to investigate given the expected changes in the magnitude and frequency of droughts resulting from global climate change (Vicente-Serrano et al., 2020).
Evidence for the SGH along water stress gradients is mixed, because several studies have found shifts from facilitative to competitive interactions under extreme drought stress (Berdugo et al., 2019;Maestre & Cortina, 2004;Soliveres et al., 2011;Soliveres & Maestre, 2014). This has led to refinements of the original formulation of the SGH, in which facilitative interactions are most prevalent in intermediately stressful conditions (Maestre et al., 2009;Michalet et al., 2014). Furthermore, few studies have investigated changes in interactions along large-scale water stress gradients, and these have typically focused on community-level metrics, such as biomass production (Dohn et al., 2013) and overall facilitation frequency or importance (Berdugo et al., 2019;Soliveres & Maestre, 2014), or cover only a few species (Ziffer-Berger et al., 2014). However, community-level responses to wide stress gradients can arise from the turnover of species with contrasting stress tolerances (Berdugo et al., 2019;Liancourt et al., 2017). Additionally, species with different life-history strategies can show opposite responses to their neighbours, which might obscure community-level patterns (Graff & Aguiar, 2017;Michalet et al., 2015). Hence, it remains unclear whether the large-scale distributions of plant species are influenced by drought-modulated species interactions and whether these follow expectations of the SGH.
Here, we tested whether shifts in plant species associations along the growing season climatic water deficit (CWD GS ) gradient are consistent with the predictions of the SGH. We focused on dry grasslands in Europe, which encompass a wide water deficit gradient and are rich in species, hence including many possible species associations (Wilson et al., 2012). We fitted a context-dependent joint species distribution model (JSDM; Tikhonov et al., 2017) to infer the pairwise associations of 161 species along the CWD GS . With this model, we identified residual associations between species, that is, those that cannot be explained directly by the abiotic environment, as a function of the CWD GS . We then investigated possible changes in associations along the CWD GS for pairs of species and at the community level. We assessed associations for drought-tolerant and drought-sensitive species separately because a water deficit acts as a stress mainly for drought-sensitive species, but not for droughttolerant species (Liancourt et al., 2005). We, therefore, hypothesized that with increasing water deficit, drought-sensitive species would be increasingly facilitated by other plant species, resulting in more positive associations. In contrast, we expected that drought-tolerant species would not be affected by the water deficit gradient, hence their associations would not change.
The vegetation plots, also called phytosociological relevés, include records of plant taxon co-occurrence at particular sites at a high spatial resolution (generally between 4 and 25 m 2 ), hence they can be considered suitable for analysis of co-occurrence patterns in local communities. We selected plots classified as dry grasslands under the European Nature Information System (EUNIS) habitat classification. From these, we selected plots recorded between 1979 and 2013 and with a location uncertainty ≤ 1 km in order to match with the spatial and temporal resolution of the environmental variable data ( Figure 1a; for a full overview of all databases included in this study, see Appendix 1; Supporting Information Table S1). This led to a dataset of 20,722 vegetation plots located in 8,660 unique 30 arc-s grid cells (for the distribution of plots per grid cell, see Supporting Information Figure S1).

| Environmental data
We selected 11 environmental variables that are known to influence plant distributions (Austin & Van Niel, 2011). To account for soil fertility, we included soil pH and soil cation exchange capacity (as centimoles of positive charge per kilogram of soil), which affect the potentially available nutrient supply (Chapin et al., 2011).
Additionally, we included soil organic carbon content (per-mille) as an indicator of soil texture and water infiltration rates and because it is an important nutrient reservoir (Chapin et al., 2011). We also included the following temperature variables (Prentice et al., 1992): maximum temperature of the warmest month (degrees Celsius), minimum temperature of the coldest month (degrees Celsius), growing degree days (GDD; > 5°C) and freezing degree days (FDD; < 0°C).
Given that water availability for plants is not determined by total precipitation alone (Piedallu et al., 2013), we included the actual evapotranspiration (AET) in the growing season (AET GS ; in millimetres) and the precipitation seasonality (coefficient of variation of monthly precipitation). We calculated monthly AET from monthly temperatures, monthly precipitation and available soil water capacity and summed monthly AET over the growing season (May-August) to obtain the AET GS (Lutz et al., 2010; for details, see Supporting Information Appendix S1). Lastly, we included the climatic water deficit during the growing season (CWD GS ; in millimetres) as an indicator of drought in our models (Figure 1b). We calculated the CWD GS as the difference between the potential evapotranspiration during the growing season (PET GS , in milimeters) and the AET GS (Lutz et al., 2010).
Soil variables are publicly available at a 30 arc-s resolution (~1 km) at http://soilg rids.org (Hengl et al., 2014). Precipitation seasonality, temperature seasonality, maximum and minimum temperature and monthly temperature maps are also available at a 30 arc-s resolution, as mean values from 1979 to 2013, from the CHELSA Climatologies dataset v.1.4 (http://chels a-clima te.org; Karger et al., 2017). We calculated the GDD, FDD, PET GS , AET GS and CWD GS from mean monthly minimum, maximum and mean daily temperature and mean monthly precipitation averaged across 1978-2013 from the same dataset (for calculations, see Supporting Information Appendix S2).

| Context-dependent JSDM
We modelled species occurrence using a context-dependent JSDM . JSDMs are able to separate co-occurrence patterns between species into shared environmental responses and residual associations that cannot be explained by the environment Pollock et al., 2014). In original JSDMs, the estimated residual associations are static over space. However, recently, Tikhonov et al. (2017) developed a new context-dependent latent variable approach, in which the residual associations are allowed to vary across the environment. We followed this approach and modelled species occurrences (y) as a probit regression: where z ij is the latent occurrence score for species j in vegetation plot i. This is calculated as : (1) y ij = 1 z ij >0 where n c is the number of environmental predictors plus intercept, β jk are the estimated regression coefficients, x ik are the measured environmental variables, n h is the number of latent variables, λ jh are the regression coefficients to the latent variables (latent loadings) and η ih is the latent factor site score. To calculate the residual associations between two species as a function of the CWD GS (R j1j2,CWD ) we first calculated the covariance matrix of the species' latent loadings as a function of the CWD GS (Ω j1j2,CWD ) as: We then transformed this covariance matrix into a correlation matrix to arrive at the residual associations (Equation 4).
We fitted the model using a Bayesian approach with an uninformative prior distribution, as described by Ovaskainen et al. (2017) and Tikhonov et al. (2017). We ran the model with three chains using a burn-in period of 1,500,000 iterations, after which we sampled 1,250,000 iterations, of which we saved 5,000 iterations to construct the posterior of the model. We evaluated model convergence based on the maximum potential scale reduction factor and minimum effective sample size (Supporting Information Figures S3 and S4; Brooks & Gelman, 1998;Gelman & Rubin, 1992). A more detailed description of the model, the priors and the fitting procedure can be found in the Supporting Information (Appendix S3).
We accounted for collinearity between the environmental variables by excluding environmental variables with Pearson's r > .7.
In total, six of the initial environmental variables were retained: soil pH, soil cation exchange capacity, soil organic carbon content, precipitation seasonality, minimum temperature of the coldest month and CWD GS (for correlations and selection rationale, see Supporting Information Figure S2). Given that plant species might show nonlinear responses to the environment, we also included the quadratic terms of each of the environmental variables. From each 30 arc-s grid cell, we randomly selected only one plot to be included in the final model in order to prevent pseudo-replication, leading to a total of 8,660 vegetation plots. Furthermore, we included only species with ≥ 10 presences per predictor in the selected vegetation plots (Peduzzi et al., 1996; leading to a minimum of 320 presences based on 12 environmental variables +20 latent variables). The final model included a total of 161 species, of which 35 were graminoids, 122 forbs and four shrubs (Supporting Information Table S2).

| Species groups
Given that we expected a difference in association shifts between drought-tolerant and drought-sensitive species, we divided the 161 species into three groups based on their Ellenberg values for soil moisture (E m ). We classified species with E m ≤ 3 as drought tolerant (DT; 36 species) and species with E m ≥ 5 as drought sensitive (DS; 43 species) (terBraak & Gremmen, 1987). For all subsequent analyses, we split the results such that they included only pairwise associations where at least one of the species belonged to each of these two groups. We included also within-group associations, because increasing abiotic stress can decrease competitive interactions when the competitive ability of stress-sensitive dominant plant species is reduced (Liancourt et al., 2005). We considered species with 3 < E m < 5 as intermediately tolerant and did not explicitly test the SGH for these species, because we did not have an a priori hypothesis regarding their responses. Note that the associations of intermediately tolerant species with drought-tolerant or drought-sensitive species were included in their respective analy- ses. Additionally, we tested a more stringent classification of species by considering only drought specialists (E m ≤ 2.5; 13 species) and highly drought-sensitive species (E m ≥ 5.5; 25 species). We retrieved the Ellenberg values for soil moisture from a list compiled by Louette et al. (2010). To test whether the CWD GS gradient is indeed a stress gradient for the drought-sensitive but not for the drought- (2) were considered to be indifferent to changes in the CWD GS . The choice to compare the associations at the 25th and 75th percentiles of the CWD GS value can be considered conservative, because this corresponds to a relatively short stress gradient (He et al., 2013; but see also Soliveres & Maestre, 2014, who found that longer stress gradients are more conservative that shorter ones). To test the sensitivity of our results to this choice, we also performed the same analysis using the 5th (246 mm) and 95th (819 mm) percentiles of the CWD GS gradient.

| Community associations
In addition to testing the shift in association strength for each pair of species, we analysed whether there is a higher overall proportion of positive compared with negative interactions with increasing CWD GS , summarized as the community association (CA): where N positive is the number of positive associations and N negative is the number of negative associations at a given CWD. For N positive and N negative , we included only associations for which the 95% credible interval did not overlap with zero. The value of CA CWD ranges from minus one, where all significant associations are negative, to plus one, where all significant associations are positive. We did not use the association strengths retrieved from the models directly, because these have been shown to be biased by species prevalence and cannot be used to compare interaction strengths among species pairs (Zurell et al., 2018

| Pairwise associations
In accordance with our expectations, the pairwise associations of drought-sensitive species more often increased (21%) than were indifferent to changes in CWD. When considering a wider CWD GS gradient (5th-95th instead of 25th-75th percentiles), we found a slightly higher proportion of drought-sensitive species with a positive response (23%), but otherwise the results were similar (Supporting Information Figure S5). A more stringent classification of species into highly drought sensitive and drought specialists yielded similar but more pronounced patterns with fewer negative shifts for highly drought-sensitive species (15%) and fewer positive and more negative shifts for drought specialists (12% and 32%, respectively; Supporting Information Figure S6).

| Community associations
We found a reduction in the number of drought-sensitive species and an increase in the number of drought-tolerant species with increasing CWD GS (Figure 4a). The CA for drought-sensitive species increased with the CWD GS in the second half of the gradient, that is, there were more positive compared with negative associations in dry than in wetter environments, which is in line with the expectations from the SGH (Figure 4b). This trend was mainly attributable to associations between drought-sensitive and other species, rather than associations within the drought-sensitive species group (Figure 4c). In contrast, the CA for drought-tolerant species decreased with the CWD GS , whereas we expected the CA for this group to be indifferent to drought stress (Figure 4b). This trend reflected associations with species of intermediate sensitiv- ity to drought in addition to associations within the species group itself (Figure 4c).
Associations between sensitive and tolerant species decreased with increasing drought along the first half of the gradient. Using the alternative Ellenberg indicator cut-off values to classify the species into drought specialists and highly drought-sensitive species revealed broadly similar patterns for sensitive species (Supporting Information Figure S7). The overall CA of drought specialists, in addition to the CA of drought specialists with intermediately tolerant species, followed a steeper decline over the gradient than those for drought-tolerant species. Additionally, while associations within the drought-tolerant group declined, associations within the drought specialist group did not change with the CWD GS (Supporting Information Figure S7). In general, we found more positive than negative associations along the entire climatic water deficit gradient, that is, the CA was always positive, except for associations between drought-sensitive and drought-tolerant species (Figure 4b Figure S8).

| D ISCUSS I ON
Traditionally, species distributions have been considered to be determined mostly by abiotic factors, whereas our results suggest that biotic interactions also play a role and that these interactions might change along environmental gradients. More specifically,

F I G U R E 3
Mean residual association of each species pair involving one drought-sensitive species (left panels) or one droughttolerant species (right panels) as a function of the growing season climatic water deficit (CWD GS ) for three different response types (see Figure 2). Top panels show associations that are more positive at the 75th percentile (665 mm) of the CWD GS than at the 25th percentile (420 mm). Middle panels show associations that become more negative with increasing CWD GS . Bottom panels show associations that do not change significantly along the CWD GS gradient. Percentages indicate the percentage of associations per response type per species group. The CWD GS gradient ranges from low (0 mm) to high water deficit (> 1,100 mm) we showed that the species-to-species associations estimated from the distribution of plant species across European dry grasslands seem to be modulated by the CWD GS . We found contrasting patterns among groups of species with different drought tolerances. Drought-sensitive species became more positively associated with intermediately drought-tolerant species along the gradient, but not with drought-tolerant species. In turn, droughttolerant species co-occurred less with other species with increasing drought. Furthermore, individual species pairs often deviated from the SGH, because the majority of species-to-species associations did not shift along the gradient, indicating the importance of species-specific factors other than stress tolerance. Our results conflict with the initial formulation of the SGH, which predicts a monotonic increase of positive interactions with increasing environmental severity (Bertness & Callaway, 1994). Instead, we found that changes in associations along the drought gradient depend on the stress tolerance of species, which is rather in line with refinements to the SHG (Graff & Aguiar, 2017;Maestre et al., 2009;Michalet et al., 2014).
We found a positive relationship between the predicted number of drought-tolerant species and the CWD GS , indicating that increasing drought does not represent a stress gradient to these species. Furthermore, the proportion of positive associations for drought-tolerant species, as indicated by the CA, decreased with increasing CWD GS . This might indicate an increased competition for water along the gradient, while these species do not benefit from facilitation. Similar results were found for stress-tolerant shrubs and grasses in an experimental field study in an arid steppe (Graff & Aguiar, 2017). However, our results also indicate that drought-tolerant species often do not act as nurses for drought- associations between drought-sensitive species themselves and between highly drought-sensitive and drought-specialist species did not vary along the gradient (Supporting Information Figure S7).
These findings indicate that facilitation of drought-sensitive species in dry environments is mainly governed by species with intermediate drought tolerance. Although neighbours might increase water availability in dry environments, they also compete for resources, and the outcome of an interaction is positive only when facilitative effects outweigh competition (Bimler et al., 2018;Michalet et al., 2014).  (Bimler et al., 2018). However, given that directionality cannot be estimated from co-occurrence, nor does it provide estimates of fitness, these hypotheses cannot be tested with our approach.
Our results showed a large variation in association shifts among species pairs. The majority of species associations did not change significantly along the CWD GS gradient. As such, it seems that associations of single species pairs do not necessarily comply with the SGH, although the frequency of positive associations within the community follows the SGH. This is in line with previous studies reporting that pairwise interactions are determined mainly by speciesspecific traits (He et al., 2013) or life-history stage . This has also been proposed as an important factor underlying the conflicting evidence for the SGH in the past (He et al., 2013).
An interesting problem that might be explored in future studies is to determine which species traits, besides water stress tolerance, can predict association shifts. The traits that mediate the facilitation between species vary depending on the mechanism of facili-  (Pescador et al., 2020), highlighting that there is a large variability in habits and traits of potential nurse plants.
Relating plant traits to interaction shifts is challenging, because shifts can arise not only from facilitation, but also from changes in the competitive ability of drought-intolerant dominant species (Liancourt et al., 2005). For example, rather than acting as a nurse, the graminoid Bromus erectus (49% positive shifts with droughtsensitive species) might have benefitted from reduced competition with increasing drought (Liancourt et al., 2005). Additionally, interactions also depend on the traits of the facilitated species and the similarity of traits of the interacting species, because nurses generally facilitate species with traits different from their own (He et al., 2013;Navarro-Cano et al., 2021). Lastly, this picture is complicated further by the plasticity of these plant traits to plant fitness and life stage. In very stressful environments, the nurse plants might perform too poorly to realize the microclimatic habitat amelioration (Schöb et al., 2013). A systematic analysis of the prevalence of different traits along the drought stress gradient was beyond the scope of our analysis and might require locally measured trait data, given the plasticity of traits involved in the amelioration of drought stress.
However, this would be an interesting avenue for further research.
Additionally, given that the outcome of interactions depends on both competitive and facilitative effects, traits explaining differences in competitive abilities along the gradient, such as the Ellenberg indicator values of species to soil moisture used in the present study, also need to be included (Bimler et al., 2018).
Species associations in JSDMs indicate that two species co-occur more or less frequently than expected from the abiotic environment.
We tried to overcome missing environmental covariates by including a broad set of environmental variables considered relevant for dry grassland plants. However, the resolution of these variables (30 arcs) is relatively coarse compared with the size of the vegetation plots and does not capture the local variability in topography and soil characteristics that affect soil moisture (Kemppinen et al., 2018).
Unfortunately, we did not have high-precision vegetation plot coordinates, and higher-resolution environmental data would require locally acquired measurements per vegetation plot, which are not available. Furthermore, interactions between species generally act at the scale of individuals, and their impacts might vanish with increasing plot size (Thuiller et al., 2015), especially for competition (Rajaniemi et al., 2006). Although the vegetation plots used in the present study were designed to capture the local plant communities and were generally small in size, they might, nevertheless, be too coarse for certain interactions to be discernible. All in all, inferences on species interactions from species associations must be made with caution (Blanchet et al., 2020).
In the present study, we showed that the large-scale distributions of plant species are mediated by shifts in the associations of species along a drought stress gradient. However, changes in species associations can also result from other stress gradients, such as nutrient availability, vapour pressure deficit and herbivory (Guignabert et al., 2020;He et al., 2013 does not allow us to determine the directionality of potential interactions, and the estimated associations between species cannot be used for predictions of future species distributions. We, therefore, conclude that stress-modulated species interactions, such as those predicted by the SGH, might have unforeseen consequences for future species distributions in response to climate change.

ACK N OWLED G M ENTS
The authors thank Helge Bruelheide and Ute Jandt for their feed-

S U PP O RTI N G I N FO R M ATI O N
Additional Supporting Information may be found online in the Supporting Information section.