Intra-seasonal rainfall variability and herbivory affect the interaction outcome of two dryland plant species

. Increases in drought frequency in combination with overgrazing may result in degradation of (semi-) arid ecosystems. Facilitative interactions between plants are a key mechanism in preventing degradation, but it is poorly understood how they respond to increased stress by combined drought and herbivory. In this study, we used an ecohydrological model, to simulate the plant growth of two plant species interacting with each other under different rainfall and herbivory pressure scenarios. The functional traits of the two modeled plants were based on a prior ﬁ eld experiment in southeastern Spain, in which an unpalatable “ nurse ” species protected a palatable prot ´ eg ´ e species from herbivory. Moreover, the nurse species was more drought-resistant; that is, it had a lower wilting point, whereas the prot ´ eg ´ e species had a higher optimal growth rate. Firstly, we investigated the coexistence of the two plant species growing under a single limiting resource, focusing on the effect of intra-seasonal rainfall variability. We found that longer periods without rainfall within the wet season resulted in stable coexistence, whereas nearly constant rainfall led to competitive exclusion of the prot ´ eg ´ e by the nurse species. Secondly, we investigated how plant interactions varied along our studied gradients. Using the neighbor effect intensity and importance indices, we found that competitive effects increased with more constant rainfall. Moreover, higher herbivory rates resulted in increased facilitative effects of the nurse on the prot ´ eg ´ e species, but facilitative effects could only prevail over competitive effects under currently observed or higher intra-seasonal rainfall variability. This study highlights the relevance of intra-seasonal rainfall variability in explaining coexistence of species in dryland ecosystems and shows that increasing intra-seasonal rainfall variability or herbivory pressure can result in more facilitative effects from a nurse species. This information is crucial to obtain a better insight into the long-term coexistence of species, and the resulting stability of dryland ecosystems in response to future climate change. scenarios (from low to high intermittency, IM = 1, 8, or 11 d) and three herbivory scenarios (from low to high herbivory intensity, h = 0, 0.1, 0.3 per yr). When calculating the indices, biomass values below 1g/m 2 were considered as equal to zero to avoid numerical arti-facts. A value of 1 for NInt C or NImp C represents maximum facilitative effects; a value of − 1 represents maximum competitive effects. An increase in intermittency resulted in a decrease in competitive effects. A combination of high herbivory and currently observed or high intermittency resulted in facilitative effects. a herbivory gradient (with steps of 0.1 per yr) for three intermittency scenarios (IM = 1, 8, or 11 d) and three annual rainfall scenarios (100, 300, or 500 mm/yr). A value of 1 represents maximum facilitative effects; a value of − 1 represents maximum competitive effects. When calculating the indices, biomass values below 1 g/m 2 were considered as equal to zero to avoid numerical artifacts. An increase in herbivory resulted in higher nurse effect intensity or importance.

(semi-) arid ecosystems. Facilitative interactions between plants are a key mechanism in preventing degradation, but it is poorly understood how they respond to increased stress by combined drought and herbivory. In this study, we used an ecohydrological model, to simulate the plant growth of two plant species interacting with each other under different rainfall and herbivory pressure scenarios. The functional traits of the two modeled plants were based on a prior field experiment in southeastern Spain, in which an unpalatable "nurse" species protected a palatable protégé species from herbivory. Moreover, the nurse species was more drought-resistant; that is, it had a lower wilting point, whereas the protégé species had a higher optimal growth rate. Firstly, we investigated the coexistence of the two plant species growing under a single limiting resource, focusing on the effect of intra-seasonal rainfall variability. We found that longer periods without rainfall within the wet season resulted in stable coexistence, whereas nearly constant rainfall led to competitive exclusion of the protégé by the nurse species. Secondly, we investigated how plant interactions varied along our studied gradients. Using the neighbor effect intensity and importance indices, we found that competitive effects increased with more constant rainfall. Moreover, higher herbivory rates resulted in increased facilitative effects of the nurse on the protégé species, but facilitative effects could only prevail over competitive effects under currently observed or higher intra-seasonal rainfall variability. This study highlights the relevance of intra-seasonal rainfall variability in explaining coexistence of species in dryland ecosystems and shows that increasing intra-seasonal rainfall variability or herbivory pressure can result in more facilitative effects from a nurse species. This information is crucial to obtain a better insight into the long-term coexistence of species, and the resulting stability of dryland ecosystems in response to future climate change.

INTRODUCTION
For the semi-arid Mediterranean area, climate models project an increase in high-temperature extremes (high confidence), an increase in meteorological drought frequency (medium confidence), and an increase in heavy precipitation events (high confidence) (IPCC 2014). Predicted climatic changes may result in widespread degradation of ecosystems in the Mediterranean basin (Guiot and Cramer 2016). Increased drought stress in combination with overgrazing can result in rapid degradation of (semi-) arid ecosystems (Kéfi et al. 2007), and facilitation between plants is suggested to be an important process keeping an ecosystem in a vegetated state (Verwijmeren et al. 2013, Xu et al. 2015, Kéfi et al. 2016. Therefore, there is an urgent need to understand how increased drought stress and herbivory pressure differentially and simultaneously impact changes in species interactions, thus influencing their coexistence. In semi-arid ecosystems, shrubs, trees, and other perennial or annual plant species compete for water, which is considered the main growthlimiting resource. Shrubs, however, often also provide positive (i.e., facilitative) effects by relieving drought stress for less drought-tolerant neighboring plants, for example by shading or by increasing water availability (via increased infiltration or hydraulic lift) within their direct vicinity (Pugnaire et al. 2011). Next to that, shrubs can lower consumer pressure (i.e., herbivory) by protecting neighboring individuals against herbivore damage, a process known as associational resistance (Hay 1986). The total net effect of one plant (nurse plant) on the other plant (protégé plant) is a trade-off between competitive and facilitative effects, and a crucial question to be answered is how a combination of different stressors (e.g., drought and grazing) influences the direction and strength of plant-plant interactions in semi-arid ecosystems (Soliveres et al. 2015).
Early conceptual models of plant-plant interactions hypothesized that the net outcome of plant-plant interactions shifts from competition toward facilitation with increasing drought stress or grazing pressure (Bertness and Callaway 1994, Callaway 1995, Callaway and Walker 1997. Meta-analysis indeed showed that at the global scale a shift toward more facilitative plant-plant interactions is observed as stress levels increase (He et al. 2013, Cavieres et al. 2014. However, recent studies question whether positive species interactions can be expected under very severe drought stress (Butterfield et al. 2016, Metz andTielbörger 2016), as competitive effects may become more intense during severe dry periods. In addition, studies from grazed ecosystems show that under severe grazing stress, plant-plant interaction wane from facilitation to neutral, as plants that provide benefits for neighbors lose their ability to do so under very high consumer pressure, for example, due to increased biomass loss and associated lower shading capacity but equal competition for water (Smit et al. 2007, Graff andAguiar 2011). So far, only very few studies empirically tested the impact of combined effects of drought and grazing (but see Maalouf et al. 2012, Verwijmeren et al. 2014, and changes in plant-plant interactions along combined stress gradients are still not well understood (Grant et al. 2014, Soliveres et al. 2015. Up until now, there are no two-species mechanistic model studies that investigated the simultaneous impact of these two stressors on changes in plant-plant interactions. A mechanistic modeling approach can be very insightful, as experimental or observational studies are often not able to assess plant-plant interactions along wide controlled gradients consisting of multiple stressors, or for extreme stress situations that are expected in the near future. Previous ecohydrological models investigated the balance between plant interactions with positive and negative effects (e.g., Gilad et al. 2007, Díaz-Sierra et al. 2010, Synodinos et al. 2015, by modeling the competitive water uptake by plants in combination with an increased infiltration of water in the soil due to increased biomass. These studies assessed which factors tilt the balance between the facilitative and competitive effects along a drought gradient. Positive interactions between plants have also been studied in a mechanistic modeling study (Gross 2008), in which grazing removal was made dependent on the biomass of a neighboring species, thus representing grazing protection. This mechanism can result in coexistence of multiple species that are all limited by the same single resource. However, it has not yet been studied with mechanistic v www.esajournals.org models how the joint effects of grazing and drought-highly realistic in (semi-) arid ecosystems-shape the net interactions between plants. Moreover, previous models used constant rainfall, ignoring the important aspect that rainfall in semi-arid ecosystems is highly intermittent; that is, it occurs in pulses (Chesson et al. 2004, but see Mathias and Chesson, 2013).
In the Mediterranean region, rainfall shows a high temporal variability; that is, it is characterized by pulse events. Several ecohydrological model studies addressed the role of stochastic and intermittent rainfall on the coupled dynamics of vegetation and soil moisture in dryland ecosystems (Rodriguez-Iturbe 2000, Baudena and Provenzale 2008, Kletter et al. 2009, Baudena et al. 2013, Siteur et al. 2014, D'Onofrio et al. 2015. Most of these studies, however, considered only one type of vegetation, and so far, only a few studies have considered the role of intra-seasonal rainfall variability (i.e., rainfall intermittency) on species coexistence (D'Onofrio et al. 2015), while, to our knowledge, none focused specifically on species interactions. Conceptually, it has been hypothesized that the temporal stochasticity in rainfall pulses is an important driver that can explain the coexistence of plants competing for a single resource by the so-called storage effect or by relative non-linearity in growth response to a limiting resource (Chesson 2000, Chesson et al. 2004, Adler et al. 2013, Barabás et al. 2018. Under the storage effect, species may have equal responses to a single limiting resource, but each species is favored at a different period of time. Relative non-linearity implies that plants have differences in growth rate responses to a single limiting fluctuating resource. This process allows one species to benefit from temporarily high resource availability, while another plant species is benefitting in times of temporarily lower resource availability, thus favoring long-term coexistence. In this study, we used a mechanistic two-species ecohydrological model to investigate how plant coexistence and plant-plant interactions are dependent on different rainfall amounts, rainfall intermittency, and herbivory pressure scenarios. Model assumptions and parameter settings are based on the experimental outcomes of a plant-plant interactions setup at a field site in southeastern Spain (Murcia region 37°57 0 28.37″ N, 1°0 0 16.14″ W), where a relatively drought-resistant nurse plant protected a protégé plant-that had a higher optimal growth rate-from being grazed (Verwijmeren et al. 2019). We modeled realistic (stochastic) rainfall scenarios to explore the effect of rainfall intermittency on the coexistence of the two species, which differed in their optimal growth rate and wilting points, and on the net effect of their interactions. Also, we explored the shifts in nurse effect intensity and importance along gradients of combined mean annual rainfall and herbivory under the different intermittency scenarios.

Model description
In order to study the interacting effects of annual rainfall amount, rainfall intermittency, and herbivory on plant coexistence and on the competitive or facilitative effects of a nurse species on a protégé species, we introduced a mechanistic twospecies model coupled to a hydrological model of a single soil layer. The model adopted here describes the coupled dynamics of vegetation and soil moisture, and it is a combination and extension of the models presented by Baudena et al. (2007), Laio et al. (2001), and Díaz-Sierra et al. (2010). We used a one-layer bucket model, as in our field site we observed relatively shallow soil depths (20-30 cm) accessible for both plant types root systems (M. Verwijmeren, personal observation). Water input in the model consisted of stochastic rainfall events based on statistics of historical data for the yearly amount and timing of rainfall. Parameter settings for rainfall and plant growth were chosen from field data as described in more detail below. The system dynamics was modeled by using three coupled ordinary differential equations (ODEs), for the soil water (s) dynamics, and for the nurse (N) and the protégé (P) plant growth dynamics.
Soil water dynamics.-The soil water (s) dynamics was a function of the infiltration, I, transpiration by the nurse species, T n , transpiration by the protégé species, T p , evaporation, E, and leakage, L (Eq. 1, Table 1). ds dt ¼ 1 nz ½IðN, P, r, tÞ À EðsÞ À T n ðN, sÞ À T p ðP, sÞ À LðsÞ where r is the rainfall, n is the soil porosity, and z is the soil depth.
Infiltration, I, was modeled as a function of the rainfall amount over time, and the protégé and nurse biomasses: The parameters i p and i n determined the positive effect of the protégé and nurse species, respectively, on the water infiltration in the soil layer. They represent the inverse of the half-saturation constant, implying that at our chosen values of i p = i n = 0.01, the biomass-dependent increase in infiltration is half of the maximum increase in biomass-dependent infiltration at a value of 100 g/m 2 . For s 0, we used a value of 0.5, implying that half of the rainwater infiltrates regardless of the positive effect of biomass on increased infiltration (Díaz-Sierra et al. 2010). If the soil layer was saturated, and the infiltration from Eq. 2 exceeded the available water storage in the soil, the excess was converted into surface runoff (Baudena et al. 2007).
Rainfall, r(t), was a stochastic function over time, in which a rainfall value was generated per day given a certain amount of annual rainfall and the inter-arrival time between rainfall events. We modeled rainfall as stochastic Poisson events, with exponential distributions for inter-arrival time (i.e., time in between rainfall events), and for mean daily rainfall intensity (calculated from the mean annual rainfall) (Laio et al. 2001). We also included a dry season without any rain, occurring once every year, to simulate the summer dry season that is characteristic for the climate of the field site. We calculated realistic values for mean annual rainfall, rainfall inter-arrival time, and length of drought season, based on 72 yr of rainfall records for the Alcantarilla weather station nearby our field site (Agencia Estatal de Meteorología, AEMET). The annual average rainfall was 300 mm per year. The average duration of the dry season was 61 d. The average time interval between rainfall events outside the dry season was 7.8 d (hereafter referred to as the current rainfall intermittency, IM = 8 d). We varied rainfall values and intermittency values in subsequent model runs, with a higher intermittency resulting in higher rainfall inter-arrival time and thus also increased mean daily rainfall intensity (see Appendix S1: Fig. S1 for an overview of the stochastic rainfall, soil water, and plant growth dynamics).
Evaporation, E, was modeled as a function of soil moisture content: (3)  Evaporation from the soil layer was zero if the volumetric soil water content dropped below the hygroscopic point s h of the soil (Baudena et al. 2012). When soil water content was above the hygroscopic point, it increased linearly as a function of the soil water content, up to a maximum value at saturation, E o (Kim et al. 1996, Baudena andProvenzale 2008).
Plant transpiration (T n , T p ) was modeled as a function of the soil moisture content and includes plant-specific wilting point (s wi ) optimum uptake point (s *i ), and maximum water uptake rate C i , with i = p or n for the protégé or nurse plant, respectively (as in Rodriguez-Iturbe 2000, Laio et al. 2001, Baudena et al. 2012. When the water content was below the wilting point, transpiration was assumed to be zero, as plants fully close their stomata below wilting point conditions. When soil moisture was above the wilting point, transpiration increased linearly as a function of the soil water content as plants open their stomata. When the soil water content was above s *i , transpiration reached an optimum (plants are assumed to have fully open stomata) at a constant rate C n N or C p P.
Leakage losses from the soil layer were modeled using a power law: where k s is the saturated hydraulic conductivity (Kim et al. 1996, Baudena et al. 2012). The exponent is related to the pore size distribution index and we used a factor 4, as recommended for loamy soils (Kim et al. 1996). Plant dynamics.-We modeled two plant types. Based on previous models Van De Koppel 1997, Díaz-Sierra et al. 2010), plant growth was modeled as proportional to transpiration, a function of the soil water content, with a proportionality constant that would determine the maximum growth rate, g maxi . The "nurse" species (N) had a baseline mortality rate, m n , but did not suffer from herbivory biomass removal: The second plant type, the "protégé" species (P), also had a baseline mortality, m p , on top of which a grazing mortality was implemented (third term on the right-hand side of Eq. 8 below). The biomass removal per year was proportional to its own biomass, with herbivory rate h. Herbivory damage was reduced depending on the ratio of nurse biomass over protégé biomass: The nurse species decreased the amount of herbivory-induced mortality for the protégé species. We followed the approach by Gross (2008), but while that model uses the herbivory protection as a function of the neighboring species alone, we choose to model herbivory protection as a function of the ratio of nurse biomass over protégé biomass, to account for size dependence (a small nurse plant cannot protect a larger protégé plant). Herbivory protection was modeled using a Holling type III function, where the ratio of the nurse biomass over the protégé biomass determined the amount of reduction in herbivory biomass removal. The parameter a determines how effectively the nurse can protect the protégé. More specifically, the square root of a is the half-saturation constant for herbivory protection along an axis of the ratio of nurse biomass over protégé biomass. For a, we choose a value of 1, so that the protective effect is 0.5 when both plants have a similar size (i.e., the ratio of N over P is one), as we assumed that half of the protégé plant in that case is protected by the nurse. When the ratio of N over P becomes higher than 4, the Holling type III function will approach 1, assuming that a nurse of greater than five times the biomass of the protégé will provide almost full protection, and the herbivory will not affect the protégé biomass.

Selection of parameters for the different plant types
In our model study, we distinguished two distinct plant functional types, a nurse (N) and a protégé (P) species. Both are woody perennial species that were selected in a parallel experiment in which plant growth of protected and unprotected planted saplings of Anthyllis cytisoides was monitored (Verwijmeren et al. 2019). In line with this experiment, we used Artemisia herba-alba as nurse plant in this study and modeled its growth. Artemisia spp. is not preferred by goats and has been found to be spatially associated with A. cytisoides in previous studies (Haase et al. 1996, Verwijmeren et al. 2014. Also, in line with field observations and parallel experiments (Verwijmeren et al. 2019), we used Anthyllis cytisoides as protégé species in our model. A. cytisoides is a drought-deciduous shrub from the Fabacaea family, and it is highly palatable for both goats and rabbits. A. cytisoides has been found to constitute 41% of livestock goat diet and is thus considered as highly preferred food source for goats (Barroso et al. 1995). As an optimal growth rate, C p g maxp , for A. cytisoides, we used five times the average yearly growth rate as measured in our experimental setup (M. Verwijmeren, personal observation), and for A. herba-alba, we used a lower optimal growth rate (M. Verwijmeren, personal observation). We used a difference of 0.01% between the wilting points of the two species in our model, namely a wilting point of 0.06 (as percentage of volumetric water content) for A. herba-Alba and a wilting point of 0.07 for A. cytisoides. Moreover, the effect of varying parameter settings for the two-species wilting points and growth rates s explored in Appendix S1 (Figs. S2-S4). In the default parameter setting, the two species are characterized by a trade-off: Although A. herba-alba has the benefit of being able to grow under lower soil moisture levels because of its lower wilting point, A. cytisoides has a higher growth rate under more benign moist conditions. A similar trade-off between drought tolerance and optimal growth rate has been reported in several studies in dryland ecosystems and has been proposed as a possible mechanism promoting plant coexistence (Chesson et al. 2004, Angert et al. 2009). As A. cytisoides is drought deciduous, and the study site was not grazed during the dry season, the herbivory rate on the protégé was set to zero in cases of a dry season that lasted for 30 or more days. Average life span was set to 30 yr for the protégé species and 25 yr for the nurse species, as A. cytisoides has been reported to have a longer life span than A. herbaalba (Haase et al. 1997).

Calculating nurse effects on the protégé
To measure the nurse effect on the protégé along gradients, we calculated the nurse effect intensity and the nurse effect importance, formerly called interaction intensity and importance (Armas et al. 2004, Brooker andKikividze 2008). The nurse effect intensity calculates the absolute effect of the nurse on the protégé biomass, whereas the nurse effect importance calculates the relative effect of the nurse on the protégé biomass, with respect to the effect that stress has on the protégé's biomass. As recent studies showed that the widely used indices RII (Armas et al. 2004) and I imp (Seifan et al. 2010)  Commutative intensity index.-For intensity, we used the Neighbor Effect Intensity index with commutative (and multiplicative) symmetry, NInt C (Eq. 9). This index weights the difference between the protégé's biomass with and without a nurse, with respect to the sum of biomasses of the protégé with and without a nurse. Being an intensity index, it standardizes the effect that the nurse exerts on the protégé's biomass at any point of the stress gradient, with respect to the sum of the protégé with and without a nurse at that same point of the stress gradient. This index has equal boundaries ranging from −1, indicating maximal competitive effects, to +1, indicating maximal facilitative effects. The index shows commutative symmetry, meaning that a twofold increase or decrease in the protégé biomass due to the nurse presence will result in equal deviations from 0 in the index value. This index is recommended for cases where facilitative effects and competitive effects are not in the same order of magnitude (Díaz-Sierra et al. 2017). We also calculated Neighbor Effect Intensity index with additive symmetry, NInt A , as provided in Appendix S1. Both intensity indices resulted in very similar patterns across our investigated gradients (Appendix S1: Figs. S5-S8).
where P sum = P −N + P +N is the sum of the performances of the protégé species without, P −N , and with a neighboring nurse, P +N , that is used to standardize for the protégé size, and where ΔP = P +N − P −N is the absolute impact of the nurse species on the protégé species.
Commutative importance index.-For importance, we used the Neighbor Effect Importance index with commutative symmetry: NImp C (Eq. 10). This index weights the difference between the protégé's biomass with and without a nurse, both with respect to the size of the sum of the protégé with and without a nurse species, and with respect to the effect that stress has on the protégé's biomass over the entire stress gradient. NImp C thus reflects the ratio between the neighbor-driven change in performance, standardized for size, and weighted for the change in performance driven by "all the factors in the environment that influence plant success" following the definition of interaction importance by Brooker and Kikividze (2008), but corrected to include standardization for size. This index has the same boundaries and symmetry as the above-described interaction index NInt C (Eq. 9).
where P sum and ΔP are the same as defined above, and MP sum is defined as the maximum value of the sum of the performances of the protégé species with and without neighbors at any point along the combined gradient. We also calculated the Neighbor Effect Importance index with additive symmetry, NImp A , as provided in Appendix S1. Both importance indices resulted in very similar patterns across our investigated gradients (Appendix S1: Figs. S5, S6).

Scenarios
To assess changes in nurse effect intensity and nurse effect importance along our herbivory, annual rainfall, and intermittency scenarios, we performed different model runs of 1000 yr, varying herbivory rate, the annual rainfall, and intermittency, respectively. We varied annual rainfall over a gradient ranging from 50 mm/yr up to rainfall of 500 mm/yr with steps of 50 mm/yr. Moreover, we varied the herbivory rate ranging from no herbivory (h = 0 per yr) to high herbivory (h = 0.3 per yr) with steps of 0.1 per yr. We varied the intermittency from 1 d, as a nearly constant control treatment, to 8 d, as based on current observed climatological values, to 11 d, as a high intermittency treatment, to test a possible increase in the intermittency as future climate change scenarios predict.
To simulate our experimental setup, that consisted of planted protégé sapling either sole standing or growing under the canopy of a A. herba-alba shrub, we ran all scenario's one time with the nurse species and the protégé species coexisting (both at initial values of 100 g/m 2 ), and one time with the nurse biomass set to zero (keeping the initial protégé biomass at 100 g/m 2 ). Our model reached a stochastic equilibrium after 100-200 yr, depending mainly on intermittency, among other parameter settings. To avoid including transient dynamics, we ran the simulations for 1000 yr and averaged the biomass density values over the final 200 yr, and used these for our comparisons. We also varied initial biomass of both the nurse and protégé species to start with either 1 or 1000 g/m 2 to check for possible multistability in the nurse or protégé biomass. We defined coexistence as a situation in which both species had non-negligible biomass (chosen as 10 g/m 2 ), and we qualitatively assessed coexistence based on figures showing biomass over gradients of rainfall or herbivory.

RESULTS
Nurse species average biomass values ranged from 0 till about 3400 g/m 2 at 500 mm of yearly rainfall. Protégé species average biomass values ranged from 0 till about 3200 g/m 2 at 500 mm of yearly rainfall. Biomass of nurse and protégé species and the competitive outcome were heavily dependent on the rainfall intermittency, annual rainfall amount, and herbivory rate variations (Figs 1, 2). We did not find any effects of varying the initial conditions of nurse or protégé biomass on the final biomass values, implying the modeled system did not display any multistability.

Effect of intermittency on competitive outcome
In the low rainfall intermittency scenario (IM = 1 d) and without herbivory (Fig. 1a), the v www.esajournals.org nurse species (black solid line) completely excluded the protégé species (black dashed line) along the whole rainfall gradient. With currently observed rainfall intermittency (IM = 8 d, darkgray lines), however, the protégé species was only completely excluded by the nurse at very low rainfall levels (<ca. 100 mm/yr). For high rainfall levels, we found stable coexistence between the nurse and the protégé species, with the protégé becoming more dominant at the higher end of the rainfall gradient, because of its higher growth rate under optimal conditions. With a further increase in rainfall intermittency (IM = 11 d, light-gray lines), the same pattern arose, but with slightly higher protégé biomass at low values of rainfall and a slightly lower protégé biomass at high rainfall levels, compared to protégé biomass at the currently observed intermittency. With an increase in herbivory to 0.1 per yr (Fig. 1b), the protégé could not become dominant at any point along the rainfall gradient, as the higher growth rate under optimal conditions of the protégé was counteracted by the biomass removal due to herbivory. With low intermittency, the protégé species was completely excluded. With an increase in intermittency, the protégé could persist, but at relatively low biomass values.
We varied the wilting point and optimal growth rate of the protégé species to investigate the effect on coexistence under different intermittency scenarios (see Appendix S1: Figs. S2-S4). This showed that coexistence could only occur when there was a trade-off between a lower wilting point (more drought resistance) and a lower optimal growth rate. We found two limiting cases for this where the nurse had an advantage but the protégé could still persist. When the optimal growth rates for the nurse and protégé were equal but the wilting point of the nurse was lower (Appendix S1: Fig. S3c), or with equal wilting points but higher growth rate for the nurse (Appendix S1: Figs. S4b), the protégé could still coexist with the nurse at high rainfall and intermittency conditions, because the mortality rate of the nurse was higher than the mortality rate of the protégé (m n > m p ).

Nurse effect intensity and importance along a rainfall gradient
In the low intermittency scenario (IM = 1 d; Fig. 2a, d, g), the protégé species without a neighboring nurse (P − , dashed blue line) performed better than the protégé species with a neighboring nurse (P + , dashed red line) along the entire rainfall gradient. This indicates that with low intermittency, competition for water was the dominant driver of the interaction outcome, even when herbivory pressure was intermediate or high. For currently observed rainfall intermittency and without herbivory (IM = 8 d, h = 0; Fig. 2b), the interaction remained competitive but the nurse was only dominant over the protégé species at low rainfall levels and the protégé could outperform the nurse at higher rainfall. With an increase in herbivory (IM = 8 d, and h = 0.1, 0.3 per yr; Fig. 2e, h), the protégé species with a neighbor could perform equally well or better than the protégé without a neighbor over the entire rainfall gradient. With high intermittency (IM = 11 d; Fig. 2c, f, i), we found similar results, but with a slightly better performance of the protégé with a nurse as herbivory increased. With high herbivory, the protégé without a nurse could not persist over the entire rainfall gradient (Fig. 2i).
Overall, an increase in intermittency resulted in a decrease in competitive effects of the nurse on the protégé. A combination of high herbivory and currently observed or high intermittency resulted in facilitative effects (Fig. 3). In the low intermittency scenario (Fig. 3, left panel) and without herbivory, both nurse effect intensity, NInt C , and nurse effect importance, NImp C , were close to zero or negative along the whole rainfall gradient, indicating competitive effects to be prevailing. Moreover, we found an increase in competition importance with increasing rainfall. Overall, with currently observed intermittency (Fig. 3, middle panel) nurse effect intensity and nurse effect importance increased compared to the low intermittency case, indicating competition to become less dominant. With no herbivory (solid light green line), the nurse effect intensity remained negative along the whole gradient, but competition intensity decreased with an increase in rainfall. With low herbivory intensity (dashed green line), effect intensity was negative at high rainfall but shifted to positive at intermediate and low rainfall levels. With high herbivory (solid dark green line), we found values of nurse effect intensity equal to 1 over the whole gradient, indicating obligate facilitation; that is, the protégé could only survive with a nurse. Nurse effect importance was also positive but close to zero, indicating facilitation to be prevalent but resulting in relatively small increases in biomass compared to optimal biomass values.

Nurse effect intensity and importance along a herbivory gradient
Overall, an increase in herbivory resulted in higher intensity or importance of facilitative effects. In the low intermittency scenario (Fig. 4a,  d, g), the nurse species (black solid line) was constant over the herbivory gradient, as the protégé with a nurse (dashed red line) was completely excluded along the entire gradient. The protégé without a nurse (blue dashed line) could persist when being grazed, but only for high rainfall values or low herbivory rates. For the currently observed intermittency scenario (Fig. 4b, e, h), the protégé with a nurse could persist at higher rainfall values and became dominant over the nurse when there was no herbivory. For the high intermittency scenario (Fig. 4c, f, i), biomass Fig. 3. The effect intensity NInt C and the effect importance NImp C of the nurse on the protégé biomass along a rainfall gradient (with steps of 50 mm/yr) for three intermittency scenarios (from low to high intermittency, IM = 1, 8, or 11 d) and three herbivory scenarios (from low to high herbivory intensity, h = 0, 0.1, 0.3 per yr). When calculating the indices, biomass values below 1g/m 2 were considered as equal to zero to avoid numerical artifacts. A value of 1 for NInt C or NImp C represents maximum facilitative effects; a value of −1 represents maximum competitive effects. An increase in intermittency resulted in a decrease in competitive effects. A combination of high herbivory and currently observed or high intermittency resulted in facilitative effects.
values showed a similar pattern, but with a smaller difference between biomass values for P + and P − under low herbivory and high rainfall combinations.
Overall, increased herbivory resulted in a shift from negative to positive (or zero) values both for the nurse effect intensity and the nurse effect importance (Fig. 5), indicating a shift from competitive to facilitative effects being dominant. In the low intermittency scenario (Fig. 5, left panels), the intensity shifted from −1 to 0 with increased herbivory, and importance changes from competitive to facilitative but with values close to zero. In the currently observed intermittency scenario (Fig. 5, middle panel), with low herbivory (h up to 0.1 per yr) we observed negatives values for the nurse effect intensity, whereas with an increase in herbivory, the intensity shifted to positive and gradually increased toward 1, indicating that the protégé could only survive when it was protected by a nurse species. The trends in importance were roughly equal for currently observed and high intermittency scenarios (Fig. 5, right panels). The nurse effect importance shifted from negative to positive with increased herbivory, but stayed at very low positive values for high herbivory scenarios.

DISCUSSION
Our mechanistic two-species ecohydrological model results showed that intra-seasonal rainfall variations need to be taken into account as they influence also the competition and facilitation intensity and importance, which was not yet considered in conceptual models on species interactions along stress gradients (Maestre et al. 2009, Soliveres et al. 2015. Also, we showed that coexistence between plant species growing under a single limiting resource can arise under stochastic rainfall due to differences in growth rates as a function of the temporal variability in soil moisture. The two species could coexist because of their opposite water strategies, with a trade-off between growth rate for scarce and optimal water availability, which in our model corresponded to the wilting point being lower for the one plant species that also had a lower optimal growth rate under high resource availability then the other plant species. In addition to trade-offs in favorable conditions, relative non-linearity in growth rates could contribute to coexistence with water as the sole limiting resource under stochastic, realistic rainfall. Due to its lower wilting point, the nurse species will take up water during the moments in time when soil water content exceeds its own wilting point. This way, with low intermittency or (nearly) constant rainfall, the nurse species prevents soil water content to exceed the wilting point of the protégé species, thereby preventing the protégé species from increasing in biomass. For this reason, under more constant rainfall scenarios the nurse is dominant, as it has strong competitive effects on the protégé species. For higher intermittency or higher annual rainfall, stable coexistence between the nurse and the protégé species is possible. Under more intermittent rainfall, rainfall is accumulated in fewer events, resulting in increased mean daily rainfall intensity. Because of this, the soil water will reach more often values above the wilting point of the protégé, enabling the protégé to coexist next to the nurse species.
Our results show that considering rainfall intermittency is essential to make mechanistically sound predictions on how plant-plant interactions vary along drought gradients, as not only the total annual rainfall amount affects the competitive outcome, but also the temporal distribution in rainfall events. As the nurse species was a Fig. 5. The effect intensity NInt C and the effect importance NImp C of the nurse on the protégé biomass along a herbivory gradient (with steps of 0.1 per yr) for three intermittency scenarios (IM = 1, 8, or 11 d) and three annual rainfall scenarios (100, 300, or 500 mm/yr). A value of 1 represents maximum facilitative effects; a value of −1 represents maximum competitive effects. When calculating the indices, biomass values below 1 g/m 2 were considered as equal to zero to avoid numerical artifacts. An increase in herbivory resulted in higher nurse effect intensity or importance.
better competitor under low resource availability due to its lower wilting point, the competitive effect was most pronounced under low intermittency in rainfall. The protégé could only effectively survive under higher intermittency scenarios because in those cases the soil water level would more often be higher than the wilting point of the protégé, allowing that species to grow. We note that this result of our model study is due to the specific trait combination of the nurse species having a lower wilting point and the protégé having a higher optimal growth rate, but such trade-offs between drought tolerance and growth capacity are seemingly common (Reich 2014) and can even occur within plant species that are considered of the same functional type (Angert et al. 2009). This underlines the importance of taking into account explicitly the plant traits related to the mechanisms that can result into competitive or facilitative effects, when making predictions on plant-plant interactions along stress gradients.
The protégé is not benefited by nurse biomass, and competition for resources is the driving mechanism when there is no herbivory, as both species in our model increased the infiltration rate equally. The protégé is only being facilitated effectively when it is being grazed, and with increased herbivory, it receives more facilitative effects from the nurse. Under very severe herbivore pressure, plant-plant interactions have been observed to wane (Smit et al. 2007, Smit et al. 2009, Graff and Aguiar 2011. The loss of net facilitation under extreme levels of consumer pressure can occur either if the nurse species reaches a limit in the ability to reduce herbivory damage, or because the positive effects from herbivory protection do not compensate the costs of sharing resources between neighboring plants (Michalet et al. 2006, Verwijmeren et al. 2013). Our results did not show a clear waning of interactions as herbivory in our model only affected the protégé species and not the nurse species. Also, the protégé received almost full protection against herbivory damage if the nurse biomass was more than five times the protégé biomass. Further exploration of the ability of nurse shrubs to give either full or partial protection against herbivory related to their size and traits is needed to make modeling studies on plant-plant interactions more realistic.
Previous studies showed that interaction intensity and importance do not need to be correlated along a stress gradient (Brooker et al. 2005, Kunstler et al. 2011. Our results show that very small changes in biomass levels can result in drastic shifts in the sign of the nurse effect intensity that are not observed in the nurse effect importance. Any intensity index by definition will amplify tiny differences in biomass if the protégé with or without nurse tends to zero, which gives constraints in using this index for low values of biomass (Díaz-Sierra et al. 2017). The importance index is less sensitive than intensity to minor changes in the parameter values (not shown). Thus, generally importance seems to be a more appropriate index for this type of studies if the full gradient of a species is considered. Our results also showed that nurse effect intensity and nurse effect importance can show contrasting trends along the same stress gradient. In the intermediate intermittency case, the nurse effect intensity under ungrazed conditions increased with rainfall annual amounts, whereas the nurse effect importance decreased. This difference in trend can be explained by the different standardization for intensity and importance indexes. As the intensity index is standardized for size, the standardized difference between P +N and P −N decreases for larger plants at the more productive, wet part of the gradient; that is, the impact of standardizing for size increases with increased rainfall, moving the index value closer to zero. As the importance index is standardized for size, but is also weighted for the effect of stress over the total gradient, this index moves toward zero with increasing stress because the standardized difference between P +N and P −N decreases with respect to the effect of stress when moving to the more stressed dry end of the gradient, that is, the impact of standardizing for stress makes the importance index tend to zero with increased stress.
Along the herbivory gradient, the trends in importance were roughly equal for the three intermittency scenarios and Nimp C remained at very low positive values for high herbivory scenarios. This low importance is due to the very low biomass at high rates of herbivory, compared to the optimum biomass at no herbivory. This shows that when the most pronounced facilitative effects of the nurse on the protégé biomass v www.esajournals.org occur at the most stressed point of the gradient, the importance will show very low values. For this reason, the NImp C shows an opposite trend along the herbivory gradient than along the rainfall gradient, where the positive effects of the nurse are largest at the least stressed, high rainfall end of the gradient.

Synthesis
Our model study shows the relevance of rainfall intra-seasonal variability in determining the plant-plant interaction outcome between two species. Moreover, our results show that if a nurse species has a lower wilting point than the coexisting protégé species, competitive effects can be expected to increase with decreases in intermittency. Facilitative effects from the nurse did increase with increased herbivory pressure or increased intermittency. This work shows the value of using mechanistic two-species ecohydrological models to evaluate the joint effects of different stressors on plant-plant interactions, which can be useful for making future predictions of plant-plant interactions and stability of semi-arid ecosystems.