The evolutionary stability of plant antagonistic facilitation across environmental gradients and its ecological consequences: soil resource engineering as a case study

Plant interactions, understood as the net effect of an individual on the fitness of a neighbor, vary in strength and can shift from negative to positive as the environmental conditions change in time and space. Evolutionary theory questions the stability of non-reciprocal interactions in which one plant has a positive net effect on a neighbor, which in return has a negative net impact on its benefactor. This type of interaction is known as antagonistic facilitation. We develop a spatially explicit consumer-resource model for below-ground plant competition, including plants able to mine resources and make them available for any other plant in the community, termed ecosystem engineers. We use the model to assess whether and under which environmental conditions antagonistic facilitation via soil resource engineering is evolutionarily stable. We find that antagonistic facilitation is stable in highly stressful conditions, which supports the theory of ecosystem engineers as drivers of primary succession and provides a theoretical ground to investigate facilitation mechanistically in the context of the stress gradient hypothesis. Among all potential causes of stress considered in the model, the key environmental parameter driving changes in the interaction between plants is the proportion of the limiting resource available to plants without mining. This finding represents a challenge for empirical studies, which usually measure the resource input or loss in the system as a proxy for stress. We also find that the total root biomass and its spatial allocation through the root system, often used to measure the nature of the interaction between plants, do not predict facilitation reliably. Synthesis. Antagonistic facilitation established between an ecosystem engineer nurse plant and neighbor opportunistic individuals can be evolutionarily stable in stressful environments where ecosystem engineers’ self-benefits from mining resources outweigh the competition with opportunistic neighbors. These results align with theories of primary succession and the stress gradient hypothesis as they show that antagonistic facilitation is stable under environmental stress, but it evolves into mutual interference in milder environments. However, using inaccurate parameters to measure facilitation and stress gradients in empirical studies might mask these patterns.

which in return has a negative net impact on its benefactor. This type of interaction is known as antagonistic facilitation.
2. We develop a spatially explicit consumer-resource model for below-ground plant competition, including plants able to mine resources and make them available for any 30 other plant in the community, termed ecosystem engineers. We use the model to assess whether and under which environmental conditions antagonistic facilitation via soil resource engineering is evolutionarily stable.
3. We find that antagonistic facilitation is stable in highly stressful conditions, which supports the theory of ecosystem engineers as drivers of primary succession and 35 provides a theoretical ground to investigate facilitation mechanistically in the context of the stress gradient hypothesis.
4. Among all potential causes of stress considered in the model, the key environmental parameter driving changes in the interaction between plants is the proportion of the limiting resource available to plants without mining. This finding represents a 40 challenge for empirical studies, which usually measure the resource input or loss in the system as a proxy for stress. We also find that the total root biomass and its spatial allocation through the root system, often used to measure the nature of the interaction between plants, do not predict facilitation reliably. 5. Synthesis. Antagonistic facilitation established between an ecosystem engineer nurse 45 plant and neighbor opportunistic individuals can be evolutionarily stable in stressful environments where ecosystem engineers' self-benefits from mining resources outweigh the competition with opportunistic neighbors. These results align with theories of primary succession and the stress gradient hypothesis as they show that antagonistic facilitation is stable under environmental stress, but it evolves into 50 mutual interference in milder environments. However, using inaccurate parameters to measure facilitation and stress gradients in empirical studies might mask these patterns.
Once the ecosystem engineer has established, opportunistic individuals that could not colonize the landscape by themselves can establish and proliferate (Verdú, Gómez, Valiente-Banuet, & Schöb, 2021). Following this proliferation, both types of individuals might interact via a non-reciprocal cooperative interaction in which the ecosystem engineer increases the 70 fitness of the opportunist (facilitation), and the opportunist reduces the fitness of its benefactor (interference) via competition for limited resources. This type of interaction is known as antagonistic facilitation (Schöb, Prieto, Armas, & Pugnaire, 2014).
Antagonistic facilitation is frequent (Brooker et al., 2008;Soliveres, Smit, & Maestre, 2015) and can determine the evolution of plant communities (Michalet et al., 2011). 75 However, the evolution and stability of antagonistic facilitation pose an evolutionary dilemma, which makes its ubiquity in nature puzzling. Some models have shown that facilitation among plants only stabilizes if ecosystem engineers cluster and engage locally in mutualistic interactions, hence getting a payback via direct reciprocity sensu Nowak, 2006(Kéfi, Van Baalen, Rietkerk, & Loreau, 2008Travis, Brooker, Clark, & Dytham, 2006). 80 These results imply that plant facilitation is unstable when ecosystem engineers are surrounded by opportunistic neighbors. Other studies conclude that ecosystem engineers must escape antagonistic facilitation by evolving mechanisms to avoid helping their neighbors (Bronstein, 2009;Brooker et al., 2008). Finally, a third group of studies have suggested that facilitation can only persist when the competitive effect of opportunistic plants 85 is so limited and context-dependent that it has no impact on the evolutionary trajectory of ecosystem engineers (Schöb, Callaway, et al., 2014). In this limit, the net effect of the opportunistic plant on the facilitator individual is practically neutral and hence not antagonistic. All these three families of studies question the evolutionary persistence of plant antagonistic facilitation.

90
Understanding the diversity of mechanisms by which a plant can facilitate neighbors' success is crucial to determining the evolutionary fate of antagonistic facilitation. Several, independent mechanisms may override resource competition and lead to net positive effects of one plant over another (Hunter & Aarssen, 1988). Facilitation can occur in any situation where plants experience strain and are, therefore, far from the optimal physiological 95 conditions for their development (Liancourt & Dolezal, 2021). A plant can reduce strain for its own benefit and the benefit of neighbors by driving any environmental parameter value toward its physiological optimum. Plants can, for example, engineer their environment to increase the availability of a scarce resource, such as water or nutrients (Lozano, Armas, Hortal, Casanoves, & Pugnaire, 2017;Francisco I. Pugnaire, Armas, & Valladares, 2004).
The root exudation of organic acids, chemical signals, and enzymes can improve soil 110 quality in several ways. Root exudations can dissolve unavailable soil nutrients such as phosphorus (P), iron (Fe), or potassium (K), avoid the toxic effects of hypoxia or an excess of aluminum (Al), and allow the establishment of symbiotic relations with nitrogen (N)fixing microorganisms (Dakora & Phillips, 2002;Hinsinger et al., 2011;D. L. Jones, 1998).
Because N and P are essential resources for plants and are typically scarce in soils, plants 115 that can boost the availability of these macronutrients have received considerable attention (Koffel, Boudsocq, Loeuille, & Daufresne, 2018;Lambers, Raven, Shaver, & Smith, 2008).
Plants adapted to grow in P-deficient soils can develop expensive structures, known as cluster roots, that allow them to release the P that has been sorbed to soil particles and make it available for locally foraging roots (Britto Costa, Staudinger, Veneklaas, Oliveira, & 120 Lambers, 2021;Raven, Lambers, Smith, & Westoby, 2018). Some plants growing in Nlimited soils can also develop root nodules and establish mutualistic relations with N-fixing microorganisms that trade N for photosynthates (Sprent, 1989;van Velzen, Doyle, & Geurts, 2019). Moreover, some plants can increase soil moisture by having root systems with specific characteristics. For example, plants with roots able to break superficial soil compaction can 125 increase local water infiltration to the soil (Bromley, Brouwer, Barker, Gaze, & Valentin, 1997;Montaña, 1992), and plants with deep tap roots reaching the water table can increase water availability in superficial soil layers by hydraulic lift (Caldwell & Richards, 1989;Prieto, Armas, & Pugnaire, 2012;Zapater et al., 2011). In all these situations, plant ecosystem engineers have specific root traits that allow them to locally ameliorate soil conditions at a 130 cost to themselves, creating islands of fertility where opportunistic plants can potentially grow.
Determining the environmental conditions and underlying biophysical mechanisms that turn antagonistic facilitation evolutionarily stable is crucial for supporting ecological theories that strongly depend on it, such as the stress gradient hypothesis (SGH). Some 135 ecosystems, such as deserts, Mediterranean shrublands, or tropical savannas, are more stressful and poorer in resources than others, such as temperate or tropical forests, even when successionally mature. The SGH predicts that facilitation must dominate in more stressful habitats while interference must prevail in mild conditions (Bertness & Callaway, 1994).
Because meta-analyses of studies reporting facilitation over stress gradients did not find 140 consistent support for this theory (Maestre, Valladares, & Reynolds, 2005), several authors suggested that stress could be classified as resource-stress and non-resource-stress and that only the former leads to the SGH predictions (He, Bertness, & Altieri, 2013;Maestre, Callaway, Valladares, & Lortie, 2009). Alternatively, the humped-back SGH suggests that plant facilitation maximizes at moderate, rather than highest, levels of stress (Holmgren & 145 Scheffer, 2010; Maestre et al., 2009;Michalet et al., 2006). Predicting how the evolutionary stability of antagonistic facilitation based on resource dynamics changes across environmental gradients is essential to improve our theoretical understanding of the SGH.
On the other hand, validating the SGH experimentally requires controlling for gradients of environmental stress and measuring one plant's effect on a neighbor's fitness.

150
Both these tasks require using reliable surrogates measurable in nature, which might be challenging to define. First, the concept of environmental stress is complex (Körner, 2003;Phillips, 1986) and compromises the empirical study of facilitation and the SGH (He & Bertness, 2014;Liancourt, Le Bagousse-Pinguet, Rixen, & Dolezal, 2017;Richard Michalet, Schöb, Lortie, Brooker, & Callaway, 2014). Changes in resource availability can define 155 environmental gradients for plants (e.g., Laliberté, Zemimik, & Turner, 2014), but these must be measured carefully. For example, some studies use the resource input into the soil, such as fertilization experiments or rainfall gradients, as surrogates of resource availability (e.g., Dohn et al., 2013;Wilson & Tilman, 1993). Other studies, however, use the rate of resource loss, such as nutrient leaching or soil drying rates (e.g., Cabal & Rubenstein, 2018;Pugnaire 160 & Luque, 2001). Both resource input and loss rates determine resource dynamics and can explain its scarcity, but they represent different mechanisms and hence may lead to different plant foraging behaviors (Cabal, Martínez-García, De Castro, Valladares, & Pacala, 2021).
Second, given the elusive and non-tangible nature of fitness (Costa, 2013;Primack & Hyesoon Kang, 1989), empirical studies need to rely on indirect measures to quantify how a 165 focal plant influences the fitness of a neighbor. The most common fitness surrogates used by plant biologists are survival rates, reproductive allocation, or biomass production, to name the most common. In particular, many studies consider that biomass production is a good proxy for fitness (Younginger, Sirová, Cruzan, & Ballhorn, 2017) and compare the total biomass (e.g., Dohn et al., 2013) or the spatial distribution of root biomass (e.g., Britto Costa 170 et al., 2021) growing alone vs. interacting to quantify fitness consequences of plant interactions. These studies assume that more biomass or allocation of biomass toward the neighbor when interacting as compared to alone indicates that plants are being facilitated. A mechanistic understanding of antagonistic facilitation based on resource dynamics would help to determine which of the parameters that can be measured in the field or controlled 175 experimentally are better to identify facilitation and define stress gradients.
Empirical studies reporting facilitation among plants in specific systems are copious, but we still lack a solid theoretical ground establishing general principles (Zhang & Tielbörger, 2020). In particular, we lack mechanistic models that unravel the biophysical processes sustaining net positive interactions among plant individuals and how those interactions 180 change across environmental gradients (Brooker et al., 2008;Soliveres et al., 2015). Here,

Plain text summary of the model
We model the below-ground root growth of plant individuals sharing neighborhood areas and competing for a single limiting resource. We assume that the total resource input 195 splits into two pools: unavailable and available resource, and consider two types of plants: soil ecosystem engineers and opportunistic individuals. Both types of plants forage from the available resource pool (e.g., water that infiltrates in the soil), and ecosystem engineers can mine from the unavailable pool and render resources available for all plants (Figure 1a). We further assume that the fraction of unavailable resource (e.g., runoff water) is lost if not mined 200 instantaneously and therefore does not accumulate over time (e.g., plant-triggered increase in infiltration). The available resource also leaves the soil at a constant rate due to physical losses (e.g., percolation, evaporation). Given this resource dynamics, plants spread their root system (root biomass per volume unit of soil) to maximize their net resource gain by balancing the benefits obtained from the resources they can forage and the root production 205 costs. The amount of resource plants can forage depends on the input and physical loss of available resource, total mining intensity, and foraging by other plants. Root production costs increase with the distance to the plant stem for any type of plant (representing the need to grow longer and thicker transportation roots) and, for ecosystem engineers, with mining intensity.

210
Complete solutions of this model provide both the optimal root distribution for single or multiple plants and, for ecosystem engineers, the optimal value of the resource mining trait. We calculate both quantities at equilibrium, when plants are fully grown and producing additional roots would decrease their net resource gain. We first obtain the optimal spatial distribution of root density that maximizes the net resource gain at each soil location. Next, 215 we use this root density distribution to calculate the net resource gain for a plant at each point of soil. Third, integrating the spatial distribution of root density and the net resource gain across the extension of the root system, we calculate the plant-level root biomass and plantlevel resource net gain, respectively. Lastly, we calculate the optimal mining intensity as the trait value that maximizes the plant-level net resource gain for the root density distribution 220 obtained in the first step.
We can obtain these model outputs for different setups and parameterizations. For example, we run the model using both a single plant and a pair of plants and compare the plant-level net resource gains of a chosen focal plant in these two scenarios. Based on this comparison between plant-level net resource gains in different growing environments, we 225 determine whether plants interact via facilitation (the plant-level net resource gain of the focal plant is higher when it grows paired than when it grows alone) or interference (otherwise). Additionally, for the ecosystem engineers, we assess the evolution of the mining trait under different environmental conditions and for engineer plants growing alone or in the presence of opportunistic plants. Lastly, we explore how root spatial distributions and 230 optimal mining intensities vary with the different environmental parameters: total resource input, proportion of such resource that is available spontaneously without mining, and physical rate at which the resource leaves the soil. Different environmental parameterizations represent different scenarios that can be interpreted as stress gradients in space or across successional time. Therefore, through its various variables and parameters, our model 235 accounts for processes occurring at different time scales: plants' phenotypically plastic root distribution, the evolution of the ecosystem engineer trait, and changes in the abiotic conditions across an environmental gradient (that can be defined in space or time) (Figure   1b). 240 We extend the model for resource competition introduced in Cabal, Martínez-García, De Castro Aguilar, Valladares, & Pacala, 2020 to investigate the evolution of plant antagonistic facilitation. Because we model individual root growth and spatial allocation, we use a spatially explicit modeling approach. The model assumes a series of processes governing the dynamics of a soil resource W and provides the net resource gain for each plant 245 due to the balance between how much resources it uptakes at a soil location, modeled via a classical consumer-resource interaction term, and the cost of foraging in that soil patch. It therefore consists of a series of consumer-resource models at each soil patch, with roots from different plants competing for a single limiting resource. We consider two types of plants:

Mathematical model formulation
Opportunistic plants, that can only forage soil resources, and plant ecosystem engineers that 250 consume soil resource and also its availability by evolving a resource mining trait.
Resource dynamics. We model resource dynamics in each soil point as a combination of three processes: input at rate I, abiotic loss at rate , and resource uptake by foraging roots R at a per capita rate  (see Table 1 for a summary of the environmental parameters and their dimensions). For mathematical simplicity, we consider a linear, non-saturating resource 255 uptake term and that both ecosystem engineer and opportunistic plants have the same per capita uptake rate. More complex scenarios, which are not the focus of this study, could account for plant type-specific per capita uptake rates and nonlinear, saturating uptake terms, including tradeoffs in resource uptakes between plant types (Grover, 1991;Tilman, 1982;Vincent, Scheel, & Brown, 1996). To incorporate the effect of ecosystem engineers, we model the resource input as a combination of spontaneously available resource and unavailable resource that is temporarily available for mining. In soils lacking ecosystem engineering roots, only a fraction b of all the potential resource,  is spontaneously available to roots (0  b  1). Ecosystem engineer roots RE increase the fraction of  that is available for plant foraging due to a per capita mining intensity . We model this positive feedback as which saturates at  when the mining intensity tends to infinity to account for the finiteness of resources. Putting all these terms together, the resource concentration changes at each soil location according to where we have dropped the dependence on space and time from the terms on the right side of Eq.
(2) and ℓ is the spatial coordinate of the soil location measured from the insertion point 270 of the focal plant's stem to the soil surface.
Net resource gain. At each soil location, the net resource gain of an individual plant, j, at time t is the balance between the resource uptake and the cost associated with growing and maintaining the roots needed to forage in that soil location, where R is the foraging root density (units of fine root biomass per unit of soil volume), C is 275 the cost of producing and maintaining such root density, and WUE is the resource use efficiency that represents the conversion factor from harvested resource to new root biomass.
Because the cost function depends on the distance between the soil patch and the plant insertion point to the soil surface, the model is spatially explicit. We account for three contributions to the cost function: the cost of producing and maintaining fine roots that absorb 280 resources from the soil cr, the cost of building and maintaining transportation roots that connect the fine roots to the plant stem ct, and the cost per density of fine roots associated to a unit increase in the mining intensity ce For opportunistic plants,  = 0 and the cost function in Eq. (4) has only two contributions (see Table 2 for a summary of the plant parameters) (Cabal et al., 2020).

Two-step model analysis.
Three parameters determine soil resource physical dynamics in our model and therefore define the environmental stress: the potential resource input , the fraction of such resource that is spontaneously available to plant roots b, and the resource physical loss rate . For example, if we consider that water is the limiting resource,  could represent 290 precipitation at a given location; b and 1-b, the fraction of the water input that infiltrates in bare soil and is lost by runoff, respectively; and  the loss rate due to percolation to the subterranean water We conducted a two-step computational analysis for different values of these parameters. In the first step, we obtained the resource mining trait value that optimizes the plant-level net resource gain and is thus expected to have evolved in those environmental conditions. We used this result to parameterize the model using a resource mining intensity To perform all the computational analysis, because the resource dynamics is much 310 faster than root growth, we assumed that the resource concentration is in quasi-equilibrium in each step of plant growth. This approximation allows us to write a closed expression for the net resource gain of each plant at a given soil location (see SM for a complete derivation), where RE is the total root biomass density of the ecosystem engineer at location ℓ. To optimize G in a given environmental scenario, we look for combinations of root density profiles and 315 resource mining intensity that optimize the net resource gain in Eq. (5).

First step: model parameterization.
To mimic natural conditions in which ecosystem engineer plants might have evolved the resource mining trait, we first allowed to evolve either assuming that ecosystem engineers grew alone, I, or competing with opportunistic plants, II. We obtained these values of  for 320 a range of environmental stresses, e.g., from b = 0 (high stress) to b = 1 (low stress). This first analysis yields the evolutionarily stable values of the resource mining trait for an ecosystem engineer plant at a given environmental stress level and assuming that the ecosystem engineer adapted to growing alone or surrounded by opportunistic plants.
Modeling a community considering several individual opportunistic and ecosystem 325 engineer plants in interaction with one another would make the model mathematically intractable and require many arbitrary choices about community size and spatial structure.
Instead, to get a general competitive background for the evolution of II, we modeled the opportunistic individual as a spreading plant that grows vegetatively covering the soil surface (de Kroons & Hutchings, 1995). These spreading opportunistic plants distribute their stem 330 everywhere on the soil surface, like a stolon or a rhizome, so their roots grow across the canopy area of the ecosystem engineer without any spreading costs (Dolezal et al., 2019).
Mathematically, this growth strategy implies removing the costs associated with transportation roots, ct = 0. We consider this assumption a more realistic and tractable representation of the average background community to which ecosystem engineers may 335 have adapted. We provide a detailed description of how this analysis was performed in the Supplementary Material.

Second step: simulated experiments.
In a second step, we considered individual plants from the above-parameterized conditions, i.e., using the values I or II expected to evolve at given environmental where X is  or depending on whether we measure changes in plant-level root biomass or 360 net resource gain, respectively. Xinteract is the value obtained for  or when the opportunistic plant grows close to an ecosystem engineer and Xsolo is the value obtained for solitary plants.
RII can vary between -1 and 1. RII informs about the net interaction, with RII = −1 indicating that the plant cannot coexist with the neighbor (competitive exclusion); −1 < RII < 0 net interference; RII = 0 no effect; 0 < RII < 1, net facilitation; and RII = 1 indicating 365 that the plant cannot exist without the neighbor (obligatory facilitation). For root biomass measurements, RII < 0 (RII > 0) indicates that opportunistic plants under-proliferate (overproliferate) roots when interacting with an ecosystem engineer relative to plants growing alone.
We hypothesize that the net plant interaction changes more substantially over an 370 environmental gradient defined by changes in b because this parameter controls how much of the total resource input becomes readily available for plants and how much goes to the temporarily mineable resource pool and is accessible only to soil engineer roots (Fig. 1a).
Hence, b-stress causes a strain that can be reduced by soil engineer roots. On the other hand,  and  determine the total input and loss rate, and create a strain that ecosystem-engineer 375 roots cannot mitigate. Therefore, we first analyze how the net plant interaction changes across an environmental gradient defined by changes in b and then conduct a more extensive scanning of the three-dimensional parameter space also to quantify the role of  and  in shaping the pairwise plant interaction.

Evolution of the mining trait and resulting net interactions with changes in b
Only ecosystem engineers with a mining intensity above a given threshold can survive alone in stressful habitats where the proportion of resource that is spontaneously available, b, is below a critical value bc-o. Neither ecosystem engineers investing in resource mining below that threshold nor opportunistic plants, for which the mining trait =0, can survive because their root production and maintenance costs outweigh resource uptake.
Opportunistic plants can invade the ecosystem in this regime only at the shelter of ecosystem engineers, hence establishing an obligatory interaction with them (Figure 2, I As increases and more resource is spontaneously available for plants, the resource 425 that ecosystem engineers can mine becomes progressively lower and II decreases. This decrease in the optimal mining intensity also reduces the difference between the plant-level net resource gain of solitary ecosystem engineer and solitary opportunistic plants ( Figure   SM2). When b reaches a critical value bc-e, the optimal mining intensity becomes II = 0, indicating that ecosystem engineers evolve to lose the resource mining trait and become

The different effects of environmental parameters on plant interactions and root proliferation
Through the three-dimensional parameter space defined by b, the physical loss of However, II changes abruptly from 0 at high levels of and -stress, to a positive value II 440 > 0 when stress levels in both parameters go below a threshold (Fig 4a, b).
The RII coefficient for the plant-level net resource gain of opportunistic plants RII shows that the interaction transitions gradually from facilitation at low values of b (high bstress) to interference at higher values of b. In contrast, at very high and -stresses not even engineers can survive (Fig 4c, d) because the stress is induced by resource loss that 445 cannot be mined. As the stress induced by these two parameters decrease, there is a first environmental threshold (solid green lines in Fig. 4c, d) Fig. 4c, d). After passing that second threshold, the net interaction, measured by RII , depends strongly on b but barely changes with  and . Only for a very narrow range of intermediate values of b, the interaction changes slightly from facilitation to 455 competition as -stress decreases, but from competition to facilitation as -stress decreases.
We observe that the fraction of resource that is available without the intervention of mining, , rather than the potential resource input  or the physical decay rate of the resource , accounts for most of the variation in the sign and the strength of the biotic interaction in the environmental region when pairs of ecosystem engineer/opportunistic plants interact.

460
The analysis of root production using RIIR shows that the transition from facilitation to competition (dashed golden lines in (active mining of the ecosystem engineer) or for II = 0 (both plants behave as purely exploitative, i.e., opportunistic). In addition to that, we found that a skewed distribution of 470 the opportunistic plant's root density in space toward the ecosystem engineer can also be observed in both scenarios, facilitation and interference (Fig SM3).

Successional changes in plant biotic interactions
For facilitation via soil resource engineering, the amount of resource available to plants without the intervention of root mining is the leading force driving the change in plant interactions across environmental gradients. Only ecosystem engineers can colonize highly stressed ecosystems, which could represent the initial stages of succession at which opportunistic plants cannot survive on their own. Nevertheless, after ecosystem engineers colonize the system, opportunistic plants can establish an obligatory interaction with their 480 benefactors and invade the ecosystem. Previous work has suggested that ecosystem engineers may lose their mining trait over the course of evolution due to the presence of opportunistic neighbors that benefit from mined resources at no cost (Koffel et al., 2018;Song, Yu, Jiang, Korpelainen, & Li, 2019;Walker & Chapin III, 1987). Our results, however, suggest that antagonistic facilitation is evolutionarily stable at high environmental stresses. In this regime,

Consequences for the stress gradient hypothesis
Our model further dissects stress induced by resource availability into three parameters; the total resource input  (e.g., levels of precipitation), the fraction of such resource that becomes available to plants spontaneously b (e.g., how much of the rainwater 500 infiltrates in bare soil and becomes available to roots without the intervention of ecosystem engineers), and the rate of physical resource decay  (e.g., how fast the infiltrated water leaves the soil due to evaporation or percolation). Facilitation can occur only if there is strain, and ecosystem engineers can reduce such strain. Contrarily, if the physical environment is near optimal, or if plants cannot improve its quality, facilitation is not possible. Our results 505 support the SGH for resource stress if the stress is caused by the fraction of the total resource that is available for plants without mining b, but we did not find a significant impact of the other environmental parameters, the physical decay rate  and total resource input , on the net interaction. This is because, even though all parameters can produce strain, only under bstress resource pools can be mined by ecosystem engineers to reduce strain. As a result, b is

Consequences for empirical studies
Our results have also important consequences for the interpretation of empirical 525 studies addressing interaction shifts across environmental gradients. We have identified two major sources of potential confusion based on our model. First is the fact that b is the main environmental driver of shifts in plant interactions. This parameter, which could represent, for example, the ratio between water infiltration and runoff, is much harder to determine empirically than the other two environmental parameters in our model,  or  (e.g., rainfall 530 or soil drying rates). Many empirical studies assess the change in plant interaction across environmental gradients using the latter two parameters, which our model suggests could lead to misleading conclusions. Second, our calculations of the plant-level net resource gain suggest that plant biomass is a misleading fitness surrogate to determine the net effect of one plant on another. Specifically, root over-proliferation or root density skewed in space towards 535 a neighboring ecosystem engineer might not be a univocal consequence of facilitation, as both behaviors may be observed when interacting negatively with the neighbor due to a tragedy of the commons (Cabal, 2022;Gersani, Brown, O'Brien, Maina, & Abramsky, 2001). This result represents a challenge for experimental designs, because the plant-level net resource gain is not easy to determine empirically. In cases where this measure is not 540 possible, other fitness proxies, such as resource allocation into reproductive tissues, may be better fitness indicators (Kozłowski, 1992;Younginger et al., 2017).

Conclusions
Antagonistic facilitation enables more diverse plant communities in stressful environments at the initial stages of primary succession and may accelerate the development of mature

600
The authors have no conflicts of interest to declare Author contributions.  Figure 1: a. Schematic representation of soil resource dynamics in a specific soil location in our model. This dynamics is replicated at each point of space labeled with ℓ in the model equations. We consider a constant input of potential resource () that splits into spontaneously available to plants (b) and temporarily minable (1-b). Ecosystem engineer roots can mine the latter and make it available to all plants at a rate proportional to the mining 890 intensity () and the ecosystem engineer root density (RE). Every plant can forage the available resource W (in quasi-equilibrium) at a rate proportional to their per-capita absorption rate () and root density (RE for the ecosystem engineer and RO for the opportunistic). The available resource also leaves the soil at a physical decay rate (). b. The different scales at which we assessed the emergence of facilitative interactions, including the 895 plastic response of plant root distribution R at the life span time scale (top), the adaptive response of the ecosystem engineers' facilitator trait  at the evolutionary time scale (center), and the change in the fraction of soil resource that is spontaneously available over succession b, which can change through space or at a geological time scale (bottom).

Figure 2:
Graphical summary of the changes across an environmental gradient defined by the proportion of resource that is spontaneously available to plant roots, b (see Figure SM2 for detailed results obtained numerically from model equations). Across the gradient, b varies from 0 (high stress, none of the resource is available without mining) to 1 (low stress, all the resource is spontaneously available). The figure shows the evolution of the mining trait  in