The structure of plant communities – their composition and diversity – forms the foundation of many ecosystem services on which human well-being crucially depends. These include provisioning services such as food, fresh water, wood, and fiber; regulating services such as flood regulation and water purification; cultural services such as recreation and aesthetic enjoyment; and supporting services such as soil formation, photosynthesis, and nutrient cycling (Duraiappah and Naeem, 2005). These services are at risk due to potential changes in the composition and diversity of plant communities as a result of global warming and the development of drier climates (Harrison et al., 2020). Understanding the factors that affect community structure in varying environments calls for integrated studies of mechanisms operating at different levels of organization, from phenotypic changes at the organism level, through intraspecific interactions at the population level, to interspecific interactions at the community level (Gratani, 2014; Falik et al., 2003; Bertness and Callaway, 1994; Pérez-Ramos et al., 2019). Of these mechanisms, the role of intraspecific interactions in driving community dynamics through spatial self-organization, has hardly been studied (see however Vandermeer and Yitbarek, 2012; Bonanomi et al., 2014; Cornacchia et al., 2018; Zhao et al., 2019; O'Sullivan et al., 2019; Pescador et al., 2020; Inderjit et al., 2021).

Spatial self-organization in periodic and non-periodic patterns of biota and resources, driven by intraspecific competition that leads to partial plant mortality, is widely observed in stressed environments (Rietkerk and van de Koppel, 2008). An important class of these phenomena are vegetation patterns in drylands. In sloped terrains, these patterns commonly appear as parallel vegetation stripes oriented perpendicular to the slope direction (Lefever and Lejeune, 1997; Valentin et al., 1999; Bastiaansen et al., 2018). In flat terrains, with no preferred direction imposed by slope or wind, stripe-like configurations often appear as labyrinthine patterns. However, in such terrains two additional pattern morphologies are often observed; nearly periodic patterns of bare-soil gaps in otherwise uniform vegetation, and nearly periodic patterns of vegetation spots in otherwise uniform bare soil (Rietkerk et al., 2004; Deblauwe et al., 2008; Borgogno et al., 2009; Getzin et al., 2016). Spatial self-organization may not necessarily result in periodic patterns; according to pattern-formation theory (Meron, 2015; Knobloch, 2015), it can also result in non-periodic patterns, such as single or randomly scattered vegetation spots in otherwise bare soil, randomly scattered bare-soil gaps in otherwise uniform vegetation and others (Tlidi et al., 2008; Dawes and Williams, 2016; Parra-Rivas and Fernandez-Oto, 2020; Jaïbi et al., 2020; Zelnik et al., 2015). The driving mechanisms of these self-organized vegetation patterns are scale-dependent positive feedback loops between local vegetation growth and water transport toward the growth location (Rietkerk and van de Koppel, 2008; Pringle and Tarnita, 2017; Meron, 2019).

Vegetation patterns involve not only spatial distributions of plant biomass, but also less-visible distributions of soil-water, nutrients, soil biota, and possibly toxic substances (Paz-Kagan et al., 2019; Inderjit and Duke, 2003; De Deyn et al., 2004; van der Putten et al., 2013). The various habitats that these self-organizing distributions form lead to niche differentiation and community reassembly (Weiher et al., 2011). Thus, spatial self-organization and community dynamics are intimately coupled processes that control community composition and diversity. Understanding this unexplored coupling is essential for assessing the impact of global warming and drier climates on community structure and ecosystem services.

In this paper, we incorporate spatial self-organization into community-assembly studies, using a mathematical model of dryland plant communities. Our model study provides three new insights, illustrated in Figure 1: (i) Spatial self-organization acts to reverse community-structure changes induced by environmental stress, (ii) it buffers the impact of further stress, and (iii) it offers new directions of ecosystem management that integrate the need for provisioning ecosystem services with the need to conserve community structure. More specifically, using a trait-based approach (Dı́az and Cabido, 2001), we show that drier climates shift the composition of spatially uniform communities toward stress-tolerant species, and reduce functional diversity, that is the diversity of functional traits around the most abundant one. By contrast, self-organization in spatial patterns, triggered by these droughts, shifts the composition back to less tolerant species that favor investment in growth, and increase functional diversity. Once patterns are formed, spatial self-organization provides various pathways to relax further stresses without significant changes in community composition and diversity. Furthermore, multistability ranges of uniform and patterned states open up new opportunities for grazing and foddering management by forming mixed community states of increased functional diversity.

Figure 1

Download asset Open asset

We use a continuum modeling approach where plant populations are described by above-ground biomass densities rather than by plant-number densities. This approach allows a continuum representation even of highly discrete plant populations, as often encountered in arid regions (Meron et al., 2019). The model is based on a model platform that has been introduced and discussed in earlier studies (Gilad et al., 2004; Gilad et al., 2007; Nathan et al., 2016; Meron, 2016), but involves several modifications. It consists of a system of partial differential equations (PDEs) for biomass variables, $B_i$ ($i=1,\mathrm{\dots },N$), representing the spatial densities of above-ground biomass of the $N$ functional groups, ${\chi }_1,{\chi }_2,\mathrm{\dots },{\chi }_N$, and two water variables, $W$ and $H$, representing spatial densities of below-ground and above-ground water, respectively, all in units of kg/m². We restrict ourselves in this work to one-dimensional spatial dependence of all state variables. Thus, solutions that are periodic in the $x$ direction are assumed to be independent of the orthogonal $y$ direction and represent periodic stripe patterns. The community model then reads:

where the second ‘trait derivative’, ${\partial }_{\chi }^2{B}_i\equiv N^2(B_{i+1}-2B_i+B_{i-1})$, represents mutations at a very small rate $D_{\chi }$. In these equations, the growth rate of the $i$ th functional group, ${\mathrm{\Lambda }}_i$, the infiltration rate of above-ground water into the soil, $I$, and the evaporation rate, $L$, are given by the expressions:

where $\overline{B}={\sum }_{j=1}^N{B}_j$ and $\overline{\overline{B}}={\sum }_{j=1}^N{Y}_j{B}_j$. We note that the biomass variable $B_i$ represents the above-ground biomass of all species with functional traits ${\displaystyle K_i,M_i}$ within the interval $\mathrm{\Delta }\chi$ that precedes ${\chi }_i$. In other words, while $B_i$ represents a biomass density in physical space, it does not represent a density in trait space, that is biomass per unit trait-space length. In this approach the solutions of the model Equations (1) depend on the choice of $N$. Specifically, the amplitude of the biomass distribution decreases as $N$ increases, although the metrics ${\chi }_{max},FR,FE$ are practically independent of $N$ for sufficiently large $N$ values. An alternative approach is to express the model equations in terms of the biomass densities $b_i=B_i/\mathrm{\Delta }\chi =N{B}_i$. Expressed in terms of ${\displaystyle b_i}$, also the biomass distributions are independent of $N$. That alternative model presentation, where $B_i$ is replaced by $b_i\mathrm{\Delta }\chi$, is also convenient when considering the continuum limit $N\to \mathrm{\infty }$; then $b_i\mathrm{\Delta }\chi$ becomes $b_{i}d\chi$ and the summation over the $N$ functional groups should simply be replaced by an integral.

A list of all model parameters, their descriptions, units and their values are given in Table 1. Below, we focus on the new elements in the model and on several aspects we wish to emphasize. The biomass dependence of the growth rates, ${\mathrm{\Lambda }}_i$, models competition for light and accounts for growth attenuation due to shading. That attenuation is described by the parameters $K_i$, which quantify the capacity of functional groups to capture light; high $K_i$ values represent plant species investing preferably in shoot growth that are less affected by shading. Note that the growth rate also includes attenuation due to self-shading (Noy-Meir, 1975). This attenuation form represents a relative decrease in the photosynthetic capacity of a plant, as photosynthesis is progressively limited to the upper layer of leaves. The parameter ${\mathrm{\Lambda }}_0$ represents the growth rate at low biomass values for which competition for light is absent.

The biomass dependence of the infiltration rate, $I$, is responsible for differential infiltration, quantified by the dimensionless parameter $0\le f\le 1$ and the parameter $Q$. Low $f$ values represent highly differential infiltration, low in bare soil ($I=fA$) and high ($Q$-dependent) in vegetation patches, and constitute an important element in the pattern-forming feedback associated with overland water flow toward denser vegetation patches (Meron, 2019). We assume that all species contribute to this differential infiltration, but these contributions are not necessarily equal as they depend on various factors, including the number density of individual plants and their root densities. We quantify these differences through the parameters $Y_i$.

The biomass dependence of the evaporation rate, $L$, accounts for reduced evaporation in vegetation patches due to shading, quantified by the parameter $R$. The parameter ${\displaystyle L_0}$ represents evaporate rate in bare soil. The reduced evaporation in vegetation patches is a positive feedback between biomass and water that can result in bistability of bare soil and uniform or patterned vegetation. This feedback is not scale dependent, as it does not involve water transport, and therefore cannot induce spatial patterning.

Tolerance to water stress is modeled through the mortality parameters $M_i$ – a plant investment in tolerating stress reduces the mortality rate. This parameter may lump together (as additional additive terms) other biomass-decay factors, such as grazing stress. For simplicity we choose $M_i$ to be independent of $W$ (unlike in Tzuk et al., 2019 for example). This choice can represent stress-tolerance mechanisms associated with plant architecture rather than phenotypic changes, such as hydraulically independent multiple stems that lead to a redundancy of independent conduits and higher resistance to drought (Schenk et al., 2008). We note that although $M_i$ is independent of $W$, lower values of $M_i$ do represent higher tolerance to water stress. This can be seen by considering the fitness of the $i$ th functional group ${\mathrm{\Lambda }}_{i}W-M_i$, which implies that functional groups with lower mortality tolerate lower soil-water content or higher water stress, as their fitness can remain positive.

Note that the atmosphere is considered in this model as the ecosystem’s environment, rather than a system part in feedback relationships with vegetation growth. As a consequence, it is quantified by a precipitation parameter $P$ in the equation for $H$ (1 c), assumed in this study to be constant. This simplification is justifiable in drylands where vegetation is relatively sparse and the overall transpiration is low. Also, for simplicity, we choose to describe here overland water flow as a linear diffusion problem by setting $\alpha =0$ in Equations (1). Although that process is nonlinear, the qualitative results and conclusions reported here do not depend on that choice (see Gilad et al., 2004 for the choice $\alpha =1$).

As pointed out earlier, we distinguish between different species through the different tradeoffs they make between shoot growth and tolerance to water stress. We capture this tradeoff using the parameters $K_i$ and $M_i$ through the tradeoff relations:

where $K_i=K({\chi }_i)$ and $M_i=M({\chi }_i)$. According to these relations, $\chi \to 0$ represents the functional group $(K_{max},M_{max})$ with highest investment in growth and lowest investment in tolerance (highest mortality), while $\chi =1$ represents the functional group $(K_{min},M_{min})$ with lowest investment in growth and highest investment in tolerance. This tradeoff is likely to affect the contributions, $Y_i$, of the various functional groups to the infiltration rate $I$; denser roots associated with lower-χ species (increased investment in shoot growth) make the soil more porous and increase the infiltration rate. An additional contribution to that effect is lower soil-crust coverage in patches of lower-χ species. We therefore assume the following form for $Y_i$:

The model for the single functional group $\chi =1$ used in deriving the results of Figure 4 is obtained from the community model (1) by setting all biomass variables identically to zero, apart of $B_N$ for which $\chi ={\chi }_N=N\mathrm{\Delta }\chi =1$.

We measure functional diversity using two metrics, functional richness, $FR$, and functional evenness, $FE$ (Mason et al., 2005). The first metric is given by the extent of the biomass distribution around ${\chi }_{max}$, as Figure 3 illustrates. The second metric contains information about the abundance of functional groups in the community and how even the abundance is among the groups. We use here the following analogue of the Shannon diversity index,

and the related index of Pielou for functional evenness Pielou, 1966,

We use numerical continuation methods (AUTO Doedel, 1981) to study the model equations for single functional groups, and numerical time-integration in the composed trait-space plane to study the full community model. We use periodic boundary conditions in $x$ and zero-flux conditions in χ. Initial conditions are chosen to contain all functional groups, even if at arbitrarily small biomass values. Such small values represent the presence of seeds that remain viable even when they cannot germinate (DeMalach et al., 2021). Source-code files for integrating the community and single functional-group models are available at the GitHub repository (Bera, 2021).

We wish to point out that the results presented here are not sensitive to the particular choice of the parameter values displayed in Table 1, and similar results have been obtained with other sets of parameter values. The main constraint on this choice is the need to capture the Turing instability of the uniform vegetation state, and the need to define the tradeoff relations such that the Turing threshold split the community into functional groups that form and do not form periodic patterns. Since vegetation patterns are observed on spatial scales that differ by orders of magnitude, from periodicity of tens of centimeters for herbaceous vegetation to periodicity of tens of meters for woody vegetation (Rietkerk et al., 2004), we focus on generic community aspects associated with vegetation patterning, rather than attempt to model a particular ecosystem with a specific spatial scale.

The current manuscript is a computational study, so no empirical data have been generated for this manuscript. Modelling code is uploaded to https://github.com/bidesh001/Plant-community-model (copy archived at https://archive.softwareheritage.org/swh:1:rev:976242ce55d8f501d25cb67c365004f04b895b6c).

1. 1. Angert AL
   2. Huxman TE
   3. Chesson P
   4. Venable DL

   (2009) Functional tradeoffs determine species coexistence via the storage effect

   PNAS 106:11641–11645.

   https://doi.org/10.1073/pnas.0904512106

   • PubMed
   • Google Scholar
2. 1. Bastiaansen R
   2. Jaïbi O
   3. Deblauwe V
   4. Eppinga MB
   5. Siteur K
   6. Siero E
   7. Mermoz S
   8. Bouvet A
   9. Doelman A
   10. Rietkerk M

   (2018) Multistability of model and real dryland ecosystems through spatial self-organization

   PNAS 115:11256–11261.

   https://doi.org/10.1073/pnas.1804771115

   • PubMed
   • Google Scholar
3. 1. Bastiaansen R
   2. Doelman A
   3. Eppinga MB
   4. Rietkerk M

   (2020) The effect of climate change on the resilience of ecosystems with adaptive spatial pattern formation

   Ecology Letters 23:414–429.

   https://doi.org/10.1111/ele.13449

   • PubMed
   • Google Scholar
4. 1. Bastiaansen R
   2. Doelman A

   (2019) The dynamics of disappearing pulses in a singularly perturbed reaction–diffusion system with parameters that vary in time and space

   Physica D: Nonlinear Phenomena 388:45–72.

   https://doi.org/10.1016/j.physd.2018.09.003

   • Google Scholar
5. Software

   1. Bera BK

   (2021) Python Codes for a Single Functional-Group Model and for the Plant Community Model

   Python Codes for a Single Functional-Group Model and for the Plant Community Model.

   https://github.com/bidesh001/Plant-community-model
6. 1. Bertness MD
   2. Callaway R

   (1994) Positive interactions in communities

   Trends in Ecology & Evolution 9:191–193.

   https://doi.org/10.1016/0169-5347(94)90088-4

   • PubMed
   • Google Scholar
7. 1. Bonanomi G
   2. Incerti G
   3. Stinca A
   4. Cartenì F
   5. Giannino F
   6. Mazzoleni S

   (2014) Ring formation in clonal plants

   Community Ecology 15:77–86.

   https://doi.org/10.1556/ComEc.15.2014.1.8

   • Google Scholar
8. 1. Borgogno F
   2. D'Odorico P
   3. Laio F
   4. Ridolfi L

   (2009) Mathematical models of vegetation pattern formation in ecohydrology

   Reviews of Geophysics 47:RG1005.

   https://doi.org/10.1029/2007RG000256

   • Google Scholar
9. 1. Cornacchia L
   2. van de Koppel J
   3. van der Wal D
   4. Wharton G
   5. Puijalon S
   6. Bouma TJ

   (2018) Landscapes of facilitation: how self-organized patchiness of aquatic macrophytes promotes diversity in streams

   Ecology 99:832–847.

   https://doi.org/10.1002/ecy.2177

   • PubMed
   • Google Scholar
10. 1. Dawes JH
    2. Williams JL

    (2016) Localised pattern formation in a model for dryland vegetation

    Journal of Mathematical Biology 73:63–90.

    https://doi.org/10.1007/s00285-015-0937-5

    • PubMed
    • Google Scholar
11. 1. De Deyn GB
    2. Raaijmakers CE
    3. Van der Putten WH

    (2004) Plant community development is affected by nutrients and soil biota

    Journal of Ecology 92:824–834.

    https://doi.org/10.1111/j.0022-0477.2004.00924.x

    • Google Scholar
12. 1. Deblauwe V
    2. Barbier N
    3. Couteron P
    4. Lejeune O
    5. Bogaert J

    (2008) The global biogeography of semi-arid periodic vegetation patterns

    Global Ecology and Biogeography 17:715–723.

    https://doi.org/10.1111/j.1466-8238.2008.00413.x

    • Google Scholar
13. 1. Deblauwe V
    2. Couteron P
    3. Lejeune O
    4. Bogaert J
    5. Barbier N

    (2010) Environmental modulation of self-organized periodic vegetation patterns in Sudan

    Ecography 33:1–13.

    https://doi.org/10.1111/j.1600-0587.2010.06694.x

    • Google Scholar
14. 1. Deblauwe V
    2. Couteron P
    3. Bogaert J
    4. Barbier N

    (2012) Determinants and dynamics of banded vegetation pattern migration in arid climates

    Ecological Monographs 82:3–21.

    https://doi.org/10.1890/11-0362.1

    • Google Scholar
15. 1. DeMalach N
    2. Kigel J
    3. Sternberg M

    (2021) The soil seed bank can buffer long‐term compositional changes in annual plant communities

    Journal of Ecology 109:1275–1283.

    https://doi.org/10.1111/1365-2745.13555

    • Google Scholar
16. 1. Dı́az S
    2. Cabido M

    (2001) Vive la différence: plant functional diversity matters to ecosystem processes

    Trends in Ecology & Evolution 16:646–655.

    https://doi.org/10.1016/S0169-5347(01)02283-2

    • Google Scholar
17. Conference

    1. Doedel E

    (1981)

    Auto: a program for the automatic bifurcation analysis of autonomous systems

    Congressus Numerantium. pp. 265–284.

    • Google Scholar
18. 1. Dovrat G
    2. Meron E
    3. Shachak M
    4. Golodets C
    5. Osem Y

    (2019) Plant size is related to biomass partitioning and stress resistance in water-limited annual plant communities

    Journal of Arid Environments 165:1–9.

    https://doi.org/10.1016/j.jaridenv.2019.04.006

    • Google Scholar
19. Report

    1. Duraiappah AK
    2. Naeem S

    (2005)

    Ecosystems and Human Well-Being: Biodiversity Synthesis

    Washington, United States: World Resources Institute.

    • Google Scholar
20. 1. Eppinga MB
    2. Rietkerk M
    3. Borren W
    4. Lapshina ED
    5. Bleuten W
    6. Wassen MJ

    (2008) Regular surface patterning of peatlands: confronting theory with field data

    Ecosystems 11:520–536.

    https://doi.org/10.1007/s10021-008-9138-z

    • Google Scholar
21. 1. Falik O
    2. Reides P
    3. Gersani M
    4. Novoplansky A

    (2003) Self/non-self discrimination in roots

    Journal of Ecology 91:525–531.

    https://doi.org/10.1046/j.1365-2745.2003.00795.x

    • Google Scholar
22. 1. Getzin S
    2. Yizhaq H
    3. Bell B
    4. Erickson TE
    5. Postle AC
    6. Katra I
    7. Tzuk O
    8. Zelnik YR
    9. Wiegand K
    10. Wiegand T
    11. Meron E

    (2016) Discovery of fairy circles in Australia supports self-organization theory

    PNAS 113:3551–3556.

    https://doi.org/10.1073/pnas.1522130113

    • PubMed
    • Google Scholar
23. 1. Gilad E
    2. von Hardenberg J
    3. Provenzale A
    4. Shachak M
    5. Meron E

    (2004) Ecosystem engineers: from pattern formation to habitat creation

    Physical Review Letters 93:098105.

    https://doi.org/10.1103/PhysRevLett.93.098105

    • PubMed
    • Google Scholar
24. 1. Gilad E
    2. Shachak M
    3. Meron E

    (2007) Dynamics and spatial organization of plant communities in water-limited systems

    Theoretical Population Biology 72:214–230.

    https://doi.org/10.1016/j.tpb.2007.05.002

    • PubMed
    • Google Scholar
25. 1. Gowda K
    2. Riecke H
    3. Silber M

    (2014) Transitions between patterned states in vegetation models for semiarid ecosystems

    Physical Review E 89:022701.

    https://doi.org/10.1103/PhysRevE.89.022701

    • Google Scholar
26. 1. Gratani L

    (2014) Plant phenotypic plasticity in response to environmental factors

    Advances in Botany 2014:1–17.

    https://doi.org/10.1155/2014/208747

    • Google Scholar
27. 1. Harrison S
    2. Spasojevic MJ
    3. Li D

    (2020) Climate and plant community diversity in space and time

    PNAS 117:4464–4470.

    https://doi.org/10.1073/pnas.1921724117

    • PubMed
    • Google Scholar
28. 1. Inderjit N
    2. Callaway RM
    3. Meron E

    (2021) Belowground feedbacks as drivers of spatial self-organization and community assembly

    Physics of Life Reviews 38:1–24.

    https://doi.org/10.1016/j.plrev.2021.07.002

    • PubMed
    • Google Scholar
29. 1. Inderjit N
    2. Duke SO

    (2003) Ecophysiological aspects of allelopathy

    Planta 217:529–539.

    https://doi.org/10.1007/s00425-003-1054-z

    • PubMed
    • Google Scholar
30. 1. Jaïbi O
    2. Doelman A
    3. Chirilus-Bruckner M
    4. Meron E

    (2020) The existence of localized vegetation patterns in a systematically reduced model for dryland vegetation

    Physica D: Nonlinear Phenomena 412:132637.

    https://doi.org/10.1016/j.physd.2020.132637

    • Google Scholar
31. 1. Knobloch E

    (2015) Spatial localization in dissipative systems

    Annual Review of Condensed Matter Physics 6:325–359.

    https://doi.org/10.1146/annurev-conmatphys-031214-014514

    • Google Scholar
32. 1. Lefever R
    2. Lejeune O

    (1997) On the origin of tiger bush

    Bulletin of Mathematical Biology 59:263–294.

    https://doi.org/10.1007/BF02462004

    • Google Scholar
33. 1. Lejeune O
    2. Tlidi M
    3. Lefever R

    (2004) Vegetation spots and stripes: dissipative structures in arid landscapes

    International Journal of Quantum Chemistry 98:261–271.

    https://doi.org/10.1002/qua.10878

    • Google Scholar
34. 1. Marasco A
    2. Iuorio A
    3. Cartení F
    4. Bonanomi G
    5. Tartakovsky DM
    6. Mazzoleni S
    7. Giannino F

    (2014) Vegetation pattern formation due to interactions between water availability and toxicity in plant-soil feedback

    Bulletin of Mathematical Biology 76:2866–2883.

    https://doi.org/10.1007/s11538-014-0036-6

    • PubMed
    • Google Scholar
35. 1. Mason NWH
    2. Mouillot D
    3. Lee WG
    4. Wilson JB

    (2005) Functional richness, functional evenness and functional divergence: the primary components of functional diversity

    Oikos 111:112–118.

    https://doi.org/10.1111/j.0030-1299.2005.13886.x

    • Google Scholar
36. Book

    1. Meron E

    (2015) Nonlinear Physics of Ecosystems

    Florida, United States: CRC Press.

    https://doi.org/10.1201/b18360

    • Google Scholar
37. 1. Meron E

    (2016) Pattern formation--A missing link in the study of ecosystem response to environmental changes

    Mathematical Biosciences 271:1–18.

    https://doi.org/10.1016/j.mbs.2015.10.015

    • PubMed
    • Google Scholar
38. 1. Meron E

    (2018) From patterns to function in living systems: dryland ecosystems as a case study

    Annual Review of Condensed Matter Physics 9:79–103.

    https://doi.org/10.1146/annurev-conmatphys-033117-053959

    • Google Scholar
39. 1. Meron E

    (2019) Vegetation pattern formation: the mechanisms behind the forms

    Physics Today 72:30–36.

    https://doi.org/10.1063/PT.3.4340

    • Google Scholar
40. 1. Meron E
    2. Bennett JJR
    3. Fernandez-Oto C
    4. Tzuk O
    5. Zelnik YR
    6. Grafi G

    (2019) Continuum modeling of discrete plant communities: why does it work and why is it advantageous?

    Mathematics 7:987.

    https://doi.org/10.3390/math7100987

    • Google Scholar
41. 1. Nathan J
    2. Osem Y
    3. Shachak M
    4. Meron E

    (2016) Linking functional diversity to resource availability and disturbance: a mechanistic approach for water-limited plant communities

    Journal of Ecology 104:419–429.

    https://doi.org/10.1111/1365-2745.12525

    • Google Scholar
42. 1. Noy-Meir I

    (1975) Stability of grazing systems: an application of Predator-Prey graphs

    The Journal of Ecology 63:459–481.

    https://doi.org/10.2307/2258730

    • Google Scholar
43. 1. O'Sullivan JD
    2. Knell RJ
    3. Rossberg AG

    (2019) Metacommunity-scale biodiversity regulation and the self-organised emergence of macroecological patterns

    Ecology Letters 22:1428–1438.

    https://doi.org/10.1111/ele.13294

    • PubMed
    • Google Scholar
44. 1. Parra-Rivas P
    2. Fernandez-Oto C

    (2020) Formation of localized states in dryland vegetation: bifurcation structure and stability

    Physical Review. E 101:052214.

    https://doi.org/10.1103/PhysRevE.101.052214

    • PubMed
    • Google Scholar
45. 1. Paz-Kagan T
    2. DeMalach N
    3. Zaady E
    4. Shachak M

    (2019) Resource redistribution effects on annual plant communities in a runoff harvesting system in dryland

    Journal of Arid Environments 171:103984.

    https://doi.org/10.1016/j.jaridenv.2019.05.012

    • Google Scholar
46. 1. Pérez-Ramos IM
    2. Matías L
    3. Gómez-Aparicio L
    4. Godoy Ó

    (2019) Functional traits and phenotypic plasticity modulate species coexistence across contrasting climatic conditions

    Nature Communications 10:2555.

    https://doi.org/10.1038/s41467-019-10453-0

    • PubMed
    • Google Scholar
47. 1. Pescador DS
    2. Iriondo JM
    3. Losapio G
    4. Escudero A

    (2020) The assembly of plant–patch networks in Mediterranean alpine grasslands

    Journal of Plant Ecology 13:273–280.

    https://doi.org/10.1093/jpe/rtaa011

    • Google Scholar
48. 1. Pielou EC

    (1966) The measurement of diversity in different types of biological collections

    Journal of Theoretical Biology 13:131–144.

    https://doi.org/10.1016/0022-5193(66)90013-0

    • Google Scholar
49. 1. Pomeau Y

    (1986) Front motion, metastability and subcritical bifurcations in hydrodynamics

    Physica D: Nonlinear Phenomena 23:3–11.

    https://doi.org/10.1016/0167-2789(86)90104-1

    • Google Scholar
50. 1. Pringle RM
    2. Tarnita CE

    (2017) Spatial Self-Organization of ecosystems: integrating multiple mechanisms of Regular-Pattern formation

    Annual Review of Entomology 62:359–377.

    https://doi.org/10.1146/annurev-ento-031616-035413

    • PubMed
    • Google Scholar
51. 1. Rietkerk M
    2. Boerlijst MC
    3. van Langevelde F
    4. Hillerislambers R
    5. de Koppel J
    6. Kumar L
    7. Prins HH
    8. de Roos AM

    (2002) Self-organization of vegetation in arid ecosystems

    The American Naturalist 160:524–530.

    https://doi.org/10.1086/342078

    • PubMed
    • Google Scholar
52. 1. Rietkerk M
    2. Dekker SC
    3. de Ruiter PC
    4. van de Koppel J

    (2004) Self-organized patchiness and catastrophic shifts in ecosystems

    Science 305:1926–1929.

    https://doi.org/10.1126/science.1101867

    • PubMed
    • Google Scholar
53. 1. Rietkerk M
    2. van de Koppel J

    (2008) Regular pattern formation in real ecosystems

    Trends in Ecology & Evolution 23:169–175.

    https://doi.org/10.1016/j.tree.2007.10.013

    • PubMed
    • Google Scholar
54. 1. Ruiz-Reynés D
    2. Gomila D
    3. Sintes T
    4. Hernández-García E
    5. Marbà N
    6. Duarte CM

    (2017) Fairy circle landscapes under the sea

    Science Advances 3:e1603262.

    https://doi.org/10.1126/sciadv.1603262

    • PubMed
    • Google Scholar
55. 1. Schenk HJ
    2. Espino S
    3. Goedhart CM
    4. Nordenstahl M
    5. Cabrera HI
    6. Jones CS

    (2008) Hydraulic integration and shrub growth form linked across continental aridity gradients

    PNAS 105:11248–11253.

    https://doi.org/10.1073/pnas.0804294105

    • PubMed
    • Google Scholar
56. 1. Sherratt JA

    (2013) Pattern solutions of the klausmeier model for banded vegetation in semiarid environments V: the transition from patterns to desert

    SIAM Journal on Applied Mathematics 73:1347–1367.

    https://doi.org/10.1137/120899510

    • Google Scholar
57. 1. Siteur K
    2. Siero E
    3. Eppinga MB
    4. Rademacher JDM
    5. Doelman A
    6. Rietkerk M

    (2014) Beyond turing: the response of patterned ecosystems to environmental change

    Ecological Complexity 20:81–96.

    https://doi.org/10.1016/j.ecocom.2014.09.002

    • Google Scholar
58. 1. Suding KN
    2. Lavorel S
    3. Chapin FS
    4. Cornelissen JHC
    5. Diaz S
    6. Garnier E
    7. Goldberg D
    8. Hooper DU
    9. Jackson ST
    10. Navas M-L

    (2008) Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants

    Global Change Biology 14:1125–1140.

    https://doi.org/10.1111/j.1365-2486.2008.01557.x

    • Google Scholar
59. Book

    1. Tilman D

    (1982)

    monographs in population biology

    In: Tilman D, editors. Resource Competition and Community Structure. New Jersey, United states: Princeton University Press. pp. 1–12.

    • Google Scholar
60. Book

    1. Tlidi M
    2. Lefever R
    3. Vladimirov A

    (2008) Vegetation Clustering, Localized Bare Soil Spots and Fairy Circles

    In: Lefever R, editors. Dissipative Solitons: From Optics to Biology and Medicine. Berlin, Germany: Springer. pp. 11–20.

    https://doi.org/10.1007/978-3-540-78217-9

    • Google Scholar
61. 1. Tzuk O
    2. Ujjwal SR
    3. Fernandez-Oto C
    4. Seifan M
    5. Meron E

    (2019) Period doubling as an Indicator for ecosystem sensitivity to climate extremes

    Scientific Reports 9:19577.

    https://doi.org/10.1038/s41598-019-56080-z

    • PubMed
    • Google Scholar
62. 1. Valentin C
    2. d'Herbès JM
    3. Poesen J

    (1999) Soil and water components of banded vegetation patterns

    Catena 37:1–24.

    https://doi.org/10.1016/S0341-8162(99)00053-3

    • Google Scholar
63. 1. van der Putten WH
    2. Bardgett RD
    3. Bever JD
    4. Bezemer TM
    5. Casper BB
    6. Fukami T
    7. Kardol P
    8. Klironomos JN
    9. Kulmatiski A
    10. Schweitzer JA
    11. Suding KN
    12. Van de Voorde TFJ
    13. Wardle DA

    (2013) Plant-soil feedbacks: the past, the present and future challenges

    Journal of Ecology 101:265–276.

    https://doi.org/10.1111/1365-2745.12054

    • Google Scholar
64. 1. Vandermeer J
    2. Yitbarek S

    (2012) Linking spatial self-organization to community assembly and biodiversity

    Journal of Theoretical Biology 300:48–56.

    https://doi.org/10.1016/j.jtbi.2012.01.005

    • Google Scholar
65. 1. von Hardenberg J
    2. Meron E
    3. Shachak M
    4. Zarmi Y

    (2001) Diversity of vegitation patterns and desertification

    Physical Review Letters 89:198101.

    https://doi.org/10.1103/PhysRevLett.87.198101

    • Google Scholar
66. 1. Weiher E
    2. Freund D
    3. Bunton T
    4. Stefanski A
    5. Lee T
    6. Bentivenga S

    (2011) Advances, challenges and a developing synthesis of ecological community assembly theory

    Philosophical Transactions of the Royal Society B: Biological Sciences 366:2403–2413.

    https://doi.org/10.1098/rstb.2011.0056

    • Google Scholar
67. 1. Weltzin JF
    2. Pastor J
    3. Harth C
    4. Bridgham SD
    5. Updegraff K
    6. Chapin CT

    (2000) Response of bog and fen plant communities to warming and water-table manipulations

    Ecology 81:3464–3478.

    https://doi.org/10.1890/0012-9658(2000)081[3464:ROBAFP]2.0.CO;2

    • Google Scholar
68. 1. Yizhaq H
    2. Gilad E
    3. Meron E

    (2005) Banded vegetation: biological productivity and resilience

    Physica A: Statistical Mechanics and Its Applications 356:139–144.

    https://doi.org/10.1016/j.physa.2005.05.026

    • Google Scholar
69. 1. Yizhaq H
    2. Shachak M
    3. Meron E

    (2020) A model study of terraced riverbeds as novel ecosystems

    Scientific Reports 10:3782.

    https://doi.org/10.1038/s41598-020-60706-y

    • PubMed
    • Google Scholar
70. 1. Zelnik YR
    2. Kinast S
    3. Yizhaq H
    4. Bel G
    5. Meron E

    (2013) Regime shifts in models of dryland vegetation

    Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371:20120358.

    https://doi.org/10.1098/rsta.2012.0358

    • Google Scholar
71. 1. Zelnik YR
    2. Meron E
    3. Bel G

    (2015) Gradual regime shifts in fairy circles

    PNAS 112:12327–12331.

    https://doi.org/10.1073/pnas.1504289112

    • PubMed
    • Google Scholar
72. 1. Zhao L-X
    2. Xu C
    3. Ge Z-M
    4. van de Koppel J
    5. Liu Q-X

    (2019) The shaping role of self-organization: linking vegetation patterning, plant traits and ecosystem functioning

    Proceedings of the Royal Society B: Biological Sciences 286:20182859.

    https://doi.org/10.1098/rspb.2018.2859

    • Google Scholar
73. 1. Zhao L-X
    2. Zhang K
    3. Siteur K
    4. Li X-Z
    5. Liu Q-X
    6. van de Koppel J

    (2021) Fairy circles reveal the resilience of self-organized salt marshes

    Science Advances 7:abe1100.

    https://doi.org/10.1126/sciadv.abe1100

    • Google Scholar

Acceptance summary:

In this paper, the authors use a trait-based model of plant growth and water flow in drylands to show that under increasing water shortage, spatial self-organization can help plant communities to maintain biodiversity and thus ecosystem functioning. Spatially heterogeneous ecosystem management may support these processes.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting the paper "Linking spatial self-organization to community assembly and biodiversity" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Gianalberto Losapio (Reviewer #2); Mara Baudena (Reviewer #3).

While we cannot go forward with the current version of the manuscript, since revision would almost certainly take substantially longer than allowed per eLife policy, we remain very interested in the work. We would therefore welcome a substantially revised version, which we would treat as a new submission, but aim to recruit the same editors and reviewers.

Reviewer #1:

Bera et al. study the response of vegetation in water-limited ecosystems to changes in the precipitation regime. Previous studies have shown that spatial processes, in particular the redistribution of (soil and surface) water, may play an important role in mediating the ecosystem response. An important consequence of this redistribution is the spatial self-organization of vegetation into regular spatial patterns, consisting of vegetation patches that act as sinks for (surface) water, and surrounding areas of bare soil that act as water sources. At the ecosystem level, the additional water input in vegetation patches may enable vegetation to persist at precipitation levels that are too low to sustain a spatially uniform cover.

While most model studies of spatial self-organization and pattern formation describe vegetation dynamics through 1-2 biomass variables, the current study extends this previous work by considering a trait diversity gradient, considering a large number (N=128) of discrete trait classes that range from stress-tolerant to fast-growing characteristics. The results show that in the absence of spatial pattern formation, a decrease in precipitation leads to a shift in the biomass distribution toward the more stress-tolerant trait classes. At the onset of pattern formation, however, soil water availability increases at the locations where vegetation patches form, enabling the more fast-growing trait classes to increase in biomass, and this shift is accompanied by an increase in functional diversity of trait classes as well. It is also shown that once these patterned ecosystem states are formed, the main adaptation to further decreases in precipitation occurs either through shrinking the size of existing patches, or by reducing the number of patches; in contrast, biomass and community composition of the patches remains relatively stable. Finally, it is shown that for certain precipitation conditions, functional diversity is maximized when the ecosystem is in a hybrid state, where part of the landscape has a spatially uniform vegetation cover, and part of the landscape is in a patterned state.

A potential strength of this paper is that the community assembly and biodiversity perspective on spatial self-organization may highlight the relevance of pattern formation in ecosystems more clearly to a broad audience. The formulation of a trait/strategy gradient of discrete classes is certainly an interesting suggestion to connect the typical single/few biomass variable(s) approach to a community-level approach. The community assembly process is modelled in a very specific way, and the manuscript would benefit from an expanded ecological motivation of the processes that are being mimicked, and thereby explain more clearly what taxonomic level of organization is being considered. In addition, it would be useful if the authors could provide further clarification as to what extent the community diversity dynamics can be separated from total biomass dynamics of patterned water-limited ecosystems given the current approach. These points are explained in further detail below.

• First, it was not entirely clear to this reviewer how the reaction parts of the model equations determine the optimal trait value χ, and how this value varies as a function of precipitation. Assuming a single trait class, and plotting the relevant equilibrium values of the three state variables shed some light on this issue. [Unfortunately, there does not seem to be a possibility to attach the figure with these plots to this review report]. Assuming the non-spatial equilibrium solution was derived correctly , the optimum biomass (B) value shifts across the trait spectrum with changing precipitation (in the non-spatial model version, solving the surface water equation for equilibrium will always yield that all precipitation infiltrates, i.e. regardless of the values of surface water, H, and χ). The equilibrium of soil water availability (W), which is the growth limiting resource of the vegetation, shows an inverse pattern with biomass. This result is in line with a classical results (e.g. Tilman 1982), in that the most successful strategy is the one that is able to reduce the limiting resource to the lowest equilibrium value. With all trait classes competing for the soil water resource, however, it is then not immediately clear why the most successful trait class is not outcompeting the other classes. This leads to a second point, about the way in which community trait adaptation is modelled.

• The authors model trait adaptation through a diffusion approximation between trait classes. That is, every timestep, a small amount of biomass flows from the class with higher biomass to the neighboring trait class with lower biomass. From an ecological point of view, it seems that this process is describing adaptation of vegetation that is already present, so this process seems to be limited to intraspecific phenotypic plasticity. From the text, however, it seems that the trait classes correspond to higher taxonomic levels of organization, when describing shifts from fast growing to stress-tolerant species, for example. It is not entirely clear, however, how biomass flows as assumed in the model could occur at these higher levels of organization.

• Combining the observations from the previous two points, there is a concern that for a given level of precipitation, there is a single trait class with optimal biomass/lowest soil water level that is dominant, with the neighboring trait classes being sustained by the diffusion of biomass from the optimal class to neighboring inferior classes. This would seem a bit problematic, as it would mean that most classes are not a true fit for the environment, and only persist due to the continuous inflow of biomass. Taking a clue from the previous papers of the authors, it seems this may not be the case, though. Specifically, in the paper by Nathan et al. (2016) it seems that all trait classes are started at low initial biomass density, and the resulting steady state (in the absence of biomass flows between classes) seems to show similar biomass profiles as shown in Figures 4,5 and 7 of the current paper. While the current model formulation seems slightly different, similar results may apply here. Indeed, keeping all trait classes at non-zero (but low) density, and when the (abiotic and biotic) environment permits, let each class increase in biomass seems like the most straightforward approach to model community assembly dynamics. Given the above discussion about these trait classes competing for a single resource (soil water), and one trait class being able to drive this resource availability to the lowest level, it would then be useful to readers to explain why multiple trait classes can coexist here, and how (for spatial uniform solutions) the equilibrium soil water level with multiple trait classes present compares to the equilibrium soil water level when only the optimal trait class is present. Furthermore, if results as presented in Nathan et al. (2016) indeed hold in the current case, perhaps it means that the biomass profile responses as shown in e.g. Figure 5 would also occur if there was no biomass flow between trait classes included, but that the time needed to adjust the profile would take much longer as compared to when the drift term/second trait derivative is included. In summary, further clarification of what the biomass flows between classes represent, and the role it plays in driving the presented results would be useful for readers.

• In addition, it would be useful for readers to understand to what extent the shifts in average trait values and functional diversity can be decoupled from the biomass and soil water responses to changes in precipitation that would occur in a model with only a single biomass variable. For example, early studies on self-organization in semi-arid ecosystems already showed that the shift toward a patterned state involved the formation of patches with higher biomass, and higher soil water availability, as compared to the preceding spatially uniform state, and that the biomass in these patches remains relatively stable under decreasing rainfall, while their geometry changes (e.g. Rietkerk et al. 2002). It has also been observed that for a given environmental condition, biomass in vegetation patches tends to increase with pattern wavelength (e.g. Bastiaansen and Doelman 2018; Bastiaansen et al. 2018). Given the model formulation, one wonders whether higher biomass in the single variable model is not automatically corresponding to higher abundance of faster growing species and a higher functional diversity (as the diffusion of biomass can cover a broader range when starting from higher mass in the optimal trait class). There are some indications in the current work that the linkage is more complicated, for example, the biomass peak in Figure 7c is lower, but also broader as compared to the distribution of Figure 7b, but it is currently not entirely clear how this result can be explained (for example, it might be the case that in the spatially patterned states, the biomass profiles also vary in space).

• The possibility of hybrid states, where part of the landscape is in a spatially uniform state, while the other part of the landscape is in a patterned state, is quite interesting. To better understand how such states could be leveraged in management strategies, it would be useful if a bit more information could be provided on how these hybrid states emerge, and whether one can anticipate whether a perturbation will grow until a fully patterned state, or whether the expansion will halt at some point, yielding the hybrid state. It seems that being able to distinguish these case would be necessary in the design of planning and management strategies. Also, in Figure 3a, the region of parameter space in which hybrid states occur is not very large; it is not entirely clear whether the full range of hybrid states is left out here for visual considerations, or whether these states only occur within this narrow range in the vicinity of the Turing instability point.

References:

Bastiaansen R, Doelman A. 2018. The dynamics of disappearing pulses in a singularly perturbed reaction-diffusion system with parameters that vary in time and space. Physica D 388: 45-72.

Bastiaansen R, Jaïbi O, Deblauwe V, Eppinga MB, Siteur K, Siero E, Mermoz S, Bouvet A, Doelman A, Rietkerk M. 2018. Multistability of model and real dryland ecosystems through spatial self-organization. Proceedings of the National Academy of Sciences USA 115:11256-11261.

Nathan J, Osem Y, Shachak M, Meron E. 2016. Linking functional diversity to resource availability and disturbance: a mechanistic approach for water limited plant communities. Journal of Ecology 104: 419-429.

Rietkerk M, Boerlijst MC, van Langevelde F, HilleRisLambers R, van de Koppel J, Kumar L, Prins HHT, De Roos AM. 2002. Self-organization of vegetation in arid ecosystems. American Naturalist 160: 524-530.

Tilman D. 1982. Resource competition and community structure. Princeton University Press, Princeton, NJ, USA.

Comments for the authors:

• Line 17: the term "re-patterning" may read as a non-patterned state becoming patterned again, whereas here it seems to refer to a spatial rearrangement of an existing patterned state.

• Line 39: resources (i.e. plural)?

• Lines 80-99: This is an introduction to the model description, rather than a result, as the header suggests.

• Lines 100-164: This is the model description, which seems to be part of the material and methods rather than a result.

• Line 179: when χ increase from 0.95 to 1.00 however, it seems that the Turing threshold start to increase, how can this reversal be explained?

• Lines 302-310: this explanation is clear, but it is an example that can also be explained by the biomass dynamics of a single variable model.

• Line 329: this is case where it would be useful for readers to understand how one can anticipate the formation of either hybrid or fully patterned states, and how this relates to the particular perturbation(s) imposed.

• Figures: Why are the biomass values in Figure 4,5 and 7 about an order of magnitude higher than in Figure 3?

Reviewer #2:

Conspicuous, repetitive patterns such as spots and stripes can be observed in every biome throughout the world. This work provides a new theoretical model for understanding self-organization of vegetation patterns in arid ecosystems and their response to climate (precipitation) change. Processes of spatial self-organization underlying the development of vegetation patterns have been studied for decades, with roots in the work of the great scientist Alan Turing. Ecologists use the Turing reaction-diffusion theory that builds on positive feedback relations between two variables, namely vegetation growth and water transport. Yet, it has been difficult to include multiple, different species as in real-world vegetations.

This paper addresses such shortcoming and extends previous vegetation pattern formation models by including different plant types. It provides a general framework that builds on the resource allocation tradeoff between growth versus stress-tolerance. Authors show when and how vegetation is robust to changes in precipitation via spatial self-organization and selection (differential plant mortality) along the growth-tolerance tradeoff. With increasing aridity, the ecosystem shifts from spatial uniform vegetation to patterned one, such as stripes, and, with further drought, to bare ground. Notably, self-organizing processes mitigate the impact of drought on ecosystem functioning and services by allowing fast-growing, productive species to persist in drier climate. This framework and associated results have important implications for the conservation and management of arid ecosystems and rangelands.

The conclusions of this paper are mostly well supported by data, but some aspects of model presentation, parameter choices, and data interpretation need to be clarified and extended.

1) Model presentation. It would be better to explain the model in ecological terms first, clarifying parameter biological meaning and justifying their choice. In doing so, creating a specific 'Methods' section, which now is lacking, would be of help too. Authors should clarify whether and how the model follows the conservation of mass principle involving precipitation and evapotranspiration. Are root growth and seed dispersal included for this purpose? Why are they not referred to any further in the analysis and discussion? Why a specific term for plant transpiration is not included, or is to somehow phenomenologically incorporated into the growth-tolerance tradeoff? In doing so, authors should also pay attention to water balance as above (H) and below (W) ground water are not independent from each other.

Another unclear point is that growth rates for the same plant functional groups are assumed to be constant among different species within the same group and are confounded by biomass production. Why is that the case? Furthermore, how many different species are characterizing each functional group? How are interspecific interactions accounted for (more specifically, see comment below)?

Finally, stress tolerance is purely phenomenological. There is no actual mechanism/parameter describing it. Rather, it "simply" appears as low/high mortality, which in turn is said to be due to high/low tolerance. This leads to a sort of circularity between mortality and tolerance. Yet, mortality can occur due to other biophysical factors (e.g. disturbance, fire, herbivory, pathogens). A drawback of this assumption is that a mechanism of drought tolerance is often to invest in belowground organs, including roots. However, according to the proposed model, it turns out that fast growing species with low investment in tolerance also have high investment in roots; vice versa, tolerant species have low investment in roots. This is a bit counterintuitive and not well biologically supported.

2) Parameter choice.

N = 128 is an extremely high number for plant functional groups. It is even quite unrealistic to have 128 species per square meter, so this value is not very reasonable. Please run the model and report results with more realistic N (e.g. from 4-64) as well as with different sets of N values keeping all other parameters constant.

Gamma (rate of water uptake by plants' roots): why is it in that unit of m²/kg * y? Why are you now considering the area (and not the volume) per biomass unit?

A is not defined in the text.

M min: why 0.5 mortality? Having M max set to 0.9, please consider a lower mortality value set to 0.1, and please report evidence (hopefully) demonstrating the robustness of results to such change.

K_min and K_max are in two different units, and should both be kg/m².

Values of precipitation (P, mean annual precipitation) are not reported.

3) Results presentation and interpretation.

Parameter range of precipitation in figure 3 is odd. Why in one case precipitation ranges from 0 to 160 while in another it is only 60-120? Furthermore, in paragraph 198-213 and associated results in Figure 5. the Choice of precipitation values is somehow discordant from the previous model. Please provide motivation for this choice, clarify and uniformize it.

Throughout the text, authors claim to address plant-plant interactions, particularly intra and interspecific competition. However, it is not clear how competition was modelled neither whether it was included in the model. In its current state, it is just an assumption pulled out when discussing results – a classic 'passepartout' used by ecologists. Furthermore, why only competition is invoked in interpreting results when facilitation is known to be much more relevant in pattern formation and biodiversity maintenance in arid systems?

Finally, authors seem to create confusion around community composition, which is defined as the (taxonomic) identity of all different species inhabiting a community. Notably, it is remarkably different from the x_max parameter used in the model, which as a matter of facts is just the value of the most productive (notably, not necessarily the most abundant) functional group.

Reviewer #3:

In this paper, the authors use a mathematical model of plant and water dynamics in drylands to show that drylands adaptive capacity to respond to changes, via spatial self-organization in space has also beneficial effects in preserving its biodiversity and ecosystem functions.

The model is an extension of a large body of previous, well-established works on plant self-organisation in drylands. The model is well described and motivated (with one main exception, see below), the analyses are robust and the results are very convincingly supporting the conclusions. I however have an issue with one of the assumptions in the model equations. The authors included a term for "mutations" in traits that (1) is not introduced or motivated (2) its effects/importance are not highlighted by specific analyses (3) the possible implications or limitations connected to it are not discussed. To my knowledge, this term is also not based on earlier work. All these elements need to be included, as at the moment is for example unclear what the authors intended to represent by including the mutation term (evolutionary time scales? Or adaptation?). Also, it would be especially good to include an analysis of how influential this term is for the final results.

Assuming the authors can address this one concern, the results are surely important as they connect for the first time plant spatial self-organization to its biodiversity preservation, in the face of future expected climatic changes and probable land degradation. These findings, although theoretical, have the potential to be useful also for guiding adaptive and dynamic land management, as they underline the importance of taking into account spatial vegetation distribution in drylands management.

Besides the major point about the mutation term, I list here two other important points:

– The authors state that they represent highly tolerant plants by representing the plants with a small mortality. However, in their model, plant mortality does not depend on soil water levels. How can the authors reconcile these two aspects? Also, one could argue that mortality is related to the average life span, not specifically to tolerance to highly stressful condition. The authors should better justify this point and discuss the implication of this assumption.

– In the model, there is shading feedback too, not only infiltration feedback. However the authors state there's only infiltration feedback in l. 84, could they please explain?

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Reviewer #1:

[…]

A potential strength of this paper is that the community assembly and biodiversity perspective on spatial self-organization may highlight the relevance of pattern formation in ecosystems more clearly to a broad audience. The formulation of a trait/strategy gradient of discrete classes is certainly an interesting suggestion to connect the typical single/few biomass variable(s) approach to a community-level approach. The community assembly process is modelled in a very specific way, and the manuscript would benefit from an expanded ecological motivation of the processes that are being mimicked, and thereby explain more clearly what taxonomic level of organization is being considered.

We follow the more recent trait-based approach that shifts the focus from species (and the many traits by which they differ from one another) to groups of species that share the same values of selected functional traits. Since the general context is ecosystem response to drier climates, we choose the functional traits to include a response trait associated with stress tolerance and an effect trait associated with biomass production. We further assume a tradeoff between the two traits which is well supported by earlier studies (see e.g. Angert et al. 2009, https://doi.org/10.1073/pnas.0904512106). So, indeed, the choice we make in characterizing the community is quite specific, but it is highly relevant to the ecological context considered of dryland plant communities where plants compete primarily for water and light. The taxonomic level we consider is species except that we group them in a manner that is more transparent to questions of ecosystem function, ignoring differences between species that are not significant to these questions.

We expanded considerably the text in the section “Modeling spatial assembly of dryland plant communities” to clarify the ecological motivation of the processes we model.

In addition, it would be useful if the authors could provide further clarification as to what extent the community diversity dynamics can be separated from total biomass dynamics of patterned water-limited ecosystems given the current approach. These points are explained in further detail below.

The model describes the dynamics of all functional groups, which provides the biomass distribution 𝐵 = 𝐵(𝜒) in trait space (in the case of patterned states we first integrate over space). That distribution contains information about various community-level properties, including functional diversity (richness, evenness) as figure 3 in the revised manuscript illustrates, and total biomass, which is the area below the distribution curve. The two types of dynamics are tightly connected and cannot be separated, but in principle the approach can be used to study the relationships between diversity and total biomass by calculating biomass distributions along the rainfall gradient and extracting the two properties from the distributions.

We added in the section “Modeling spatial assembly of dryland plant communities” the information that the biomass distribution also contains information about the total biomass.

• First, it was not entirely clear to this reviewer how the reaction parts of the model equations determine the optimal trait value χ, and how this value varies as a function of precipitation.

The ‘optimal’ trait value 𝜒_𝑚𝑎𝑥 is determined by the interspecific interactions that the model captures, which divide into ‘direct’ and ‘indirect’ interactions. The direct interactions are captured by the dependence of the growth rate Λ_𝑖 of the ith functional group (see Equation 1a) on the aboveground biomass values of all functional groups, Λ_𝑖 = Λ_𝑖(𝐵₁,𝐵₂,… , 𝐵_𝑁) (see Equation 2). This dependence represents competition for light (taller plants are better competitors) and includes the effect of self-shading. The indirect interactions are through the water uptake term in the soil-water Equation 1b (2^nd term from right) and the water dependence of the biomass growth term in Equation 1a. These terms represent competition for water. For a given precipitation value 𝑃 the net effect of these interspecific interactions result in a particular functional group 𝜒_𝑚𝑎𝑥 which is most abundant. For spatially uniform vegetation, as 𝑃 is increased 𝜒_𝑚𝑎𝑥 moves to lower values. The precipitation increases surface water (Equation 1c) and consequently the amount of water 𝐼𝐻 infiltrating into the soil. The increased soil water gives competitive advantage to species investing in growth, mainly because they better compete for light as they grow taller, and therefore 𝜒_𝑚𝑎𝑥 decreases.

Assuming a single trait class, and plotting the relevant equilibrium values of the three state variables shed some light on this issue. [Unfortunately, there does not seem to be a possibility to attach the figure with these plots to this review report]. Assuming the non-spatial equilibrium solution was derived correctly , the optimum biomass (B) value shifts across the trait spectrum with changing precipitation (in the non-spatial model version, solving the surface water equation for equilibrium will always yield that all precipitation infiltrates, i.e. regardless of the values of surface water, H, and χ). The equilibrium of soil water availability (W), which is the growth limiting resource of the vegetation, shows an inverse pattern with biomass. This result is in line with a classical results (e.g. Tilman 1982), in that the most successful strategy is the one that is able to reduce the limiting resource to the lowest equilibrium value. With all trait classes competing for the soil water resource, however, it is then not immediately clear why the most successful trait class is not outcompeting the other classes. This leads to a second point, about the way in which community trait adaptation is modelled.

With the current model and parameters set the most successful trait does eventually outcompete all other traits, when trait diffusion is set to zero, 𝐷_𝜒 = 0. This is, however, a very long process because the most successful trait suffers from self-shading at late growth stages, which slows down its growth and allows nearby traits to survive for a long time. Choosing a finite but very small 𝐷_𝜒 values that represent mutations occurring on evolutionarily long times counteracts the exclusion process and results in a stationary asymptotic community, as Figure 3 in the revised manuscript shows (this behavior is reminiscent of optical solitons, where self-focusing instability is balanced by dispersion). We note that modeling stronger growth-inhibiting factors, such as pathogens, by including a factor of the form (1 − 𝐵_𝑖/𝐾) to the growth rate, results in an asymptotic stationary community also for 𝐷_𝜒 = 0 (see also earlier studies Nathan et al. 2016, Yizhaq et al. 2020).

We revised original Figure 4 (now Figure 3) by adding a new part (Figure 3a) that shows the exclusion process for 𝐷_𝜒 = 0, and the effect of the counter-acting process of trait diffusion, which results in an asymptotic distribution of finite width (Figure 3b) from which community level properties such as functional diversity can be derived. We also extended the text in section “Modeling spatial assembly of dryland plant communities” (last paragraph) to clarify the two counter-acting processes of exclusion because of interspecific competition for water and light, and trait diffusion driven by mutations, which together culminate in an asymptotic biomass distribution along the 𝜒 axis of finite width.

• The authors model trait adaptation through a diffusion approximation between trait classes. That is, every timestep, a small amount of biomass flows from the class with higher biomass to the neighboring trait class with lower biomass. From an ecological point of view, it seems that this process is describing adaptation of vegetation that is already present, so this process seems to be limited to intraspecific phenotypic plasticity. From the text, however, it seems that the trait classes correspond to higher taxonomic levels of organization, when describing shifts from fast growing to stress-tolerant species, for example. It is not entirely clear, however, how biomass flows as assumed in the model could occur at these higher levels of organization.

We do not study in this work adaptation through diffusion in trait space. That kind of adaptive dynamics can indeed be studied with the current model, but with different initial conditions, namely, initial conditions corresponding to a single resident trait where the biomass of all other traits is zero. The resulting dynamics of mutations and succession are then very slow, occurring on evolutionarily long time scales set by the small value of 𝐷_𝜒 (e.g. 10⁻⁶). In this study the initial conditions represent the presence of all traits, even if at very low biomass values that may represent a pool of seeds that germinate once environmental conditions allow. For a given precipitation value 𝑃, the functional traits we consider determine which functional groups (of species) overcome environmental filtering and grow, and which of the growing traits survive the competition for water and light. These are relatively fast processes, occurring on ecological time scales, which determine the emerging community. At longer times this community is further shaped by slow processes of interspecific competition among species of similar traits and by trait diffusion (mutations). A final remark about phenotypic changes: although in general 𝜒 can be interpreted as representing different phenotypes, the choice of very small values for 𝐷_𝜒 cannot represent relatively fast phenotypic changes and restricts the context to mutations at the taxonomic level of species.

We added an explanation in the 3^rd paragraph of the section “Modeling spatial assembly of dryland plant communities” of the need to consider mutations and the role they play in our study.

• Combining the observations from the previous two points, there is a concern that for a given level of precipitation, there is a single trait class with optimal biomass/lowest soil water level that is dominant, with the neighboring trait classes being sustained by the diffusion of biomass from the optimal class to neighboring inferior classes. This would seem a bit problematic, as it would mean that most classes are not a true fit for the environment, and only persist due to the continuous inflow of biomass. Taking a clue from the previous papers of the authors, it seems this may not be the case, though. Specifically, in the paper by Nathan et al. (2016) it seems that all trait classes are started at low initial biomass density, and the resulting steady state (in the absence of biomass flows between classes) seems to show similar biomass profiles as shown in Figures 4,5 and 7 of the current paper. While the current model formulation seems slightly different, similar results may apply here. Indeed, keeping all trait classes at non-zero (but low) density, and when the (abiotic and biotic) environment permits, let each class increase in biomass seems like the most straightforward approach to model community assembly dynamics. Given the above discussion about these trait classes competing for a single resource (soil water), and one trait class being able to drive this resource availability to the lowest level, it would then be useful to readers to explain why multiple trait classes can coexist here, and how (for spatial uniform solutions) the equilibrium soil water level with multiple trait classes present compares to the equilibrium soil water level when only the optimal trait class is present. Furthermore, if results as presented in Nathan et al. (2016) indeed hold in the current case, perhaps it means that the biomass profile responses as shown in e.g. Figure 5 would also occur if there was no biomass flow between trait classes included, but that the time needed to adjust the profile would take much longer as compared to when the drift term/second trait derivative is included. In summary, further clarification of what the biomass flows between classes represent, and the role it plays in driving the presented results would be useful for readers.

As explained in the reply to previous comments the asymptotic community is tuned by a balance between two slow counter-acting processes, interspecific competition among similar traits and mutations over evolutionarily long time scales. However, the community structure is largely determined by much faster processes of environmental filtering and interspecific competition among widely distinct traits, as all traits are initially present. Indeed, comparing the biomass distributions in new Figure 3, with and without trait diffusion indicates that the community composition, as measured by 𝜒_𝑚𝑎𝑥, is the same. Trait diffusion, however, does affect functional diversity, along with environmental factors. In that sense the emerging community is a true fit for the environment.

We thank the reviewer for these thoughtful comments, which helped us realize that our presentation of these issues was too concise and unclear. We believe that the new extended section on modeling spatial assembly of dryland plant communities, and the new figure 3a clarify these issues.

• In addition, it would be useful for readers to understand to what extent the shifts in average trait values and functional diversity can be decoupled from the biomass and soil water responses to changes in precipitation that would occur in a model with only a single biomass variable. For example, early studies on self-organization in semi-arid ecosystems already showed that the shift toward a patterned state involved the formation of patches with higher biomass, and higher soil water availability, as compared to the preceding spatially uniform state, and that the biomass in these patches remains relatively stable under decreasing rainfall, while their geometry changes (e.g. Rietkerk et al. 2002). It has also been observed that for a given environmental condition, biomass in vegetation patches tends to increase with pattern wavelength (e.g. Bastiaansen and Doelman 2018; Bastiaansen et al. 2018). Given the model formulation, one wonders whether higher biomass in the single variable model is not automatically corresponding to higher abundance of faster growing species and a higher functional diversity (as the diffusion of biomass can cover a broader range when starting from higher mass in the optimal trait class). There are some indications in the current work that the linkage is more complicated, for example, the biomass peak in Figure 7c is lower, but also broader as compared to the distribution of Figure 7b, but it is currently not entirely clear how this result can be explained (for example, it might be the case that in the spatially patterned states, the biomass profiles also vary in space).

We are not sure we understand what the reviewer means by “decoupled”, but much insight indeed can be gained from a study of a model for a single functional group (trait) and observing the behaviors described by the reviewer. In fact, these behaviors, which some of us are familiar with from numerical studies, motivated parts of the current study. Higher biomass in vegetation patches (compared to uniform vegetation) in the single trait model does not automatically imply a shift to faster growing species; in principle the stress-tolerant species that already reside in the system when uniform vegetation destabilizes to a periodic pattern can simply grow denser. To answer this and additional questions we need to take into account interspecific interactions by studying the full community model. As to Figure 7b,c, the behavior appears to be opposite to that described by the reviewer: the biomass pick in Figure 7c is higher and narrower than that in Figure 7b, not lower and broader. This is because of the much larger domain of the patterned state as compared with that of the uniform state, which increases the abundance of low-𝜒 species, i.e. species investing in growth.

The increase of biomass in vegetation patches with pattern wavelength for given environmental conditions, as observed by Bastiaansen et al. 2018, is actually another mechanism for increasing functional diversity. This is because the water stress at the patch center is higher than that in the outer patch areas and thus forms favorable conditions for stress tolerant species while the outer areas form favorable conditions for fast growing species.

We added a new paragraph in the Discussion and Conclusion section (last paragraph in the subsection Insight III) where we discuss the effect of coexisting periodic patterns of different wavelengths on functional diversity and ecosystem management. We also added citations to the references the reviewer mentioned.

• The possibility of hybrid states, where part of the landscape is in a spatially uniform state, while the other part of the landscape is in a patterned state, is quite interesting. To better understand how such states could be leveraged in management strategies, it would be useful if a bit more information could be provided on how these hybrid states emerge, and whether one can anticipate whether a perturbation will grow until a fully patterned state, or whether the expansion will halt at some point, yielding the hybrid state. It seems that being able to distinguish these case would be necessary in the design of planning and management strategies.

The hybrid states appear in the bistability range of the uniform and patterned vegetation states, and typically occupy most of this range. Their appearance is related to the behavior of ‘front pinning’ in bistability ranges of uniform and patterned states in general. Front pinning refers to fronts that separate a uniform domain and a periodic-pattern domain, which remain stationary in a range of a control parameter (precipitation in our case). This is unlike fronts that separate two uniform states, which always propagate in one direction or another and can be stationary only at a single parameter value – the Maxwell point. Thus, an indication that a given landscape may have the whole multitude of hybrid states is the presence of a front (ecotones) that separates uniform and patterned vegetation. If that front appears stationary over long period of times (on average), this is a strong indication.

We added a new paragraph in the subsection Insight III of the Discussion and conclusion section to clarify this point.

Also, in Figure 3a, the region of parameter space in which hybrid states occur is not very large; it is not entirely clear whether the full range of hybrid states is left out here for visual considerations, or whether these states only occur within this narrow range in the vicinity of the Turing instability point.

As pointed out in the reply to the previous comment the hybrid states are limited to the bistability range of uniform and patterned vegetation, which is not wide. However, this should not necessarily restrictma nagement of ecosystem services by nonuniform biomass removal, as such management will have similar effects on community structure also outside the bistability range where front propagate slowly.

The new paragraph we added also addresses this point.

Comments for the authors:

• Line 17: the term "re-patterning" may read as a non-patterned state becoming patterned again, whereas here it seems to refer to a spatial rearrangement of an existing patterned state.

We change “spatial re-patterning” to “spatial self-organization”

• Line 39: resources (i.e. plural)?

Fixed.

• Lines 80-99: This is an introduction to the model description, rather than a result, as the header suggests.

• Lines 100-164: This is the model description, which seems to be part of the material and methods rather than a result.

We added a new section “Methods” and moved the more technical part of the model description to that section.

• Line 179: when χ increase from 0.95 to 1.00 however, it seems that the Turing threshold start to increase, how can this reversal be explained?

This is a good question, which we are not sure we have a definite answer to, because the analytical expression for the Turing threshold is too complicated to be helpful in this regard. A possible explanation, however, is the following. Spatial patterning in this model is driven by a high infiltration contrast between vegetation patches and bare soil, as this increases the flow speed of overland water. Thus, high infiltration contrast (obtained for patches of high biomass density) favors higher thresholds 𝑃_𝑇. This explains why low 𝜒 species, which invest in growth, have higher 𝑃_𝑇. However, the shift to higher 𝜒 species at low P not only implies less investment in growth and thus lower biomass density, but also weaker water uptake which leaves more soil-water content for growth. In any event the non-monotonic behavior of 𝑃_𝑇 at the high 𝜒 range does not affect the results we present as these are all related to behaviors occurring for 𝜒 values well below that range.

• Lines 302-310: this explanation is clear, but it is an example that can also be explained by the biomass dynamics of a single variable model.

This explanation relies indeed on the dynamics of a single functional-group model and is used to argue that similar results about community dynamics should be obtained in 2d where morphological transitions can take place, a case we do not study in this work.

• Line 329: this is case where it would be useful for readers to understand how one can anticipate the formation of either hybrid or fully patterned states, and how this relates to the particular perturbation(s) imposed.

This point is discussed in a new paragraph we added in this section (paragraph next to last).

• Figures: Why are the biomass values in Figure 4,5 and 7 about an order of magnitude higher than in Figure 3?

We thank the reviewer for raising this point. In fact the biomass values in Figure 4 (now 3), 5 and 7 should be an order of magnitude lower than in in Figure 3 (now 4). The reason is that there are about 30 functional groups in the community that share the water content and therefore the biomass of each functional group should be much smaller compared to the case of a single functional group that benefits from the whole amount of water (Figure 3 in original manuscript). This confusion is a consequence of a different definition of the biomass variables used in the community model vs. that used in the single functional-group model. A biomass variable 𝐵_𝑖 in Equation 1 represents biomass density in trait space, which we now denote as 𝑏_𝑖, rather than the biomass of the i^th functional group, which is 𝐵_𝑖 = 𝑏_iΔ𝜒 = 𝑏_𝑖/𝑁, where 𝑁=128 in our study.

We rewrote the community model – Equation 1 (appearing now in the new Methods section) in terms of the new biomass variables 𝐵_𝑖 = (𝑏_iΔ𝜒) and updated the vertical biomass axes in Figures 4 (now 3), 5 and 7. We also added an explanation that the single functional-group model with 𝜒 = 1 used to produce Figure 3 (now 4) follows from the community model by setting all biomass variables identically to zero apart of 𝐵_𝑁 which corresponds to 𝜒_𝑁 = 𝑁Δ𝜒 = 1 (and dropping the trait-diffusion term which does not have a meaning for a single functional-group model).

Reviewer #2:

[…]

The conclusions of this paper are mostly well supported by data, but some aspects of model presentation, parameter choices, and data interpretation need to be clarified and extended.

1) Model presentation. It would be better to explain the model in ecological terms first, clarifying parameter biological meaning and justifying their choice. In doing so, creating a specific 'Methods' section, which now is lacking, would be of help too. Authors should clarify whether and how the model follows the conservation of mass principle involving precipitation and evapotranspiration. Are root growth and seed dispersal included for this purpose? Why are they not referred to any further in the analysis and discussion? Why a specific term for plant transpiration is not included, or is to somehow phenomenologically incorporated into the growth-tolerance tradeoff? In doing so, authors should also pay attention to water balance as above (H) and below (W) ground water are not independent from each other.

We added a Methods section, which in eLife is placed at the end of the manuscript. The section includes the model equations and more detailed explanations in ecological terms of various parts of the model. We also added Table 1 with a list of all model parameters, their descriptions, units and numerical values used in the simulations. Presenting the model at the end of the manuscript suits more technical information about the model, but not essential information that is needed for understanding the results. We therefore kept the subsection “A model for spatial assembly of dryland plant communities” in the Results section, where we present that information.

There is no conservation of mass in the model (and all other models of this kind) simply because the system that we consider is open. In particular, it does not include the atmosphere, which constitute part of the system’s environment. Including the atmosphere as additional state variables in the model, capturing the feedback of evapotranspiration on the atmosphere, would make the model too complicated for the kind of analysis we perform. So, although the model contains parts that represent mass conservation such as the terms describing below- and above-ground water transport, water mass is not conserved. The biomass variables represent aboveground biomass of living plants or plant parts and are not conserved either as biomass production involve biochemical reactions that convert inorganic substances coming from the system’s environment (atmosphere and the soil) into organic ones, while plant mortality involves organic matter that leaves the system.

Roots in the model platform we consider are modeled indirectly through their relation to aboveground biomass. That relation constitutes one of the scale-dependent feedbacks that produce a Turing instability to vegetation patterns, the so-called root-augmentation feedback (see Meron 2019, Physics Today), but in this particular study we eliminate this feedback for simplicity. The scale-dependent feedback that we do consider is the so-called infiltration feedback, associated with biomass-dependent infiltration rate that produces overland water flow towards vegetation patches, as explained in the subsection “A model for spatial assembly of dryland plant communities”. It will be interesting indeed to extend the study in the future to include also the root-augmentation feedback.

We assume short-range seed dispersal and take it into account through biomass “diffusion” terms (obtained as approximations of dispersal kernels assuming narrow kernels). These terms play important roles in the scale-dependent feedback that induces the Turing instability, as is explained in earlier papers which we cite. Plant transpiration is modeled through the water uptake term in the equation for the soilwater 𝑊. Indeed above-ground water 𝐻 and below-ground water 𝑊 are not independent; the infiltration term IH in the equations for both state variables account for this dependence in a unidirectional manner (loss of 𝐻 and gain of 𝑊). As we do not include the atmosphere in the model the other direction, namely, evapotranspiration that increases air humidity and affects rainfall, is not accounted for. The neglect of this effect can be justified for sparse dryland vegetation.

These good points have already been discussed in many earlier papers as well as in the book Nonlinear Physics of Ecosystems (Meron 2015), and we cannot address them all in this paper. We did however add several clarifications in the section Modeling spatial assembly of dryland plant communities and in the new Methods section, including the consideration of the atmosphere as the system’s environment quantified by the precipitation parameter 𝑃.

Another unclear point is that growth rates for the same plant functional groups are assumed to be constant among different species within the same group and are confounded by biomass production. Why is that the case? Furthermore, how many different species are characterizing each functional group? How are interspecific interactions accounted for (more specifically, see comment below)?

In the trait-based approach we focus on just two functional traits, related to growth rate and tolerance to water stress, ignoring differences in other traits that distinguish species. That is, a given functional group consists of species that share the same values of the two selected functional traits (to a given precision determined by 𝑁), taking all other traits represented in the model to be equal. In this approach we do not care about how many species belong to each functional group, only their total biomass. We wish to add that simplifying assumptions of this kind are necessary if we want the model to be mathematically tractable and capable of providing deep insights by mathematical analysis.

We expanded the discussion of the trait-based approach in the section Modeling spatial assembly of dryland plant communities and added relevant references (second paragraph).

Finally, stress tolerance is purely phenomenological. There is no actual mechanism/parameter describing it. Rather, it "simply" appears as low/high mortality, which in turn is said to be due to high/low tolerance. This leads to a sort of circularity between mortality and tolerance. Yet, mortality can occur due to other biophysical factors (e.g. disturbance, fire, herbivory, pathogens). A drawback of this assumption is that a mechanism of drought tolerance is often to invest in belowground organs, including roots. However, according to the proposed model, it turns out that fast growing species with low investment in tolerance also have high investment in roots; viceversa, tolerant species have low investment in roots. This is a bit counterintuitive and not well biologically supported.

First, we agree with the reviewer that our approach is purely phenomenological, as we model tolerance to water stress by a single parameter that lumps together the effects of various physiological mechanisms. That parameter can be distinguished from other factors affecting mortality by regarding the constant 𝑀_𝑚𝑎𝑥 in Equation 3 as representing several contributions. Since we do not study the effects of these other factors we can absorb them in 𝑀_𝑚𝑎𝑥 for mathematical simplicity.

Tolerance to water stress is not necessarily associated with roots. Plants can better tolerate water stress by reducing transpiration through stomatal closure, regulating leaf water potential, or develop hydraulically independent multiple stems that lead to a redundancy of independent conduits and higher resistance to drought (see Schenk et al. 2008 – https://doi.org/10.1073/pnas.0804294105).

We added a discussion in the Methods section (5^th paragraph, “Tolerance to water stress …”) of the simple form by which we model tolerance to water stress through the mortality parameter.

2) Parameter choice.

N = 128 is an extremely high number for plant functional groups. It is even quite unrealistic to have 128 species per square meter, so this value is not very reasonable. Please run the model and report results with more realistic N (e.g. from 4-64) as well as with different sets of N values keeping all other parameters constant.

We wish to clarify two points: (1) N=128 does not imply 128 functional groups per square meter; the emerging community has much lower functional richness (FR) as the average FR is around 0.25, meaning only 128 × 0.25 = 32 functional groups. (2) The model results, as reflected by the key metrics 𝜒_𝑚𝑎𝑥, 𝐹𝑅, and 𝐹𝐸, are independent of the particular value of N (for N values sufficiently large), as Author response image 1A and B show. The biomass 𝐵_𝑖 of each functional group, however, does change (Author response image 1A) because by changing N we change the range of traits Δ𝜒 = 1/𝑁 that belong to a given functional group. But if we look at the biomass density in trait space 𝑏_𝑖, related to 𝐵_𝑖 through the relation 𝐵_𝑖 = 𝑏_𝑖Δ𝜒, then also the biomass density is independent of 𝑁 as Author response image 1B shows. So, even if in practice there are less functional groups and thus species as considered in the model studies, the results are not affected by that. On the other hand, choosing higher 𝑁 values provides smoother curves and nicer presentation of our results.

We added a discussion of this issue in the Methods section after Equation 2.

Author response image 1

Download asset Open asset

Gamma (rate of water uptake by plants' roots): why is it in that unit of m²/kg * y? Why are you now considering the area (and not the volume) per biomass unit?

The vegetation pattern formation model we study, like most other models of this kind, does not explicitly capture the soil depth dimension. Accordingly, W is interpreted as the soil-water content in the soil volume below a unit ground area within the reach of the plant roots. In practice W has units kg/m², like B, and since Γ𝑊𝐵 should have the same units as 𝜕𝑊/𝜕𝑡 see Equation 1b, Γ must have the units of (𝐵𝑡)⁻¹.

A is not defined in the text.

We now define it in Table 1 (see Methods section).

M min: why 0.5 mortality? Having M max set to 0.9, please consider a lower mortality value set to 0.1, and please report evidence (hopefully) demonstrating the robustness of results to such change.

The results are robust to the particular values of 𝑀_𝑚𝑖𝑛 and 𝑀_𝑚𝑎𝑥, except that there are combinations of these two parameters for which the biomass distributions are pushed towards the edge of the 𝜒 domain, which make the presentation of the results less clear. Author response image 2 shows results of recalculations of the distribution 𝐵 = 𝐵(𝜒) for 𝑀_𝑚𝑖𝑛 = 0.1, as requested (using 𝑀_𝑚𝑎𝑥 = 0.15) for 3 different precipitation values. As the reviewer can see there’s no qualitative change in the results: lower precipitation push a uniform community to stress tolerant species (higher 𝜒), while the formation of patterns at yet lower precipitation push the community back to fast growing species (low 𝜒).

Author response image 2

Download asset Open asset

K_min and K_max are in two different units, and should both be kg/m².

Thanks, we fixed this typo in Table 1.

Values of precipitation (P, mean annual precipitation) are not reported.

The precipitation parameter is variable, as is now stated in Table 1, and therefore was not include it in the list of parameters’ values used. Whenever a particular precipitation value has been used our intention was to state it in the caption of the corresponding figure. This was done in Figures 5,6,7, but indeed not in Figure 4 (Figure 3 in revised manuscript). The insets on the right side of Figure 3 (Figure 4 in revised manuscript) where also calculated for particular precipitation values, but that information is not essential as the intention is to show typical forms of the various solution branches, which do not qualitatively change along the branches (i.e. at different P values).

We added the precipitation value (P=180mm/y) at which all the biomass distributions shown in new Figure 3 (Figure 4 in original manuscript) were calculated.

3) Results presentation and interpretation.

Parameter range of precipitation in figure 3 is odd. Why in one case precipitation ranges from 0 to 160 while in another it is only 60-120? Furthermore, in paragraph 198-213 and associated results in Figure 5. the Choice of precipitation values is somehow discordant from the previous model. Please provide motivation for this choice, clarify and uniformize it.

In Figure 3b (Figure 4b in revised manuscript) we restricted the precipitation range to 60-120 as the curves, which are limited to 0 < 𝜒 < 1 (by the definition of 𝜒), do not extend to 𝑃 < 60 and to 𝑃 > 120. Extending the range to 0 < 𝑃 < 160 would make the figure less compact and nice as it will contain blank parts with no information.

We are not sure we understand what the reviewer means by “is somehow discordant from the previous model”. The motivation of the choices we made for the precipitation values P=150, 100 and 80 was to show the shift of a spatially uniform community to a higher 𝜒 value as the precipitation is decreased to a lower value (from 150 to 100), and the shift back to a lower 𝜒 value at yet lower precipitation (80) past the Turing instability.

Throughout the text, authors claim to address plant-plant interactions, particularly intra and interspecific competition. However, it is not clear how competition was modelled neither whether it was included in the model. In its current state, it is just an assumption pulled out when discussing results – a classic 'passepartout' used by ecologists. Furthermore, why only competition is invoked in interpreting results when facilitation is known to be much more relevant in pattern formation and biodiversity maintenance in arid systems?

Finally, authors seem to create confusion around community composition, which is defined as the (taxonomic) identity of all different species inhabiting a community. Notably, it is remarkably different from the x_max parameter used in the model, which as a matter of facts is just the value of the most productive (notably, not necessarily the most abundant) functional group.

We thank the reviewer for this comment. Since all the emerging communities in the model studies are pretty localized around the value of 𝜒_𝑚𝑎𝑥, that value does contain information about the identity of other functional groups in the community when complemented by FR (functional richness) and FE (functional evenness). More significantly to our study, shifts in 𝜒_𝑚𝑎𝑥 represent the shifts in community composition we focus on in this study, i.e. shifts towards fast growing species or towards stress-tolerant species.

We modified the description of the community-level properties that can be derived from the biomass distribution in trait space (see modified text towards the end of the section “Modeling spatial assembly …” and also the caption of Figure 3b), explaining that both functional diversity and community composition can be described by several metrics, and clarifying the significance of 𝜒_𝑚𝑎𝑥 in describing community composition shifts.

Reviewer #3:

In this paper, the authors use a mathematical model of plant and water dynamics in drylands to show that drylands adaptive capacity to respond to changes, via spatial self-organization in space has also beneficial effects in preserving its biodiversity and ecosystem functions.

The model is an extension of a large body of previous, well-established works on plant self-organisation in drylands. The model is well described and motivated (with one main exception, see below), the analyses are robust and the results are very convincingly supporting the conclusions. I however have an issue with one of the assumptions in the model equations. The authors included a term for "mutations" in traits that (1) is not introduced or motivated (2) its effects/importance are not highlighted by specific analyses (3) the possible implications or limitations connected to it are not discussed. To my knowledge, this term is also not based on earlier work. All these elements need to be included, as at the moment is for example unclear what the authors intended to represent by including the mutation term (evolutionary time scales? Or adaptation?). Also, it would be especially good to include an analysis of how influential this term is for the final results.

Assuming the authors can address this one concern, the results are surely important as they connect for the first time plant spatial self-organization to its biodiversity preservation, in the face of future expected climatic changes and probable land degradation. These findings, although theoretical, have the potential to be useful also for guiding adaptive and dynamic land management, as they underline the importance of taking into account spatial vegetation distribution in drylands management.

Besides the major point about the mutation term, I list here two other important points:

We addressed this issue in our reply to the first reviewer and summarize our reply as follows. Without the mutation term interspecific competition among all functional groups (assumed to be present initially) results in the emergence of a community centered around a most abundant functional group 𝜒_𝑚𝑎𝑥 on a relatively fast ecological time scale. This functional group eventually outcompetes all other functional groups, but that competition becomes extremely slow as the competing functional groups become similar. On these long-time scales mutations is a process that cannot be ignored. Including mutations at evolutionarily slow rates (𝐷_𝜒 = 10⁻⁶) results in the emergence of an asymptotic community with a characteristic distribution of traits (functional diversity) around the same 𝜒_𝑚𝑎𝑥. These results are also displayed graphically in new Figure 3.

– The authors state that they represent highly tolerant plants by representing the plants with a small mortality. However, in their model, plant mortality does not depend on soil water levels. How can the authors reconcile these two aspects? Also, one could argue that mortality is related to the average life span, not specifically to tolerance to highly stressful condition. The authors should better justify this point and discuss the implication of this assumption.

The relation to the soil-water level is as follows. The fitness of the ith functional group is Λ_𝑖𝑊 − 𝑀_𝑖. This means that functional groups with lower mortality can survive (i.e. have positive fitness) lower soil-water contents or higher water stress. This simple modeling form applies to situations where stress tolerance is not a phenotypic trait switched on by water stress. An example would be plants with hydraulically independent multiple stems that lead to a redundancy of independent conduits and higher resistance to drought (see Schenk et al. 2008 - https://doi.org/10.1073/pnas.0804294105). We could model stress tolerance as a phenotypic trait, as was done in Tzuk et al. 2019, but that would not change the main results and message. There are obviously additional factors contributing to plant mortality besides water stress. However, since we do not study these factors in this work we could lump their contributions in the mortality parameter 𝑀_𝑚𝑎𝑥.

We added a discussion of these issues in the (new) Methods section (paragraph beginning with “Tolerance to water stress is modeled …”).

– In the model, there is shading feedback too, not only infiltration feedback. However the authors state there's only infiltration feedback in l. 84, could they please explain?

The shading feedback (reduced evaporation by shading) is not a scale-dependent feedback that leads to spatial patterning, as it does not involve water transport. It rather leads to bistability of states.

Author details

1. Bidesh K Bera

   Department of Solar Energy and Environmental Physics, BIDR, Ben-Gurion University of the Negev, Sede Boqer Campus, Israel

   Contribution

   Formal analysis, Investigation, Methodology, Project administration, Writing - review and editing

   Competing interests

   No competing interests declared

2. Omer Tzuk

   Physics Department, Ben-Gurion University of the Negev, Beer Sheva, Israel

   Present address

   Department of Industrial Engineering, Faculty of Engineering, Tel-Aviv University, Tel Aviv-Yafo, Israel

   Contribution

   Conceptualization, Formal analysis, Investigation, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-6541-3311

3. Jamie JR Bennett

   Department of Solar Energy and Environmental Physics, BIDR, Ben-Gurion University of the Negev, Sede Boqer Campus, Israel

   Contribution

   Formal analysis, Investigation, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-9748-5010

4. Ehud Meron

   1. Department of Solar Energy and Environmental Physics, BIDR, Ben-Gurion University of the Negev, Sede Boqer Campus, Israel
   2. Physics Department, Ben-Gurion University of the Negev, Beer Sheva, Israel

   Contribution

   Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Writing - review and editing

   For correspondence

   ehud@bgu.ac.il

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-3602-7411

• 891

  Page views

• 192

  Downloads

• 7

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.