Synthetic soil crusts against royalsocietypublishing.org/journal/rsos green-desert transitions: a spatial model Blai Vidiella1,2, Josep Sardanyés3,4 and Ricard V. Solé1,2,5 Research 1 ICREA-Complex Systems Lab, Universitat Pompeu Fabra, 08003 Barcelona, Spain Cite this article: Vidiella B, Sardanyés J, Solé 2 Institut de Biologia Evolutiva (CSIC-UPF), Psg. Maritim Barceloneta, 37, 08003 Barcelona, RV. 2020 Synthetic soil crusts against green- Spain 3 desert transitions: a spatial model. R. Soc. Open Centre de Recerca Matemàtica, and 4Barcelona Graduate School of Mathematics (BGSMath), Edifici C, Campus de Bellaterra, 08193 Cerdanyola del Vallès, Bellaterra, Barcelona, Spain Sci. 7: 200161. 5 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM 87501, USA http://dx.doi.org/10.1098/rsos.200161 BV, 0000-0002-4819-7047; JS, 0000-0001-7225-5158; RVS, 0000-0001-6974-1008 Semiarid ecosystems are threatened by global warming due to Received: 5 February 2020 longer dehydration times and increasing soil degradation. Accepted: 22 July 2020 Mounting evidence indicates that, given the current trends, drylands are likely to expand and possibly experience catastrophic shifts from vegetated to desert states. Here, we explore a recent suggestion based on the concept of ecosystem terraformation, where a synthetic organism is used Subject Category: to counterbalance some of the nonlinear effects causing Ecology, Conservation, and Global Change Biology the presence of such tipping points. Using an explicit spatial model incorporating facilitation and considering Subject Areas: a simplification of states found in semiarid ecosystems bioengineering/ecology/synthetic biology including vegetation, fertile and desert soil, we investigate how engineered microorganisms can shape the fate of these Keywords: ecosystems. Specifically, two different, but complementary, synthetic biology, ecological engineering, climate terraformation strategies are proposed: Cooperation-based: C-terraformation; and Dispersion-based: D-terraformation. The change, catastrophic shifts, mutualism, synthetic first strategy involves the use of soil synthetic microorganisms ecology to introduce cooperative loops (facilitation) with the vegetation. The second one involves the introduction of engineered microorganisms improving their dispersal capacity, Author for correspondence: thus facilitating the transition from desert to fertile soil. We show that small modifications enhancing cooperative loops Ricard V. Solé can effectively modify the aridity level of the critical transition e-mail: ricard.sole@upf.edu found at increasing soil degradation rates, also identifying a stronger protection against soil degradation by using the D-terraformation strategy. The same results are found in a mean-field model providing insights into the transitions and dynamics tied to these terraformation strategies. The potential consequences and extensions of these models are discussed. Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c. 5082867. © 2020 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. 1. Introduction 2 royalsocietypublishing.org/journal/rsos Global warming is changing the dynamics and resilience of ecosystems, damaging many of them and creating the conditions for widespread diversity loss [1–3]. Because of the presence of nonlinearities, many ecosystems could suffer the so-called tipping points [4–6], responsible for community collapses and for transitions towards more degraded states [7,8] (e.g. coral bleaching [9] or lake eutrophication [10]). Among these systems, drylands (which comprise arid, semi-arid and dry-subhumid ecosystems) are a specially fragile subset of major importance: they include more than 40% of terrestrial ecosystems and host a similar percentage of the current human population [11,12]. Increasing aridity is pushing these ecosystems towards serious declines in microbial diversity, land degradation and loss of multifunctionality as desert states are approached [11,13,14]. Dedicated efforts have been addressing several avenues to both understanding how transitions can be anticipated by means of warning signals [15–17] and even prevented [11,18,19]. Catastrophic transitions i.e. the transition that occurs at the tipping point, from vegetated to desert ecosystems can be really fast [19], even in large R. Soc. Open Sci. 7: 200161 ecosystems [20]. This type of sudden transitions have been identified in the evolution of the vegetated Sahara to the current desert state, which occurred less than 6000 years ago [21–23]. Drylands are characterized by the presence of organisms that have adapted to low moisture availability, damaging UV radiation, and high temperatures [24–26]. A rich community structure and the maintenance of physical soil coherence are essential to prevent drylands from degradation [27–29]. In this context, some universal types of interactions that occur among species in arid habitats inevitably lead to breakpoints associated with the existence of multiple alternative states, identified in field data [14,30]. A well-known class of these interactions takes place among vascular plants and is known as facilitation [31–33]. Facilitation consists on non-trophic interactions between individuals mediated through changes in the abiotic environment favouring individual growth and reproduction [34–36]. An example of this process is the fertile islands in the savannah ecosystems, where the presence of trees increase the soil organic carbon thus increasing soil moisture and nutrients [33]. This process facilitates the implantation and growth of new vegetation (grasses and trees) [37,38]. A common outcome of spatially extended ecosystems is the emergence of spatial patchiness e.g. peatlands (see [39] and references therein) or the Kalahari desert [38], to cite a few. These spatial patterns have been extensively studied using computational models [39,40], for instance, models only considering competition [41]. Other models including both competition and facilitation processes also typically display self-organized patchiness [38,42–45]. Such patterns are often remarkably organized in space [45]. Transitions between different ecological states, e.g. vegetated towards bare soil, can be sometimes observed in the same landscape under the presence of environmental gradients [31]. In this context, it has been conjectured that spatial correlations can be used as indicators of forthcoming green-desert shifts [15,46,47]. The presence of catastrophic transitions deeply modifies our perception of risks associated with climate change and land degradation. Once a tipping point is reached, unstoppable runaway processes are unleashed [8]. Despite much progress having been made in modelling drylands [47–49] as well as in identifying warning signals [15–17,46], feasible strategies to prevent green-desert transitions are currently scarce. In a recent paper [42], it has been shown that, by tuning some particular features of models exhibiting catastrophic shifts, one can modify their nature and location in parameter spaces. Based on a theoretical approach, the authors suggested that the amount of stochasticity (both demographic and extrinsic) or population dispersal range can play that role. Indeed, the impact of plant dispersal ranges on the nature of transitions has been recently investigated in a simple metapopulation model with facilitation [43], showing that short-range dispersal might involve continuous transitions. Examples of gradual regime shifts in spatial patterns have been also described theoretically [44]. What type of microscopic mechanisms could actually avoid green-desert transitions? In this paper, we are extending previous models of dryland dynamics [19,47] to predict the impact of bioengineering strategies on these fragile ecosystems. In this context, several bioremediation approaches have been developed in the last two decades to enhance and stabilize soil biocrusts (for an overview, see [26] and references therein). These include a diverse range of approaches addressed to improve moisture and soil texture, from enrichment of key nutrients [50] or enrichment by cyanobacteria [51–53], to large-scale straw chequerboard barriers used to stabilize sand dunes [54]. More recently, ecosystem terraformation has been suggested as a novel approach against green-desert shifts [55–59]. In a nutshell, some existing microorganisms in a degraded ecosystem, such as cyanobacteria of the soil crust, could be minimally modified by means of synthetic biology techniques to help improve soil moisture and create a cooperative feedback between vegetation or moss cover (see 3 [60]) and the soil microbiome [61]. This engineering scenario aims at building synthetic soil composed of royalsocietypublishing.org/journal/rsos synthetic strains mixed with natural microorganisms embedded, for example, in the soil crust. By doing so, potential tipping points could be made much more difficult to be achieved [58]. What might be the impact of these strategies on the large-scale, long-term dynamics of semiarid ecosystems? How can new or enhanced ecological interactions (cooperation and microbes dispersal capacities) modify the presence of tipping points? As shown below by means of both mathematical and computational models, the engineering of microorganisms present in arid and semiarid soils may easily expand the potential of ecosystems’ persistence. 2. Modelling terraformation of semiarid ecosystems In order to predict the impact of synthetic bioengineering strategies in arid and semiarid ecosystems, a R. Soc. Open Sci. 7: 200161 spatial model given by a stochastic cellular automaton (CA) is employed. Let us first consider the microscopic rules associated with the dynamics, as described by a set of probabilistic transitions among the three defined ecological states. After this general description of the transitions, we will investigate a mean-field model without stochasticity and the stochastic CA. The mean-field model will allow us to identify the nature of the tipping points and the expected equilibria (which may be useful for low stochasticity levels i.e. large system’s size). The basic transition diagram is displayed in figure 1c, which includes the three potential states characterizing a given patch, namely S ¼ {D, S, V}: These states, as mentioned, are defined by desert (D) patches, by fertile soil (S) composed by patches occupied by engineered microorganism(s) with their natural community (e.g. soil crust), and by vegetation (V), respectively. This is, of course, an oversimplification that ignores most of the complexity and diversity involved, but allows for an analysis of the dynamics considering key ecological interactions such as transitions between states and both competition and facilitation processes. The model described here is derived from the one proposed by Kéfi et al. [47], now considering that unoccupied sites (non-degraded and non-vegetated) contain soil crust plus (engineered) microorganisms. As discussed in a previous work [58], we aim at describing how an appropriate synthetic modification of microorganisms could help to maintain the stability and resilience of semiarid ecosystems. The first rule, which considers the transition from degraded to fertile soil containing synthetic strains, will take place with a probability P(D ! S) ¼ r þ f rv þ G(rs ): Here, r is the probability of spontaneous recovery of fertile soil due to external factors such as increased humidity, accumulation of organic matter, etc. Since we are interested in identifying the impact of the engineered strains in the recovery of the desert soil to the fertile one we will set r = 0. This will allow us to provide clearer results concerning soil recovery derived from terraformation. Constant f denotes the facilitation due to the influence of the surrounding vegetation. The variable ρv is the local density of vegetation (number of vegetated neighbours over the total number of neighbours). Finally, Γ(ρs) denotes the spreading capacity of the microorganisms associated with their synthetic engineered properties (see below). The second transition considers the reverse situation, namely a fertile soil-to-desert transition P(S ! D) ¼ 1 Y(rv ), where ɛ is the probability of loss of fertile soil due to increased aridity or to other soil degradation processes. In this paper, we consider (see below) a modification of the degradation rate mediated by a vegetation-dependent function Y(rv ). The third transition rule involves the colonization of available S patches by vegetation P(S ! V) ¼ (dDV þ (1  d)rv )(b  cDV ): (2:1) This term has the same form as the colonization used in [47]. The probability δ balances the influence of the local and global vegetation to produce the germination of new plants in a given site. Seeds are produced and germinate with probability b. However, they can also be degraded during the (a) (b) 4 royalsocietypublishing.org/journal/rsos (c) (d) V D-terraformation S R. Soc. Open Sci. 7: 200161 C-terraformation S D V positive interaction D negative interaction Figure 1. Ecological structure of semiarid ecosystems and modelling approach. (a) Semiarid landscape from Aranjuez (central Spain). The enlarged view shows the patchy structure of vegetation, surrounded by bare soil and soil crust (shown in (b)). Our computational model considers three well-differentiated ecological states: vegetation (V), functional soil (including soil crust, S) and bare soil (desert sites, D). The local dynamics follows a set of stochastic rules that allow transitions between the three different states, as indicated in (c) using black arrows. Moreover, the presence of vegetation cover has an impact on these transitions, as indicated by the links ending in black circles (e.g. vegetation provides organic matter to the soil which increases water retention and the microorganisms in there degrade this matter for vegetation uptake). These are the basic rules used in previous literature modelling semiarid ecosystems [19,46,47]. The proposed scenario of synthetic terraformation considered in this article is indicated using coloured links. The orange line denotes C-terraformation (a synthetic strain establishing cooperative feedbacks with vegetation), inhibiting the transition from S to D. The blue line indicates D-terraformation, allowing for increased dispersal of the engineered organisms (SYN) thus favouring the transition form D to S. (d) We use a stochastic cellular automaton including these patch types, here displayed in the lattice projected on a sphere. dispersion (with probability c). Finally, an exponential decay of vegetated patches to fertile soil occurs following the transition rule P(V ! S) ¼ m, where m is the death probability of plants, which can occur due to pathogens, grazing, or to the recollection by livestock activity. Here two different, but complementary, terraformation strategies are considered and weighted by means of constants α and β. Their role within the transition diagram is displayed in figure 1c. They involve: 1. C-terraformation: Engineering cooperative loops between vegetation and soil organisms. Here, a decreasing function of the vegetation cover, namely 1 Y(rv ) ¼ , 1 þ arv introduces a reduction in the impact of desertification due to increased soil quality (by e.g. nutrient deposition) favoured by the synthetic microbial population. The efficiency of this term is weighted by the constant α: large values of α imply a lower soil degradation rate due to the action of the synthetic microbes which are cooperatively coupled with vegetation at a local scale. 2. D-terraformation: Engineering the capacity of microbial spreading (e.g. faster replicative rates of the microbes and/or increased formation of endospores able to successfully colonize local surroundings). The dispersal is here considered proportional to the term 5 G(rs ) ¼ b rs : royalsocietypublishing.org/journal/rsos Here, ρs is the local density of fertile soil. In other words, with probability β, engineered strains improve their dispersal capabilities, thus changing the properties of desert patches by retaining humidity, depositing organic matter, etc. Setting α and β to zero, the CA model recovers the original model introduced by Kéfi et al. [47]. 3. Models and results In the next sections, we will explore the dynamics tied to the rules described above using both a well-mixed (mean-field) and a discrete spatial setting using a stochastic Cellular Automaton (CA). The first approach, developed in §3.1, employs differential equations to predict the presence of tipping points and changes in R. Soc. Open Sci. 7: 200161 their location in the parameter space due to the effects of the addition of a synthetic strain. The second approach, studied in §3.2, explicitly deals with a spatially extended population where local interactions take place on a two-dimensional lattice and the set of rules are applied in a probabilistic manner, thus considering stochastic effects. Finally, we provide insights into the nature of the tipping points found in the spatial model, focusing on the early warning signals for the system without and with the synthetic engineered microorganisms. 3.1. Mean-field model The differential equations model, which includes the interactions and processes tied to the bioengineered synthetic strains, is given by dV ¼ V((b  cV)S  m), (3:1) dt dS ¼ (fV þ bS)D  1 S Y  V((b  cV)S  m) (3:2) dt and dD ¼ 1 Y S  (r þ fV þ G)D: (3:3) dt Equation (3.1) contains a logistic term that includes variable S as a multiplicative term, indicating that plants need viable soil to persist, and an exponential decay proportional to m. Equation (3.2) includes the positive effects triggered by vegetation and the existing soil cover, as well as negative terms that are in fact symmetric to those present in the previous equation for V. Here recall that Y ¼ 1=(1 þ a V). Equation (3.3) provides the dynamics of desert patches resulting from desertification, recruitment to viable soil patch, and an exponential decay from S to D introduced with Γ = β S. Since the three classes of patches cover the entire lattice in the spatial model, it is possible to normalize and consider the states as population fractions in the mean field approach, i.e. V(t) + S(t) + D(t) = 1. This actually assumes a constant amount of available sites that can transition between them, also allowing to reduce the three-variables system to a two-variables one using the linear relation D(t) ¼ 1  V(t)  S(t), (3:4) which needs to satisfy dV/dt + dS/dt + dD/dt = 0 (recall that r = 0) since the sum of all states remains constant through time. The reduced system is thus given by dV ¼ V((b  cV)S  m) (3:5) dt and   dS 1 ¼ (fV þ bS)(1  S  V)  1 S  V((b  cV)S  m): (3:6) dt 1 þ aV Note that now the fraction D is automatically obtained from equation (3.4) once the fractions of states V and S are determined. We must note that equations (3.5) and (3.6) have been recently studied in [19] for the non-engineered system (that is, taking α = β = 0). We are specially interested in those scenarios involving a full dominance of the desert state, focusing in the transitions between states. We have identified four equilibria for equations (3.5) and (3.6), two of them implying the extinction of vegetation and another one allowing for the stable coexistence of the 6 three ecological states. The first equilibrium is at the origin, i.e. P1 ¼ (V ¼ 0, S ¼ 0), were the system royalsocietypublishing.org/journal/rsos becomes a desert (D = 1). The second equilibrium, keeps the fertile soil and the desert patches without vegetation, and is given by P2 ¼ (V ¼ 0, S ¼ 1  1=b). This state is interesting from the point of view of the terraformation strategies, since its existence depends on the chosen engineering strategy. If β = 0, when no engineered organisms are found in the soil crust, this equilibrium is not biologically meaningful (S(t → ∞) = −∞). Two more equilibrium points, labelled P3 and P4 , have been identified numerically (see electronic supplementary material, figures S1–S4 and S6–S8). The equilibrium P4 (numerical results suggest it is always stable within the simplex; see below) involves V(t → ∞) > 0, S(t → ∞) > 0 and V(t → ∞) + S(t → ∞) < 1, allowing the coexistence of the three ecological states when it is an interior equilibrium. P3 is a saddle point. See below and electronic supplementary material, figures S1–S4 for further details on the dynamics of the equilibria on the simplex (V, S). The conditions defining the local stability of each equilibrium can be obtained by means of the R. Soc. Open Sci. 7: 200161 eigenvalues of the Jacobian matrix J , given, for equations (3.5) and (3.6), by 0 1 m þ (b  2 c V)S V(b  c V) J ¼@ A: m  f(1 þ 2V þ S) þ S (b þ 2cV þ 1a(1 þ V a)2  b)  V(b þ f  c V)  1(1 þ V a)1 þ b  (V þ 2 S)b The eigenvalues of matrix J evaluated at an equilibrium provide its local (and linear) stability properties. The stability for equilibria P1,2 can be easily computed, while the stability of equilibria P3,4 will be characterized numerically. In the case of the desert equilibrium P1 ¼ (V ¼ 0, S ¼ 0), the eigenvalues are λ1 = −m and λ2 = β − ɛ. Since all the parameters are positive, this equilibrium will be stable when ɛ > β. For P2 ¼ (0, 1  1=b), its eigenvalues are λ1 = ɛ − β and λ2 = −m + b(1 − ɛ/β); see electronic supplementary material, figure S5 for a stability diagram for P2 in the parameter space (ɛ, β). Note that P1 and P2 suffer a transcritical bifurcation at β = ɛ, since the two conditions for this bifurcation are fulfilled: the two fixed points collide at the bifurcation value (at ɛ = β, P1 ¼ P2 ¼ (0, 0)), and they interchange the stability. Figure 2 displays the critical degradation rate of the fertile soil, ɛc, computed numerically in the parameter space (α, β), separating two domains with and without vegetation at equilibrium. The value of ɛc moves to higher values as either α or β are increased. The yellow region in the surface of ɛc in figure 2a denotes those values of α and β where εc ≥ 1 (although they are all set to ɛc = 1 for further comparison with the probabilistic results of §3.2). This region denotes that no critical degradation is achieved under the chosen parameter values and under the restriction 0 ≤ ɛ ≤ 1. The four panels in figure 2b (corresponding to the values of α and β indicated in figure 2a with the same letters) provide a more detailed picture of the changes in the fraction of the states due to both C- and D-terraformation strategies. Panel (b.a) shows how the states change without terraformation (α = β = 0) at increasing ɛ. Note that for ɛ > 0.2 the desert state becomes dominant. This critical degradation rate is displaced to larger values for the terraformed system, meaning that the ecosystem becomes much more resistant to soil degradation. More specifically, ɛc ≈ 0.3 when applying the C-terraformation strategy ( panel (b.b) in figure 2b). This effect is further amplified for the D-terraformation, resulting in ɛc ≈ 0.7 (see panel (b.c) in figure 2b). The combination of both terraformation strategies (see figure 2b panel (b.d)) further increases the size of non-desert phases, here with ɛc ≈ 0.8. In both cases (specially for D-terraformation), the domain of fertile soil is increased due to the presence of the synthetic strains (electronic supplementary material, figures S6–S8 show how the values of ɛc are displaced for different intensities of the proposed terraformation strategies). The dynamics tied to increases of soil degradation rate are summarized in figure 2c (see also electronic supplementary material, figures S1–S4), where different phase portraits are displayed for the values of ɛ indicated in figure 2b panel (b.c) (from (i) to (iv)). The transition from the coexistence of vegetation and fertile soil to no vegetation is catastrophic, and is due to a saddle-node bifurcation between fixed points P3 and P4 (transition between phase portraits (ii) and (iii), see also electronic supplementary material, figures S1–S4). Once this bifurcation takes place, the only stable point is P2 (with fertile soil and desert areas). Then, further increase of ɛ involves the transcritical bifurcation previously mentioned, after which the origin becomes globally stable and the desert becomes the only possible state (transitions between phase portraits (iii) and (iv) of figure 2c), see also electronic supplementary material, figure S7 where the transcritical bifurcation is shown. Finally, electronic supplementary material, figure S6 provides one-dimensional bifurcation diagrams tied to the increase of ɛ for the non-terraformed system and the three possible terraformation strategies critical rate of 7 (a) fertile soil loss, ec (b) V : vegetation royalsocietypublishing.org/journal/rsos S : fertile soil 0 0.5 1 D : desert soil a 1 b 1 d 1.0 S S fraction 0.8 D D 0.6 V V ec 0 0 c 0.4 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 c 1 d 1 0.2 D D fraction 0 S S R. Soc. Open Sci. 7: 200161 10 on 8 ati 6 er 1.0 0.8 a op V iv V 4 0 i ii iii co 0.6 0 2 in 0.4 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 ase i b incre se 0.2 0 fertile soil loss rate, e fertile soil loss rate, e ea n dis 0 cr persio in n (c) 1 (i) 1 (ii) 1 (iii) 1 (iv) S S S S 0 0 0 0 0 1 0 1 0 1 0 1 V V V V Figure 2. Dynamics of the mean-field model for the non-terraformed ecosystem and the terraformed one. (a) Critical degradation rate of fertile soil (ɛc) computed in the parameter space of engineering strategies (α, β) from equations (3.5) and (3.6). The inner surface separates the values of ɛ allowing for the presence of vegetation. The yellow area indicates those pairs of α and β giving εc ≥ 1 (those values of ɛc > 1 are set to ɛc = 1). The region with no extinction of vegetation is also displayed projected on the ground of the plot (white zone). (b) Fraction of the states at equilibrium increasing ɛ and using a full vegetated system as initial conditions V(0) = 1, S(0) = D(0) = 0. (b.a) Results for a non-engineered ecosystem (α = β = 0), and (b.b) for the engineered ecosystem incorporating cooperation loops (α = 1 and β = 0). Finally, panels (b.c) and (b.d) display, respectively, the results engineering the resilience of the soil crust (α = 0 and β = 1) and both strategies (α = β = 1). The dashed vertical lines indicate the critical value of fertile soil loss obtained for the non-engineered system εc = 0.218 [19]. (c) Phase portraits for the case α = 0 and β = 1, with: (i) ɛ = 0.0; (ii) ɛ = 0.2; (iii) ɛ = 0.6; and (iv) ɛ = 0.9. The circles indicate the fixed points (stable: solid; unstable: open). We note that in the phase portraits (i) and (ii) there exists an interior saddle (see electronic supplementary material, figures S1–S4 for the identification of the nullclines and dynamics). Green, blue and red regions of the phase portraits denote the equilibrium states of V-S coexistence, S-D coexistence, and full desert, respectively. The other parameters are r = 0, f = 0.9, δ = 0.1, b = 0.6, c = 0.3 and m = 0.15. (C-, D- and (C, D)-terraformation). Here, we want to emphasize the presence of so-called delayed transitions (also called ghosts), which arise just after a saddle-node bifurcation takes place [62]. This is a dynamical phenomenon that involves extremely long transients once the bifurcations has occurred, and the time trajectories experience a long bottleneck before rapidly achieving, in equations (3.5) and (3.6), another attractor (the full desert state in electronic supplementary material, figure S6(a) and (b)). These long transients [63] are typically found in systems with strong feedbacks, such as cooperation, catalytic processes [64,65] or metapopulations with facilitation [43]. Also, this dynamical delay tied to saddle-node bifurcations has been recently described in both deterministic and stochastic well-mixed approaches for the non-terraformed system explored in this article (see [19] for more details). 3.2. Spatial stochastic model The mean-field model studied in the previous section provides a first approximation to understand the qualitative dynamics arising from the nonlinear interactions of the studied system and the resulting tipping points, especially when including the proposed terraformation strategies. However, in order to test the robustness of these results in a more ecologically realistic setting, one has to take into account 8 both local spatial correlations and stochastic effects [66]. To do so, we build a spatially explicit royalsocietypublishing.org/journal/rsos simulation model given by a stochastic cellular automation (CA). The CA incorporates the previously studied states (V, S, D) and transitions among them, which are now probabilistic. Spatial degrees of freedom are introduced by using a L  L square lattice with periodic boundary conditions. At each time step, the following four transition probabilities defining the rate of change for each site k of the lattice (figure 1) are applied in an asynchronous manner: 1 X 1 X P(Dk ! Sk ) ¼ r þ f Vm þ b Sm , MN meM MN meM k k d P(Sk ! Dk ) ¼ P , 1 þ (a=MN ) meMk Vm ! ! d X (1  d) X bc P(Sk ! Vk ) ¼ Vm þ Vm  P , MN meM R. Soc. Open Sci. 7: 200161 L meV 2 k L2 meV Vm P(Vk ! Sk ) ¼ m: Here, Mk is the neighbourhood of site k and MN ¼ 8 the number of neighbour (using a Moore neighbourhood). The nature of local (spatial) interactions is known to largely influence ecological dynamics [40,43,67]. This is particularly important in drylands, where carbon and water limitation deeply constrains the outcome of nonlinear exchanges [14,59], leading to spatial patterning [48,68]. Also, facilitation processes (involving strong nonlinearities) are known to introduce important changes in spatial systems, as compared with well-mixed ones. In this sense, recent research has found a shift from catastrophic tipping points to continuous ones, due to local spatial processes [42,43]. Similar analyses to those presented in figure 2 for the mean-field model are displayed in figure 3 for the CA simulations. For the sake of comparison, we have kept the same parameter values of the mean- field model, implemented as probabilities in the CA, also using as initial conditions a full vegetated lattice. The critical surface separating the parameter scenarios allowing the persistence of vegetation is displayed in figure 3a. Below the surface, V persists, while above the surface it is only possible to find states S and D. Here, similarly to the mean-field model, the yellow region indicates those values of α and β for which no critical degradation rate of the fertile soil is found (i.e. ɛc = 1). Note that the surface for the spatial system differs from the one obtained with the mean-field model. This effect, taking into account that the lattice has 4 × 104 sites (large system size), is probably introduced by space more than by stochasticity. This change of the surface is especially visible for parameter α. This may be explained because soil degradation rate (ɛ) makes the vegetation decrease globally, but due to the local interactions of facilitation [33], plants remain in the ecosystem for larger rates of soil degradation. The fraction of the three states computed at increasing ɛ is displayed in figure 3b. All diagrams (b.a–b.d) display the same phenomenon: the introduction of spatial correlations shifts the extinction of the vegetation to larger values of ɛ. Consistently with mean-field results, both the C- and D-terraformation displace the value of ɛc to larger values, the (C, D)-terraformation strategy being much more efficient in doing so under the selected parameters. The spheres displayed in figure 3c show the spatial patterns at equilibrium for the regions labelled in panel (c) of figure 3b, with: (i–ii) coexistence of the three states; (iii) only fertile soil and desert patches; and (iv) the full desert scenario. In order to identify the nature of the transitions tied to the increase of ɛ, we have computed the stationary values of the fraction of areas with vegetation and fertile soil. We recall that both mean field and well-mixed stochastic simulations for the system without terraformation revealed a catastrophic transition and delayed transitions due to a ghost [19]. Also, the mean-field model studied in §3.1 has revealed the presence of saddle-node bifurcations responsible for vegetation extinctions. Figure 4 shows these results together with the spatial patterns on the square lattice and time series. Specifically, figure 4a shows results increasing ɛ for the system without terraformation. For this case ɛc ≈ 0.24. The same results are displayed by setting α = β = 1, for which the critical degradation moves to ɛc ≈ 0.85. For both cases, the discontinuity of the transition is not so evident as in the mean-field model, even looking like a continuous one. To identify the transition governing vegetation extinction, we have used an indirect method analysing the properties of the dynamics close to the transition value. The results displayed in electronic supplementary material, figure S9 indicate the presence of bottlenecking phenomena tied to the critical probability of 9 (a) (b) V : vegetation fertile soil loss, ec S : fertile soil royalsocietypublishing.org/journal/rsos D : desert soil 0 0.5 1 a 1 b 1 1.0 S S fraction 0.8 D D 0.6 V V ec 0 0 0.4 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 0.2 c 1 d 1 D D 0 S fraction 10 S 8 on R. Soc. Open Sci. 7: 200161 ati 1.0 6 er 0.8 4 a op V V 0.6 co 0.4 2 0 i ii iii iv in incre b 0.2 0 0 se ase i 0 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 ea n dis probability of probability of cr pers in ion fertile soil loss, e fertile soil loss, e (c) (i) (ii) (iii) (iv) Figure 3. Dynamics of the spatial model (using a 200 × 200 sites lattice) with and without terraformation strategies (see figure 2 for comparison). (a) Vegetation extinction surface depending on the engineering strategies α and β. The yellow region indicates that vegetation extinction cannot be achieved because ɛ = 1 (this region is also projected on the ground of the plot (white zone)). (b) Bifurcation diagram using the probability of fertile soil loss ɛ for different engineering strategies (α β). Here, the dashed vertical lines also indicate the critical value of fertile soil loss obtained with the mean-field model for the non-engineered system, εc = 0.218 · · ·. (b.a) Non-engineered ecosystem (α = β = 0), (b.b) engineering of cooperative loops between synthetic microorganisms and the vegetation (α = 1 and β = 0), (b.c) engineering resilience of the soil crust (α = 0 and β = 1); and (b.d ) engineering both terraformation strategies (α = β = 1). (c) Spatial patterns for the case α = 0 and β = 1 ( panel (b.c) above) with: (i) ɛ = 0.00; (ii) ɛ = 0.55; (iii) ɛ = 0.75; and (iv) ɛ = 0.90. The probability values have been fixed as the parameters of figure 2, also with V(0) = 1. extinction of the vegetation, meaning that the transition is governed by a saddle-node bifurcation. We must notice that we have used a lattice of side size L ¼ 50 for these analyses due to the huge computational cost of computing extinction transients in large spatial systems. Electronic supplementary material, figure S9(a) shows the transition for the non-engineered system. Specifically, panel (a.1) displays five realizations obtained with ɛ = 0.26, and the bottleneck region can be clearly seen (see electronic supplementary material, figure S6 for comparison with the mean field bottlenecks). The delaying effect of the ghost can also be seen in panels (a.2), (a.3) and (a.4), which indicate the values at which vegetation spends longer iterating before achieving the full desert state. The same results are shown for the system with (C, D)-terraformation. Here, the extinction dynamics of the vegetation (investigated setting ɛ = 0.86) is also given by a clear bottleneck. Previous research in catalytic hypercycles also revealed extinctions due to saddle-node bifurcations and the corresponding bottlenecking phenomena in a stochastic CA model [65]. As discussed above, the approach to tipping points is often associated with a qualitative change in the nature of the fluctuations exhibited by the system. How will our model behave when dealing with an extra phase associated with vegetation-free, synthetic soil crust? A divergence in a stochastic dynamical system can be detected by using standard deviation measures. We have used several measures that can detect the approach to criticality in both the non-manipulated and the terraformed scenarios (spatial variance and distribution of patch sizes). The spatial variance computed as standard (a) 0.68 10 fertile soil (S) 0.305 vegetation royalsocietypublishing.org/journal/rsos 0.30 0.67 0.295 0.6 0.66 0 5 × 102 103 iterations 0.31 fraction of area 0.58 fertile soil (S) vegetation 0.302 0.4 0.57 0.56 0.294 0 5 × 102 103 iterations 0.22 0.30 0.2 fertile soil (S) vegetation R. Soc. Open Sci. 7: 200161 0.29 0.21 d = 0.25 d = 0.15 d = 0.05 0.28 0.27 0.2 0 5 × 102 103 iterations 0 0 0.2 0.4 0.6 0.8 1.0 probability of fertile soil loss, e (b) 0.60 fertile soil (S) vegetation 0.31 0.59 0.6 0.30 0.58 0 5 × 102 103 iterations fraction of area fertile soil (S) vegetation 0.47 0.4 (C, D)-terraformation 0.31 0.46 0.45 0.30 0 5 × 102 103 iterations 0.2 0.23 fertile soil (S) vegetation 0.16 d = 0.65 d = 0.85 0.22 d = 0.3 0.15 0.21 0 0 5 × 102 103 iterations 0 0.2 0.4 0.6 0.8 1.0 probability of fertile soil loss, e Figure 4. Equilibrium population values for the vegetation (V) and the fertile soil (S) at increasing ɛ obtained from the cellular automaton (CA) using a lattice with 200 × 200 sites. For each value of ɛ, the fraction of areas are plotted taking the last 200 values from time series with 106 generations repeated over 10 realizations (to illustrate the variance of the fluctuations). The darkest dots correspond to the last population value of each realization. In (a), we display the results for the non-engineered system (with α = β = 0). Three examples of the spatial patterns at the end of a single realization are displayed for three values of ɛ, together with time series (using the same colours as in the spatial patterns) obtained at the end of the simulations. (b) The same information as in (a) for the engineered system (α = 1 and β = 1). Note the huge displacement of the critical soil degradation probabilities between panels (a) and (b). deviation σψ, with ψ referring to either V or S, is determined from !1=2 1 XL 2 hsc i ¼ (c(k, t)  hc(t)i)2 , (3:7) L2 k where we have used the mean over the entire lattice, given by 1 XL 2 hc(t)i ¼ c(k, t): (3:8) L2 k (a) ·sy Ò (c) (d) 11 1 0.05 1 0.05 0.05 0.10 non-engineered 0.10 royalsocietypublishing.org/journal/rsos 0.15 0.20 0.04 0.20 0.30 –0.75 original 0.25 –0.75 0.40 0.03 10–1 10–1 0.50 0.60 0.02 V 0.70 0.80 0.01 10–2 10–2 S count 0 (b) 0.05 original terraformed 10–3 10–3 0.04 terraformed S S 0.03 V 10–4 V 10–4 V 0.02 0.01 D D S 10–5 10–5 0 0 0.2 0.4 0.6 0.8 1.0 1 10 102 103 104 1 10 102 103 104 probability of fertile soil loss, e patch size patch size R. Soc. Open Sci. 7: 200161 Figure 5. Warning signals for the non-engineered (α = β = 0) and terraformed (α = 1 and β = 1) ecosystems taking vegetation and fertile soil as indicators at increasing ɛ. Panels (a–b) display the standard deviation 〈σψ=V,S〉 of the population fluctuations for the original and the terraformation scenarios, corresponding to the simulations shown in electronic supplementary material, figure S9. Values of 〈σψ〉 have been obtained from the last 103 values from time series with 106 iterations. (c–d ) Cumulative distributions of patch sizes for different values of ɛ (indicated with different colours inside the plots). The patch size distributions have been computed using 100 replicas from a lattice with 200 × 200 sites. The fluctuations diverge in a characteristic fashion, as figure 5a,b shows: whereas 〈σV〉 slowly grows and accelerates close to criticality, the same measure for 〈σS〉 displays a slight decline only growing very quickly close to ɛc. Despite these trends, both variables exhibit high fluctuations close to the transition to the desert state. An additional confirmation of the criticality associated with the transition to the desert state is obtained by determining the distribution of vegetation cluster sizes. In a nutshell, criticality is known to be linked to power law distributions of connected clusters. One particular instance of these scaling behaviours was found in [47]. Clusters of size S are composed by a set of S connected sites (i.e. sharing the same state and all elements in contact with another as nearest neighbours). The number of clusters of a given size and their distributions have been obtained using a burning algorithm [69]. If P(S) indicates the probability of finding a cluster of size S, most measured patch size distributions found in arid and semiarid habitats follow a general form described as a truncated power law, namely P(S)  S t exp (S=S c ). This defines a power law decay with a scaling exponent τ that characterizes the statistical organization of systems close to percolation [69] along with an upper bound introduced by a cut-off S c in the exponential decay term. The cut-off defines the crossover from a regime of critical clusters ( power-law distributed) to that of non-critical clusters. As we get closer to percolation points, the system also displays divergent correlation lengths and S c ! 1 thus rapidly increasing the cut-off value and moving closer to pure power laws. Since the lattice is finite, we will always observe a rapid decay at the end of the distribution. These distributions can be systematically estimated from vegetation data [38,70,71]. Close to critical transitions we should expect a fat-tailed decay in P(S) following the power law form P(S)  S t (when no characteristic scale is present, i.e. for S c ! 1). Instead, far below the transition point, the exponential term dominates and a characteristic scale is observed. In order to improve the smoothness of the distribution, the so-called cumulative distribution P. (S) will be used, and is given by ð Sm P. (S) ¼ P(S) dS, (3:9) S0 where S 0 and S m indicate the smallest (the one-site scale) and largest (lattice) sizes. At criticality, the cumulative form scales as a power law P. (S)  S tþ1 : (3:10) In our model, the shapes of the distributions change at increasing values of ɛ: from single-scaled to scale- free, as shown in figure 5c (original model) and figure 5d (terraformed system). In both cases, the same pattern is found, but the power-law behaviour is observed at much higher levels of soil degradation. Interestingly, the scaling exponent for both cases is the same: τ ≈ 1.75, in agreement with the results reported in Kéfi et al. [46]. The consideration of the extra processes tied to terraformation do not alter a universal behaviour pattern, thus suggesting that, despite the parameter shift in ɛc, the warning signals 12 associated with both scenarios will behave in a similar way. royalsocietypublishing.org/journal/rsos 4. Discussion The impact of anthropogenic-driven processes is rapidly increasing, leading to a degradation of extant ecosystems and, in many cases, pushing them closer to viability thresholds [1,2,72,73]. In the case of semiarid ecosystems, this degradation is produced mostly by global warming and intensive grazing [22]. The increasing pressure could end up in catastrophic transitions that are essentially irreversible. It is expected that the expansion of arid areas will increase in the next decades, with an associated increase in the likelihood of green-desert tipping points. Under this potential scenario, it is important to develop strategies of intervention aimed at avoiding these transitions. Theoretical results on the nature of these shifts and the parameter values at which they take place R. Soc. Open Sci. 7: 200161 suggest that some generic factors (such as noise and dispersal [42,43]) could be crucial. Can realistic interventions help preventing green-desert transitions? As discussed above, several bioremediation strategies to face this problem have been proposed. Some recent proposals suggest how to exploit nonlinear features after tipping points based on periodic replanting [19]. More recently, an engineering approach based on synthetic biology has also been suggested [56–58]. In this article, we investigate this latest approach for arid and semiarid ecosystems employing both mathematical and computational models. Beyond the standard vegetation-desert model, in our study we pay attention to the possibility that an engineered soil crust (a ‘synthetic’ soil) might have on the fate of fragile or endangered ecosystems. The original motivation is to describe, using a toy model, the impact of two different but complementary strategies: the so-called C-terraformation using a designed cooperative loop ( promoted by designed microbial strains); and the D-terraformation, consisting on engineered strains with increased dispersal capacities. Cooperative loops (included with the C-terraformation at the soil crust level in our modelling approach) have indeed been identified in natural ecosystems. For instance, termites activity may modify the nature of transitions to desert states by increasing water infiltration rates and soil moisture [74]. In this manuscript, we have explored the impact of these two synthetic terraformation strategies on the location and behaviour of tipping points. We show that both terraformation strategies increase the resistance to soil degradation, pushing the tipping points to higher (or even much higher) critical degradation rates (ɛc). This is consistently predicted by both mean-field and spatial models. However, when space is explicitly included, the parameter space supporting this protection against catastrophic shifts is much larger, with a fluctuation dynamics (and underlying warning signals) that behave similarly between the non- engineered and the engineered systems, i.e. they show the same universality pattern. Our theoretical predictions are limited by the simplifying assumptions that define our model. In particular, the description of soils ignores their development, spatial organization or diversity, as well as the network structure (e.g. trophic relations) of ecosystems [75]. These webs are known to display their own critical thresholds, and thus it remains open how such diverse communities may impact the results reported here. However, despite all these limitations, it is important to highlight that similar previous models of vegetation dynamics in drylands have been very successful in providing deep insight into their natural counterparts [15]. Some key universal properties might pervade the success of these models. In our context, this is likely to be related to the general patterns found close to phase transitions. Additionally, the terraformation scheme described here is likely to be effective when the cooperative interaction helps create the proper ecosystem engineering effect. By engineering ecosystem engineers, endangered arid and semiarid habitats might get protected (at least temporarily) from sudden collapse. Finally, we may emphasize that ecosystems microbiome is a complex ecosystem itself within the full ecosystem. The introduction of a synthetic species could induce secondary effects at this level that could scale up towards higher ecological (e.g. trophic) levels. However, the synthetic biology approaches suggested here and in previous research [56–58] may re-use already existing microorganisms, which could be synthetically modified. The impact of a truly new species in the ecosystem may be investigated elsewhere by e.g. invasion (adaptive dynamics) or complex networks theory. Data accessibility. The programming code for the cellular automaton model is provided as a electronic supplementary material. Authors’ contributions. R.V.S. and B.V. conceived the model. B.V. and J.S. wrote the codes of the mean-field model. B.V. 13 wrote the codes of the cellular automaton. All authors analysed the mathematical and computational models. All royalsocietypublishing.org/journal/rsos authors wrote the article and gave final approval for publication. Competing interests. Authors have no competing interests. Funding. This study was supported by an European Research Council Advanced Grant (SYNCOM), by the Botin Foundation (Banco Santander through its Santander Universities Global Division), the PR01018-EC-H2020-FET- Open MADONNA project, by the FIS2015-67616-P grant, and by the Santa Fe Institute. This work has also received the support of Secretaria d’Universitats i Recerca del Departament d’Economia i Coneixement de la Generalitat de Catalunya. J.S. has been funded by a ‘Ramón y Cajal’ contract RYC-2017-22243, and by the MINECO grant no. MTM2015-71509-C2-1-R and the Spain’s ‘Agencia Estatal de Investigación’ grant no. RTI2018- 098322-B-I00, as well as by the CERCA Programme of the Generalitat de Catalunya. Acknowledgements. The authors specially thank Sergi Valverde for the spatial spheric rendering used for the cellular automaton. We also thank Nuria Conde, Joan V., Fernando T. Maestre and Miguel Berdugo for helpful discussions. R. Soc. Open Sci. 7: 200161 References 1. Barnosky AD et al. 2012 Approaching a state of 14. Berdugo M et al. 2020 Global ecosystem 181–189. (doi:10.1890/1540- shift in Earth’s biosphere. Nature 486, 52–58. thresholds driven by aridity. Science 367, 9295(2003)001[0181:TWAYFD]2.0.CO;2) (doi:10.1038/nature11018) 787–790. (doi:10.1126/science.aay5958) 25. Pointing SB, Belnap J. 2012 Microbial 2. Rockström J et al. 2009 A safe operating space 15. Kéfi S, Guttal V, Brock WA, Carpenter SR, Ellison colonization and controls in dryland ecosystems. for humanity. Nature 461, 472–475. (doi:10. AM, Livina VN, Dakos V. 2014 Early warning Nature Rev. Microbiol. 10, 551–562. (doi:10. 1038/461472a) signals of ecological transitions: methods for 1038/nrmicro2831) 3. Rothman DV. 2017 Thresholds of catastrophe in spatial patterns. PLoS ONE 9, 0092097. (doi:10. 26. Belnap J, Lange OL eds. 2003 Biological soil the Earth system. Sci. Adv. 43, e1700906. 1371/journal.pone.0092097) crusts: structure, function and management. (doi:10.1126/sciadv.1700906) 16. Scheffer M et al. 2009 Early-warning signals for Berlin, Germany: Springer. 4. Lenton TM, Held H, Kriegler E, Hall JW, Lucht critical transitions. Nature 461, 53–59. (doi:10. 27. Moore DC, Singer MJ. 1990 Crust formation W, Rahmstorf S, Schellnhuber HJ. 2008 Tipping 1038/nature08227) effects on soil erosion processes. Soil Sci. Soc. elements in the Earth climate system. Proc. Natl 17. Sommer S, van Benthem KJ, Fontaneto D, Ozgul Am. J. 54, 1117–1123. (doi:10.2136/sssaj1990. Acad. Sci. USA 105, 1786–1793. (doi:10.1073/ A. 2017 Are generic early-warning signals 03615995005400040033x) pnas.0705414105) reliable indicators of population collapse in 28. Belnap J, Büdel B. 2016 Biological soil crusts as 5. Steffen W, Rockström J, Richardson K, Lenton rotifers? Hydrobiologia 796, 111–120. (doi:10. soil stabilizers. In Biological soil crusts: an TM, Folke C, Liverman D, Schellnhuber HJ. 2018 1007/s10750-016-2948-7) organizing principle in drylands, pp. 305–320. Trajectories of the earth system in the 18. Sardanyés J, Duarte J, Januário C, Martins N. Cham, Switzerland: Springer. anthropocene. Proc. Natl Acad. Sci. USA 2012 Controlling delayed transitions with 29. Jones CG, Lawton JCG, Shachak M. 1997 Positive 115, 8252–8259. (doi:10.1073/pnas. applications to prevent single species and negative effects of organisms as physical 1810141115) extinctions. Adv. Differ. Equ. Control Proc. 10, ecosystem engineers. Ecology 78, 1946–1957. 6. Lenton TM, Rockström J, Gaffney O, Rahmstorf 29–41. (doi:10.1890/0012-9658(1997)078[1946: S, Richardson K, Steffen W, Schellnhuber HJ. 19. Vidiella B, Sardanyés J, Solé R. 2018 Exploiting PANEOO]2.0.CO;2) 2019 Climate tipping points—too risky to bet delayed transitions to sustain semiarid 30. Berdugo M, Kéfi S, Soliveres S, Maestre FT. 2017 against. Nature 575, 592–595. (doi:10.1038/ ecosystems after catastrophic shifts. J. R. Soc. Plant spatial patterns identify alternative d41586-019-03595-0) Interface 15, 20180083. (doi:10.1098/rsif.2018. ecosystem multifunctionality states in global 7. Lenton TM. 2013 Environmental tipping points. 0083) drylands. Nat. Ecol. Evol. 1, 0003. (doi:10.1038/ Annu. Rev. Environ. Resourc. 38, 1–29. (doi:10. 20. Cooper GS, Willcock S, Dearing JA. 2020 Regime s41559-016-0003) 1146/annurev-environ-102511-084654). shifts occur disproportionately faster in larger 31. Berdugo M, Soliveres S, Kéfi S, Maestre FT. 2018 8. Scheffer M. 2009 Critical transitions in nature ecosystems. Nat. Commun. 11, 1175. (doi:10. The interplay between facilitation and habitat and society. Princeton, NJ: Princeton University 1038/s41467-020-15029-x) type drives spatial vegetation patterns in global Press. 21. deMenocal P, Ortiz J, Guilderson T, Adkins J, drylands. Ecography 42, 1–13. 9. Bellwood DR, Hughes TP, Folke C, Nyström M. Sarnthein M, Baker L, Yarusinsky M. 2000 32. Bronstein JL. 2009 The evolution of facilitation 2004 Confronting the coral reef crisis. Nature Abrupt onset and termination of the African and mutualism. J. Ecol. 97, 1160–1170. (doi:10. 429, 827–833. (doi:10.1038/nature02691) Humid Period: rapid climate responses to 1111/j.1365-2745.2009.01566.x) 10. Scheffer M, Jeppesen E. 2007 Regime shifts in gradual insolation forcing. Quat. Sci. Rev. 19, 33. Zélé F, Duncan AB. 2018 Ecology and evolution shallow lakes. Ecosystems 10, 1–3. (doi:10. 347–361. (doi:10.1016/S0277-3791(99)00081-5) of facilitation among symbionts. Nat. Commun. 1007/s10021-006-9002-y) 22. Foley JA, Coe MT, Scheffer M, Wang G. 2003 9, 4869. (doi:10.1038/s41467-018-06779-w) 11. Reynolds JF et al. 2007 Global desertification: Regime shifts in the Sahara and Sahel: 34. Brooker RW et al. 2008 Facilitation in plant building a science for dryland development. interactions between ecological and climatic communities, the past, the present, and the Science 316, 847–851. (doi:10.1126/science. systems in Northern Africa. Ecosystems 6, future. J. Ecol. 96, 18–34. (doi:10.1111/j.1365- 1131634) 524–532. (doi:10.1007/s10021-002-0227-0) 2745.2008.01373.x) 12. Cherlet M, Hutchinson C, Reynolds J, Hill J, 23. Collins JA, Prange M, Caley T, Gimeno L, 35. Callaway RM, Pugnaire FI. 1999 Facilitation in Sommer S, von Maltitz G eds. 2018 World atlas Beckmann B, Mulitza S, Skonieczny C, Roche D, plant communities. In Handbook of functional of desertification. Luxembourg: Publication Schefuß E. 2017 Rapid termination of the plant ecology (eds. FI Pugnaire, F Valladares), Office of the European Union. African Humid Period triggered by northern pp. 623–648. New York, NY: Marcel Dekker. 13. Maestre FT et al. 2016 Structure and functioning high-latitude cooling. Nat. Commun. 8, 1–11. 36. Rubio A, Escudero A. 2000 Small-scale spatial of dryland ecosystems in a changing world. (doi:10.1038/s41467-016-0009-6) soil-plant relationship in semi-arid gypsum Annu. Rev. Ecol. Evol. and Syst. 47, 215–237. 24. Belnap J. 2003 The world at your feet, desert environments. Plant Soil 220, 139–150. (doi:10. (doi:10.1146/annurev-ecolsys-121415-032311) biological soil crusts. Front. Ecol. Environm. 1, 1023/A:1004764411116) 37. Rodriguez-Iturbe I, Chen Z, Stave AC, Levin SA. 50. Maestre FT, Martín N, Díez B, Lopez-Poma R, 62. Strogatz SH. 2000 Nonlinear dynamics and 14 2019 Tree clusters in savannas result from Santos F, Luque I, Cortina J. 2006 Watering, chaos. Boulder, CO: Westview Press. islands of soil moisture. Proc. Natl Acad. Sci. USA fertilization, and slurry inoculation promote 63. Hastings A, Abbott KC, Cuddington K, Francis T, royalsocietypublishing.org/journal/rsos 116, 6679–6683. (doi:10.1073/pnas. recovery of biological crust function in degraded Gellner G, Lai YC, Morozov A, Petrovskii S, 1819389116) soils. Microb. Ecol. 52, 365–377. (doi:10.1007/ Scranton K, Zeeman ML. 2018 Transient 38. Scanlon TM, Caylor KK, Levin SA, Rodriguez- s00248-006-9017-0) phenomena in ecology. Science 361, eaat6412. Iturbe I. 2007 Positive feedbacks promote 51. Park C, Li X, Zhao Y, Jia RL, Vur J-S. 2017 Rapid (doi:10.1126/science.aat6412) power-law clustering of Kalahari vegetation. development of cyanobacterial crust in the field 64. Sardanyés J, Solé R. 2006 Ghosts in the origins Nature 449, 209–212. (doi:10.1038/ for combating desertification. PLoS ONE 12, of life? Int. J. Bifurc. Chaos 16, 2761–2765. nature06060) e0179903. (doi:10.1142/S0218127406016446) 39. Eppinga MB, Rietkerk M, Wassen MJ, De Ruiter 52. Zheng Y, Xu M, Zhao J, Bei S, Hao L. 2011 65. Sardanyés J, Solé R. 2006 Bifurcations and PC. 2009 Linking habitat modification to Effects of inoculated Microcoleus vaginatus on phase transitions in spatially-extended two- catastrophic shifts and vegetation patterns in the structure and function of biological soil member hypercycles. J. Theor. Biol. 243, bogs. Plant Ecol. 200, 53–68. (doi:10.1007/ crusts of desert. Biol. Fertil. Soils 47, 473–480. 468–482. (doi:10.1016/j.jtbi.2006.07.014) s11258-007-9309-6) (doi:10.1007/s00374-010-0521-5) 66. Marro J, Dickman R. 1999 Nonequilibrium phase 40. Solé R, Bascompte J. 1997 Modeling 53. Bowker MA. 2007 Biological soil crust transitions in lattice models. New York, NY: spatiotemporal dynamics in ecology. New York, rehabilitation in theory and practice: an Cambridge University Press. R. Soc. Open Sci. 7: 200161 NY: Springer. underexploited opportunity. Rest. Ecol. 15, 67. Peters DC, Parsons AJ, Bestelmeyer BT, Herrick 41. Marítnez-García R, Calabrese JM, Hernández- 13–23. (doi:10.1111/j.1526-100X.2006.00185.x) JE, Wainwright J, Okin GS, Fredrickson EL. 2009 García E, López C. 2013 Vegetation pattern 54. Li XR, Xiao HL, He MZ, Zhang JG. 2006 Sand barriers Do changes in connectivity explain formation in semiarid systems without of straw checkerboards for habitat restoration in desertification? BioScience 59, 237–244. (doi:10. facilitative mechanisms. Geophys. Res. Lett. 40, extremely arid desert regions. Ecol. Eng. 28, 1525/bio.2009.59.3.8) 6143–6147. (doi:10.1002/2013GL058797) 149–157. (doi:10.1016/j.ecoleng.2006.05.020) 68. Porporato A, Dódorico P, Laio F, Ridolfi L, 42. Villa Martín PV, Bonachela JA, Levin SA, Muñoz 55. Maestre FT, Solé R, Singh BK. 2017 Microbial Rodriguez-Iturbe I. 2002 Ecohydrology of water- MA. 2015 Eluding catastrophic shifts. Proc. Natl biotechnology as a tool to restore degraded controlled ecosystems. Adv. Water Res. 25, Acad. Sci. USA 112, E1828–E1836. (doi:10.1073/ drylands. Microb. Biotech. 10, 1250–1253. 1335–1348. (doi:10.1016/S0309- pnas.1414708112) (doi:10.1111/1751-7915.12832) 1708(02)00058-1) 43. Sardanyés J, Piñero J, Solé R. 2019 Habitat loss- 56. Solé R. 2015 Bioengineering the biosphere? 69. Christensen K, Moloney NR. 2005 Complexity induced tipping points in metapopulations with Ecol. Complexity 22, 40–49. (doi:10.1016/j. and criticality, vol. 1. London, UK: Imperial facilitation. Pop. Ecol. 61, 436–449. (doi:10. ecocom.2015.01.005) College Press. 1002/1438-390X.12020) 57. Solé R, Montañez R, Duran-Nebreda S. 2015 70. Staver AC, Asner GP, Rodriguez-Iturbe I, Levin 44. Bel G, Hagberg A, Meron E. 2012 Gradual Synthetic circuit designs for earth SA, Smit IP. 2019 Spatial patterning among regime shifts in spatially extended ecosystems. terraformation. Biol. Direct 10, 37. (doi:10.1186/ savanna trees in high-resolution, spatially Theor. Ecol. 5, 591–604. (doi:10.1007/s12080- s13062-015-0064-7) extensive data. Proc. Natl Acad. Sci. USA 011-0149-6) 58. Solé R, Montañez R, Duran-Nebreda S, 116, 10 681–10 685. (doi:10.1073/pnas. 45. Klausmeier CA. 1999 Regular and irregular Rodriguez-Amor D, Vidiella B, Sardanyés J. 2018 1819391116) patterns in semi-arid vegetation. Science 284, Population dynamics of synthetic terraformation 71. Martín PV, Domínguez-García V, Muñoz MA. 1826–1828. (doi:10.1126/science.284.5421. motifs. R. Soc. Open Sci. 5, 180121. (doi:10. 2019 Intermittent percolation and the scale-free 1826) 1098/rsos.180121) distribution of vegetation clusters. (http://arxiv. 46. Kéfi S, Rietkerk M, Alados CL, Pueyo Y, 59. Conde-Pueyo N, Vidiella B, Sardanyés J, org/abs/1912.08369). Papanastasis VP, ElAich A, De Ruiter PC. 2007 Berdugo M, Maestre FT, De Lorenzo V, Solé R. 72. Barnosky AD et al. 2011 Has the Earth’s sixth Spatial vegetation patterns and imminent 2020 Synthetic biology for terraformation mass extinction already arrived? Nature 471, desertification in Mediterranean arid lessons from Mars, Earth, and the microbiome. 51–57. (doi:10.1038/nature09678) ecosystems. Nature 449, 213–217. (doi:10. Life 10, 14. (doi:10.3390/life10020014) 73. Ellis EC, Kaplan JO, Fuller DQ, Vavrus S, 1038/nature06111) 60. Delgado-Baquerizo M, Maestre FT, Eldridge DJ, Goldewijk KK, Verburg PH. 2013 Used planet: 47. Kéfi S, Rietkerk M, van Baalen M, Loreau M. Bowker MA, Jeffries TC, Singh BK. 2018 a global history. Proc. Natl Acad. Sci. USA 110, 2007 Local facilitation, bistability and transitions Biocrust-forming mosses mitigate the impact of 7978–7985. (doi:10.1073/pnas.1217241110) in arid ecosystems. Theor. Pop. Biol. 71, aridity on soil microbial communities in 74. Bonachela JA, Pringle RM, Sheffer E, Coverdale 367–379. (doi:10.1016/j.tpb.2006.09.003) drylands: observational evidence from three TC, Guyton JA, Caylor KK, Levin SA, Tarnita CE. 48. D’Odorico P, Porporato A (eds). 2006 Dryland continents. New Phytol. 220, 824–835. (doi:10. 2015 Termite mounds can increase the ecohydrology. Dordrecht, The Netherlands: 1111/nph.15120) robustness of dryland ecosystems to climatic Springer. 61. Park C, Li X, Jia RL, Vur J-S. 2015 Effects of change. Science 347, 651–655. (doi:10.1126/ 49. D’Odorico P, Okin GS, Bestelmeyer BT. 2012 superabsorbent polymer on cyanobacterial science.1261487) A synthetic review of feedbacks and drivers of biological soil crust formation in laboratory. Arid 75. Montoya JM, Pimm SL, Solé R. 2006 Ecological shrub encroachment in arid grasslands. Land Res. Manage. 29, 55–71. (doi:10.1080/ networks and their fragility. Nature 442, Ecohydrology 5, 520–530. (doi:10.1002/eco.259) 15324982.2014.928835) 259–264. (doi:10.1038/nature04927)