Impact statement
Dryland ecosystems support unique biodiversity of flora and fauna and many human populations for crucial ecosystem services. Dryland vegetation can exhibit multiple stable states, such as savannas and woodlands. They can also sometimes abruptly switch between these states, even with little environmental changes, making them difficult systems to predict and manage. Mathematical models have played a crucial role in elucidating insights into dynamics of such systems, but to date, they all have ignored a key biological reality: Individuals both within and among species vary from each other in their phenotype (or traits).
In this study, we fill this major gap. Our model analyses, which incorporate diversity in plant demographic and functional traits, reveal that these variations in savanna species favour woodlands over grasslands, that is, grasslands exist over a smaller range of rainfall levels. Depending on which trait is varied, trait variations can make the transition from one state to another more, or less, abrupt. Our study also shows that changes in tree cover and community composition could be decoupled due to strong impacts of historical contingency effects, such that differing community compositions can coexist for the same tree cover. Our modelling study suggests that the management of savannas could benefit from a trait-based dynamics approach.
Introduction
Savannas are intriguing systems where grasses and trees coexist without competitively excluding each other (Scholes and Archer, Reference Scholes and Archer1997). Savannas differ from woodlands in that they have a continuous layer of grass with a discontinuous layer of trees, giving rise to a predominantly grassy and fire-prone state. In contrast, in woodlands, trees form a continuous layer alongside a discontinuous grassy layer, giving rise to a high canopy cover (as referred to in Walker and Noy-Meir, Reference Walker, Noy-Meir, Huntley and Walker1982; Wilcox et al., Reference Wilcox, Birt, Fuhlendorf and Archer2018). These states differ not only in their physiognomy but the fire regimes they experience. As a consequence of these complex interactions, savannas and woodlands can coexist as alternative stable states under the same environmental conditions (Walker and Noy-Meir, Reference Walker, Noy-Meir, Huntley and Walker1982; Dublin et al., Reference Dublin, Sinclair and McGlade1990; Van Langevelde et al., Reference Van Langevelde, Van De Vijver, Kumar, Van De Koppel, De Ridder, Van Andel, Skidmore, Hearne, Stroosnijder and Bond2003; Beckage et al., Reference Beckage, Platt and Gross2009; Warman and Moles, Reference Warman and Moles2009; Hirota et al., Reference Hirota, Holmgren, Van Nes and Scheffer2011; Staver et al., Reference Staver, Archibald and Levin2011). Systems with two alternative stable states – also known as bistable systems – show complex dynamics in response to changes in external conditions or driver (Holling, Reference Holling1973; May, Reference May1977). Their response can be nonlinear and discontinuous rather than linear and continuous. Such a response to external conditions, which may also vary stochastically with time, can result in an abrupt shift from one stable state to another (Scheffer et al., Reference Scheffer, Carpenter, Foley, Folke and Walker2001; Guttal and Jayaprakash, Reference Guttal and Jayaprakash2007; Chen et al., Reference Chen, Jayaprakash, Yu and Guttal2018). In the context of savanna-woodland ecosystems, studies show that rainfall is a key driver of ecosystem states (Sankaran et al., Reference Sankaran, Hanan, Scholes, Ratnam, Augustine, Cade, Gignoux, Higgins, Le Roux and Ludwig2005; Hirota et al., Reference Hirota, Holmgren, Van Nes and Scheffer2011; Staver et al., Reference Staver, Archibald and Levin2011; Eby et al., Reference Eby, Agrawal, Majumder, Dobson and Guttal2017). While grassland and savannas are found to occupy regions with low rainfall, woodlands occupy regions with high rainfall. Interestingly, both biome states can be found at intermediate levels of rainfall. When rainfall exceeds a certain threshold value, the system can abruptly shift from a grassy state to a woody state. Such shifts may dramatically alter biodiversity, ecosystem function and the services they provide to large human populations that depend on savanna ecosystems (Solbrig et al., Reference Solbrig, Medina and Silva1996; Jackson et al., Reference Jackson, Banner, Jobbágy, Pockman and Wall2002; Bond and Parr, Reference Bond and Parr2010). In this broad context, our goal is to investigate the role of trait diversity, an important feature of real systems, on the bistability of the savanna-woodland system.
Mathematical models (Higgins et al., Reference Higgins, Bond and Trollope2000; D’Odorico et al., Reference D’Odorico, Laio and Ridolfi2006; Meyer et al., Reference Meyer, Wiegand and Ward2009; Baudena et al., Reference Baudena, D’Andrea and Provenzale2010; Yatat et al., Reference Yatat, Dumont, Tewa, Couteron and Bowong2014; Majumder et al., Reference Majumder, Tamma, Ramaswamy and Guttal2019; Djeumen et al., Reference Djeumen, Dumont, Doizy and Couteron2021) have proved to be important tools for understanding ecosystems with alternative stable states, where manipulative experiments, such as long-term fire manipulation experiments (Biggs et al., Reference Biggs, Biggs, Dunne, Govender and Potgieter2003; Higgins et al., Reference Higgins, Bond, February, Bronn, Euston-Brown, Enslin, Govender, Rademan, O’Regan and Potgieter2007) or landscape-scale fire-exclusion experiments (Andersen et al., Reference Andersen, Cook, Corbett, Douglas, Eager, Russell-Smith, Setterfield, Williams and Woinarski2005) are difficult to test various hypotheses at large scales. While the literature is replete with models describing such complex systems and their dynamics, very little is known about how trait diversity affects the dynamics of bistable systems (Dakos et al., Reference Dakos, Matthews, Hendry, Levine, Loeuille, Norberg, Nosil, Scheffer and De Meester2019). It is well known that variation between species can be high; however, for certain traits, variation between individuals of a species can be comparable to variation among individuals of different species (Des Roches et al., Reference Des Roches, Post, Turley, Bailey, Hendry, Kinnison, Schweitzer and Palkovacs2018). For example, in a meta-analysis of various plant communities, intraspecific variation in plant height was found to be as high as 67% of the total within-community trait variance (Siefert et al., Reference Siefert, Violle, Chalmandrier, Albert, Taudiere, Fajardo, Aarssen, Baraloto, Carlucci and Cianciaruso2015). These differences between individuals can be more or less pronounced depending on the scale or level at which these comparisons are made (Messier et al., Reference Messier, McGill and Lechowicz2010). Intraspecific or interspecific trait variation can have a stronger effect on ecological responses if they influence these responses indirectly or directly, respectively (Des Roches et al., Reference Des Roches, Post, Turley, Bailey, Hendry, Kinnison, Schweitzer and Palkovacs2018). However, these trait differences are often obscured in studies that treat all individuals as identical to each other or use only the mean value of the trait (Moran et al., Reference Moran, Hartig and Bell2016). While this assumption is convenient for simplifying study design and modelling, such assumptions could potentially limit our understanding of how communities and ecosystems persist and how they respond to changes in their environment (Bolnick et al., Reference Bolnick, Amarasekare, Araújo, Bürger, Levine, Novak, Rudolf, Schreiber, Urban and Vasseur2011).
In this context, it is pertinent to ask how trait variation can affect ecosystem structure, function, and dynamics. A few experimental (Agashe, Reference Agashe2009; Becks et al., Reference Becks, Ellner, Jones and Hairston2010), meta-analytic (González-Suárez et al., Reference González-Suárez, Bacher and Jeschke2015) and theoretical studies (Hart et al., Reference Hart, Schreiber and Levine2016; Crawford et al., Reference Crawford, Jeltsch, May, Grimm and Schlägel2019; Baruah, Reference Baruah2022) have used phenomenological models to show how trait variation affects species diversity, community structure and (or) ecosystem processes. Recent articles have also speculated on the role of individual variation in vector-borne disease dynamics (Cator et al., Reference Cator, Johnson, Mordecai, Moustaid, Smallwood, LaDeau, Johansson, Hudson, Boots, Thomas, Power and Pawar2020), tipping points in ecosystems (Dakos et al., Reference Dakos, Matthews, Hendry, Levine, Loeuille, Norberg, Nosil, Scheffer and De Meester2019), and tritrophic food chains (DeAngelis, Reference DeAngelis2013). Cator et al. (Reference Cator, Johnson, Mordecai, Moustaid, Smallwood, LaDeau, Johansson, Hudson, Boots, Thomas, Power and Pawar2020) stressed the importance of using trait-based mechanistic models over phenomenological ones as the latter do not account for the underlying mechanism and can be unreliable over longer timescales. How significant are these individual differences to the various processes and patterns we see in nature is a question of paramount importance (Funk et al., Reference Funk, Larson, Ames, Butterfield, Cavender-Bares, Firn, Laughlin, Sutton-Grier, Williams and Wright2017).
In the context of savanna-woodland bistable systems, it is unclear whether trait diversity would increase or decrease the chances of abrupt shifts (or “tipping events”) and how easily these ecosystems/biome states can be restored to the previous state after a tipping event. For instance, in the case of mutualistic networks, Baruah (Reference Baruah2022) showed that trait variation may increase the occurrence of tipping events. Similarly, a model of shallow lake considers trait evolution and reveals that the resilience of the system increases through adaptive evolution (Chaparro Pedraza et al., Reference Chaparro Pedraza, Matthews, de Meester and Dakos2021). Another recent study (Limberger et al., Reference Limberger, Daugaard, Gupta, Krug, Lemmen, Van Moorsel, Suleiman, Zuppinger-Dingley and Petchey2023) showed that depending on the functional trait that was varied, the resilience of the system with alternative stable states increased, decreased or exhibited no change. However, the effect of trait variation on savanna-woodland dynamics remains understudied.
Here, we consider a well-studied model of savanna-woodland bistable ecosystem (Staver and Levin, Reference Staver and Levin2012) and incorporate trait variation to analyse its effects on the ecosystem dynamics. We introduce heritable variation in two demographic traits (the death rates of savanna saplings and the death rates of savanna trees) and in a functional trait that nonlinearly determines the transition of savanna saplings to adults, referred to as sapling resistance to fire. We look at how these trait variations affect the state of the system and the qualitative nature of these shifts at the ecosystem level, for traits that vary linearly (death rate of savanna sapling and trees) and non-linearly (sapling resistance to fire) with growth rates. We also look at the distribution of traits in the system and which individuals comprise the steady-state population.
Model and methods
Savanna-woodland model with no trait variations
We begin by introducing the model developed by Staver and Levin (Reference Staver and Levin2012) which describes the dynamics of the ecosystem consisting of grasses, savanna tree saplings and savanna tree adults. This model does not incorporate trait variations and thus serves as a baseline model for our study. Here, the proportion of grass cover, savanna saplings and adult trees are denoted by
$ G $
,
$ S $
and
$ T $
, respectively, which together add up to unity. The model incorporates how these three components of the ecosystem interact differently with fires: Fire spread and frequency of the fire depend on the grass cover; at high grass cover, the system experiences frequent fires that spread across the landscape. During fires, grasses and savanna tree saplings lose their aboveground biomass (Higgins et al., Reference Higgins, Bond and Trollope2000). However, when the saplings reach a certain height such that they can escape the flame zone, they become resistant to subsequent fires (Hoffmann et al., Reference Hoffmann, Geiger, Gotsch, Rossatto, Silva, Lau, Haridasan and Franco2012). This fire-resistant stage is termed as the savanna tree adult. In this way, fire limits the transition of savanna tree saplings to adults.
The above nonlinear/threshold behaviour of savanna saplings to fire frequency is captured by the rate constant
$ \omega \left(G;\theta \right) $
, which is a nonlinear function of the grass cover as shown in Figure 1(A). At high values of grass cover
$ G $
, due to the higher frequency and spread of fires, the transition of savanna tree saplings to adult savanna trees occurs at a very low rate, while it occurs at a higher rate for low grass cover. The value of the parameter
$ \theta $
can be used to determine the approximate value of grass cover at which nonlinearity kicks-in and can be used as a proxy to determine the resistance of the sapling to fire. We note that the fire is modelled implicitly, the model is deterministic and nonspatial (also referred to as mean-field models). We relax some of these assumptions in our second version of the trait diversity model.
Figure 1. (A) Savanna sapling to tree transition rate (
$ \omega $
) for different values of sapling resistance fire (
$ \theta $
).
$ \omega $
varies as a function of Grass cover (
$ G $
) and
$ \theta $
. At low values of grass cover and thus low fire frequency, saplings have a high rate of transition to adult trees and vice versa. Saplings with high
$ \theta $
value (
$ \theta $
= 0.8, denoted by blue solid line) transition to trees even at high grass cover while those with low
$ \theta $
(
$ \theta $
= 0.3, denoted by green dashed line) value remain at the sapling stage even at low grass cover. (B) Stability or bifurcation diagram of the savanna-woodland model with no trait variations (no variation model). The figure depicts how the steady-state grass cover changes as a function of sapling birth rate (
$ \beta $
), different stable states (grassland in green, savanna in orange and woodland in blue), threshold points (denoted by star – P1, P2, P4), the extent of bistable region (grey region) and the quantum of shift at threshold points (Inset). Parameters: sapling resistance to fire
$ \theta $
= 0.5, tree death rate
$ \nu $
= 0.1 and sapling death rate
$ \mu $
= 0.5.
Other key parameters of the model are as follows: trees produce saplings at a rate
$ \beta $
, with saplings replacing grasses as they are assumed to be competitively superior to grasses. Saplings and adult trees die at a constant intrinsic death rate
$ \mu $
and
$ \nu $
, respectively. The grasses spread instantaneously and cover the space previously occupied by saplings or trees. Putting these processes in a mathematical form yields the model described by three coupled ordinary differential equations (Staver et al., Reference Staver, Archibald and Levin2011) (see Table 1).
Table 1. Symbols and descriptions of the model and model variables and parameters
In savanna-woodland systems, rainfall is often the key environmental driver of the vegetation states. Here, we assume that the savanna sapling birth rate (
$ \beta $
) increases with rainfall (Wilson and Witkowski, Reference Wilson and Witkowski1998). Thus, a smaller value of
$ \beta $
corresponds to a higher aridity stress on the system. We characterise the bifurcation diagram which illustrates how the steady states (grass cover) change as a function of the driver value (
$ \beta $
; see Figure 1(B)), identifying their stability. In this system, we define three ecosystem states: the term grassland state refers to a system entirely occupied by grasses (~100% grass cover), whereas a system with more than 50% tree cover is referred to as the woodland state. The state with a grass cover less than 100% but greater than 50% is termed a savanna state.
Figure 1(B) shows that at low sapling birth rates (high aridity), the system exists only in the grassland state; we refer to this as the monostable grassland state. In contrast, at high sapling birth rates, the steady state of the system is woodland, also referred to as monostable woodland state. At intermediate
$ \beta $
values (denoted by the grey region,
$ 0.34<\beta <1.9 $
), we find the existence of alternative stable states, also called a bistable region. In this region, we find either the coexistence of grassland and woodland states (for
$ 0.34<\beta <1.1 $
) or savanna and woodland states (for
$ 1.1<\beta <1.9 $
). In the bistable region, depending on the initial value of the grass cover, the system will reach one of the two alternative states. The range of parameters over which the system is bistable is referred to as the extent of bistability (grey horizontal bar in Figure 1(B)).
In this system, we observe several types of transitions from one stable state to another. Figure 1(B) shows that at the bifurcation or the threshold point marked P1 in the upper branch, the system exhibits a continuous transition from a grassland to a savanna state as a function of the increasing sapling birth rate
$ \beta $
; we refer to this as a gradual regime shift. In addition, at a higher
$ \beta $
value, marked P2 in Figure 1(B), the system exhibits a discontinuous transition from a savanna state to a woodland state; we refer to such transitions as abrupt regime shifts. We see another abrupt regime shift at the threshold point P4, from a woodland to a grassland state as a function of decreasing
$ \beta $
values. We can quantify the transitions for each of these regime shifts by the difference in grass cover between the two stable states at the threshold point; we term this as the quantum of shift. For a gradual regime shift, this quantity will be zero whereas it will be a nonzero number for abrupt transitions. Larger the quantum of shift, more drastic the regime shift.
Trait variation models
In this manuscript, our goal is to understand how the dynamics at the ecosystem level and population level are affected by trait diversity. At the ecosystem level, we investigate the total vegetation cover (grass, tree or sapling) which is obtained by summing over proportions of all traits of that vegetation type. Specifically, we first obtain the bifurcation diagram which captures the relationship between the state variable of the ecosystem (e.g. grass cover or tree cover) and the birth rate of the sapling (
$ \beta $
); as remarked before, the parameter
$ \beta $
can be thought of as a proxy of annual rainfall, a key driver of vegetation in semi-arid ecosystems. This will reveal how properties of the bifurcation diagram such as the extent of bistability (i.e. grey region in Figure 1(B)), qualitative nature of transitions, as well as the quantum of shifts (inset of Figure 1(B)) at the threshold point vary when we introduce trait variations in the model. Next, we also look at the consequences of trait variation at the population level, that is, we not only focus on total grass or tree cover, but also how the population of each of the traits evolves over time, and which tree and sapling traits constitute the steady-state population.
To investigate these questions, we consider two modelling approaches. First, we use the deterministic modelling framework of Staver and Levin (Reference Staver and Levin2012) and include discrete variations in three traits of the model. Second, we use a stochastic agent-based modelling framework where the same three traits are considered, but are now assumed to vary continuously and evolve via demographic stochastic processes. Both modelling approaches assume clonal inheritance of traits.
Deterministic discrete trait variation model
We describe an extension of the model of Staver and Levin (Reference Staver and Levin2012) to include trait variations in two demographic traits (sapling death rate (
$ \mu $
), tree death rate (
$ \nu $
)) and one functional trait (sapling resistance to fire –
$ \theta $
; defined below). First, we define
$ \theta $
, which appears in the sapling-to-tree transition rate function
$ \omega \left(G;\theta \right) $
, as the trait that determines the resistance of the sapling to fire. Although the original model treats
$ \omega \left(G;\theta \right) $
as a sapling trait that is entirely dependent on the ecosystem-level properties of fire governed by the grass cover, we treat
$ \theta $
to capture how resistant individual saplings are to fires. Saplings with a high value of
$ \theta $
can transition to adult trees even when fires are frequent, whereas saplings with a low value of
$ \theta $
can transition to trees only if the frequency of the fire is low (Figure 1(A)). This trait, for example, could be correlated with the height of the saplings such that taller saplings survive topkill during fires and thus transition to adult trees at a higher rate (Higgins et al., Reference Higgins, Bond and Trollope2000, Reference Higgins, Bond, Combrink, Craine, February, Govender, Lannas, Moncreiff and Trollope2012).
We now introduce variations in demographic traits – tree death rate (
$ \mu $
) and sapling death rate (
$ \nu $
) – which result from a combination of a suite of traits rather than one specific plant trait. For example, tree death in savannas can occur due to megaherbivores such as elephants (Vanak et al., Reference Vanak, Shannon, Thaker, Page, Grant and Slotow2012), which plants may resist by investing in spiny or toxic stems (Sheil and Salim, Reference Sheil and Salim2004). Another major factor for tree or sapling deaths is drought (Fensham et al., Reference Fensham, Fairfax and Ward2009; Sankaran, Reference Sankaran2019). Trees with low drought resistance mediated by traits such as low cavitation resistance, low hydraulic safety margins and low wood density experience high mortality during droughts (Sankaran, Reference Sankaran2019) and thus have high tree death rate (
$ \nu $
) values. In addition, saplings that invest in bark thickness and stem diameter are more likely to survive in the sapling stage (Hoffmann and Solbrig, Reference Hoffmann and Solbrig2003). Some of the traits that help saplings survive from background death processes may also help survive fires, and hence be correlated with
$ \theta $
, a point we ignore for the study since we only vary one trait at a time. We emphasise a key difference between the functional and demographic traits we have chosen to vary: whereas the trait
$ \theta $
representing sapling resistance to fire affects the growth rate of trees via nonlinear frequency dependence on grass cover as shown in Figure 1(A), the sapling death rate (
$ \mu $
) and tree death rate (
$ \nu $
) are independent of frequency of any of the ecosystem states (grass cover, sapling or tree covers) and affect the corresponding growth rate term linearly.
For simplicity, we assume that each trait follows a discrete uniform distribution (see Sensitivity Analysis for other distributions), by considering a finite number
$ k $
of trait types, which could represent different tree species at one extreme or a continuous trait approximated as discrete types for mathematical convenience. We assume that the traits are clonally inherited, with no mutations. If we consider trait variations in only one of the traits keeping the other traits fixed, the corresponding mathematical equations will be in the form of
$ 2k+1 $
ordinary differential equations (Table 1).
In the Results section, we present results for two levels of trait variations – a ‘low variation’ case and a ‘high variation’ case, keeping the mean trait value the same. The low trait variation is implemented as a smaller range of the trait values (0.3–0.7), while the high trait variation corresponds to a larger range (0–1). We use a combination of analytical methods and numerical simulations, details of which are presented in Supplementary Material Section 1; see Table 1 for the description and values of parameter values.
Sensitivity analysis
To test the robustness of our results, we perform sensitivity analyses. First, we consider different number of discrete trait types (
$ k $
) to be 10, 100 or 1000. In addition, we also consider different trait distributions – unimodal beta and bimodal beta distributions, in addition to the case of uniform distribution already discussed. Finally, we also consider a model variation in which Forest trees are included, as in the full model of (Staver and Levin, Reference Staver and Levin2012) thus making the total vegetation types to four (see Supplementary Material Section 2 for model details).
Stochastic continuous trait variation model
To further test the robustness of our conclusions on our choice of discrete trait types, we employ an agent-based stochastic dynamical model. This is an alternative and complementary modelling strategy that helps us relax some major assumptions of the model in the previous section: by modelling traits as continuous (rather than discrete) and by including stochastic birth – death processes (rather than deterministic). Specifically, here too, we consider three types of individuals: grass (
$ \mathcal{G} $
), savanna saplings (
$ \mathcal{S} $
) and savanna trees (
$ \mathcal{T} $
). A tree
$ i $
with trait value
$ {x}_i $
(can be either
$ {\mu}_i $
,
$ {\nu}_i $
or
$ {\theta}_i $
) can give rise to a sapling offspring having the same trait value, that is, clonal inheritance, replacing a grass ‘individual’. A sapling individual can grow into a tree individual, maintaining the trait value. Both tree and sapling individuals can die at their respective death rates (
$ {\mu}_i $
and
$ {\nu}_i $
) and will be replaced by grass. These processes, in the notation of chemical kinetics, can be represented as:
$$ \mathcal{T}\left({x}_i\right)+\mathcal{G}\overset{\hskip1em {\beta}_i\hskip1em }{\to}\mathcal{T}\left({x}_i\right)+\mathcal{S}\left({x}_i\right) $$
$$ \mathcal{S}\left({x}_i\right)\overset{\omega \left(G,{\theta}_i\right)}{\to}\mathcal{T}\left({x}_i\right) $$
$$ \mathcal{S}\left({x}_i\right)\overset{\hskip1em {\mu}_i\hskip1em }{\to}\mathcal{G} $$
$$ \mathcal{T}\left({x}_i\right)\overset{\hskip1em {\nu}_i\hskip1em }{\to}\mathcal{G} $$
Each reaction is a random process whose rate depends on the trait values of the individuals. We highlight an important contrast with a typical reaction kinetic framework of population models, where the rates of reactions are constant. In our model, each individual
$ i $
has their own rates of transitions (
$ {\mu}_i $
,
$ {\nu}_i $
or
$ {\theta}_i $
– which we interpret as traits) and the distribution of these rates change over time. We use the Gillespie algorithm (Gillespie, Reference Gillespie1976) to simulate this process for a population size of
$ N=1000 $
and for 1,000 time units. We note that the system reaches a steady state typically within 200 time units.
Results
Ecosystem-level dynamics of the deterministic model
We first present the results of the deterministic trait variation model at the ‘ecosystem-level’. We do this by investigating the grass cover, tree cover (and sapling cover) at steady-state for various possible variations in traits. Specifically, we consider how properties of the bifurcation diagram change with trait variations (no, low and high variation, as already described). In Figure 2, we display how the bifurcation diagram changes with trait variations. We reveal four major effects: First, we find that with an increase in variation in any of the three traits, the grassland state (i.e. the state with
$ G=1 $
and
$ T=0 $
) exists for a shorter range of the sapling birth rate (
$ \beta $
); in particular, the grassland state will be found for lower values of
$ \beta $
(row A, Figure 2). In contrast, trait variations have the opposite impact on the woodland state: woodlands can occur for a larger range of sapling birth rates and even at lower values of
$ \beta $
(i.e. at higher aridity levels). In other words, trait variations seem to reduce bistability in favour of woodland state.
Figure 2. Ecosystem-level properties such as bistable region and nature of regime shifts can depend on trait variations. For each bifurcation diagram, we assume that only one trait (mentioned in the title of that diagram) exhibits variations keeping other traits constant. (Row A) shows that the extent of bistable region reduces with increasing trait variation (see the panels below each subfigure). Extent of grassland reduces, while that of woodland increases, with increasing trait variations. (Row B) shows that the quantum of abrupt shift reduces with increasing trait variation. (Row C) shows how nature of transition changes with trait variations: In (C1), whereas the no variation case (denoted by closed circles) exhibits a transition from woodland to grassland, the high variation case (denoted by open triangles) shows a transition from woodland to savanna. In (C2), the qualitative nature of transitions remains same for both no and high variations. In (C3), while the no variation case shows a transition from woodland to grassland, the high variation case shows a transition from savanna to grassland. Parameters: (A1–C1) sapling resistance to fire
$ \theta $
= 0.5 and tree death rate
$ \nu $
= 0.1. (A2–C2)
$ \theta $
= 0.5 and sapling death rate
$ \mu $
= 0.05. (A3–C3)
$ \nu $
= 0.1 and
$ \mu $
= 0.2.
Second, we find that increasing trait variations leads to a reduction in the bistable region in the system, that is, the region of parameter space where the two alternative states are stable (see the bottom panels of row A, Figure 2). For example, in the no variation case for the sapling death rate
$ \mu $
(Figure 2(A1)), the bistable region ranges between 0.34 and 1.90 units of
$ \beta $
. This reduces to a range of 0.28–0.37 units of sapling birth rate (
$ \beta $
) for the high-variation case. This is true for all the three cases of trait variation we studied (Figure 2(A1), (A2), (A3)).
Third, we observe qualitative changes in the nature of regime-shifts with increasing trait variation, but these patterns differ between traits in which variation exists (Figure 2, row C). For instance, when there is no variation in the sapling death rate (
$ \mu $
; Figure 2(C1)), we find an abrupt shift from a woodland to a grassland state as a function of decreasing
$ \beta $
value, whereas, with high trait variation, we find an abrupt shift from a woodland to a savanna state. No qualitative change in regime shifts is observed when we introduce variations in tree death rate (
$ \nu $
); the system shifts from a woodland state to a savanna state as a function of decreasing
$ \beta $
for both no and high variation cases (Figure 2(C2)). In contrast, for the high variation in the resistance of the saplings to fire (
$ \theta $
), we observe an abrupt shift from savanna to grasslands (Figure 2(C3)).
Fourth and finally, we also observe changes in the quantum of shifts in grass cover at the threshold point, with increasing trait variations (Figure 2, row B). The quantum of shift decreases across all the cases of increasing variation in
$ \mu $
,
$ \nu $
and
$ \theta $
(Figure 2(B1)–(B3)). Thus, with higher trait variation in the system, regime shifts become less drastic compared to the no-variation case.
Population-level dynamics of the deterministic model
We now consider population-level dynamics, which we define as explicit consideration and analysis of trait distributions for various scenarios (Figure 3). Specifically, we compare the initial and final distribution of the traits in each of the three regimes (monostable grassland, bistable and monostable woodland as previously described), to understand the population composition at steady state. We begin by noting that the bifurcation diagram, that is, steady states denoted by grass cover and tree cover, of the no trait variation model was independent of initial grass cover
$ {G}_0 $
or the distribution of traits in the monostable woodland regime. However, in the trait variation model, despite the ecosystem-level property being independent of initial conditions, we find that population-level states can depend on initial values: When we consider the variation in the resistance of the saplings to fire (
$ \theta $
; Figure 3, column C), we find that the proportion of each persisting tree trait depends on the initial grass cover (
$ {G}_0 $
), highlighting that the dynamics of the population level can depend on the initial conditions. Thus, woodlands with equivalent ecosystem-level tree cover can have different compositions of tree (and sapling) types in the population, depending on their initial conditions. However, we find that when the sapling or tree death rates vary (Figure 3, columns A and B, respectively) only trees and saplings with the lowest trait value survive, while others get eliminated. This is true for all initial values of grass cover (rows representing different initial grass cover; also see Supplementary Figure S4 for more values of initial grass cover, and S5 and S6 for different initial distributions). This exclusion of individuals with higher death rates results in a final population with no variation, as expected. Unlike the monostable woodland regime, results from the bistable regime do not show a strong pattern (see Supplementary Figure S7) while the monostable grassland regime has no tree population.
Figure 3. Population-level dynamics can be sensitive to initial conditions, for
$ \beta $
= 0.45 which corresponds to monostable woodland regime. Specifically, we show how the uniform initial trait distribution (shown in the top most row) evolves for different initial values of grass cover (across rows). We find that for variations in sapling death rate (
$ \mu $
; A1–A3) and tree death rate (
$ \nu $
; B1–B3), only the sapling and trees with the least death rates (thus, higher survival) survive. However, for variation in sapling resistance to fire (
$ \theta $
; C1–C3), we find that the steady-state distribution depends on the initial grass cover, with mean value of
$ \theta $
increasing as initial grass cover increases.
Our analysis of temporal evolution further reveals the importance of trait variations (Figure 4). Whereas the system with no variation case yields a grassland state (Figure 4(A)), the system with trait variation in sapling resistance to fire (
$ \theta $
) reaches a woodland state (Figure 4(B)). A closer examination of the dynamics at the trait level (Figure 4(C)–(E)) further reveals that not all the initial tree/sapling types survive. At
$ t=0 $
, trees with high
$ \theta $
values start increasing in proportion, as their saplings immediately transition to adults, resulting in a sharp decline in the corresponding sapling cover. Consequently, sapling cover of these trees increases, reducing the grass cover. With the decline in grass cover, other sapling types with lower
$ \theta $
eventually transition to trees, while types with low
$ \theta $
values decline to zero asymptotically due to death processes at sapling death rate
$ \mu $
. Trait variation also results in the system reaching a stable state earlier than when no variation is present.
Figure 4. A comparison of time series for the no variation and the trait variation model with variation in sapling resistance to fire (
$ \theta $
), at a fixed value of sapling birth rate. Lines in red and blue hues represent tree cover (
$ T $
) and sapling cover (
$ S $
), respectively, while the dark grey line represents Grass cover (
$ G $
). (A) shows that the system exists in a grassland state for the model with no variation. (B) shows that the system exists in a woodland state for the model with trait variation. (C) shows the individual proportions of different tree and sapling types for the trait variation model in (B). The red and blue bubbles above the lines determine the
$ \theta $
value of the specific tree and sapling type, respectively. (D) and (E) show a zoomed-in view of the grey region in (C), showing only tree and sapling cover, respectively. Parameter values at t = 0 are: total tree cover = total sapling cover = 0.25, while individual tree and sapling cover = 0.025.
$ \nu $
= 0.1,
$ \mu $
= 0.2,
$ \beta $
= 0.45,
$ \theta $
= [0.05, 0.15, 0.25, …, 0.95].
Robustness of results to parameter variations and model variations
Our analysis with alternative trait distributions (unimodal beta, Supplementary Figure S1; bimodal beta, Supplementary Figure S2) and increasing the number of trait types (Supplementary Figure S3) shows that our conclusions on how trait diversity affects the ecosystem level dynamics are robust. In addition, our simulations of the stochastic continuous trait variation model (Supplementary Figure S8) and the extended model that includes forest trees (Supplementary Figure S11) too confirm our findings. In all these cases, as trait variation increases, the extent of grassland decreases and that of the woodland increases. Further, the quantum of the jump in the abrupt transition also reduces as trait variation increases. In the stochastic model, we note that the bistable region naturally shrinks due to the stochasticity, as expected from previous studies on stochastic bistable systems (Guttal and Jayaprakash, Reference Guttal and Jayaprakash2007). Furthermore, we find that the population-level dynamics exhibit qualitatively similar historical contingency effect irrespective of the trait distributions chosen (Supplementary Figures S4–S6), including the case of the stochastic continuous trait model (Supplementary Figures S9 and S10) and the model that includes forest trees (Supplementary Figure S12).
Discussion
Our model study shows that trait diversity can influence the bistable savanna-woodland system in several interesting ways. We find a general pattern across three traits that we studied: the range of bistability reduces when the trait variation increases. More specifically, a higher trait variation leads to a larger parameter range where woodlands persist. Consequently, the parameter range of stability for the grassland reduces. Furthermore, we find that the qualitative nature of transitions and the quantum of transitions can change, making them less drastic, due to trait variations. At the population level, we find a novel result that the steady-state distribution of traits depends on the initial conditions of the ecosystem; this is true even when the steady-state ecosystem state (e.g. tree cover or grass cover) is identical. Broadly, this suggests that historical contingency has a stronger impact at the population level rather than at the ecosystem level when trait diversity is considered.
Different trait variations may also affect the dynamics at the population and ecosystem levels differently. When we considered trait variations in death rates, we found that the ecosystem in steady state consists of trees (or saplings) with the lowest death rates, as expected. However, in the case of sapling resistance to fire, we find that trait types that are above a threshold value survive, while others are eliminated from the population (Figure 3(C1)–(C3)). This is because saplings with
$ \theta $
values below that threshold fail to transition to trees due to fire. Therefore, saplings with a value of
$ \theta $
lower than the threshold experience strong selection, while those with
$ \theta $
above the threshold experience weak or no selection pressure, leading to high variation in the trait in steady state. A similar finding of certain traits being absent in the steady-state population has been reported by Dantas et al. (Reference Dantas, Pausas, Batalha, Paula Loiola and Cianciaruso2013) in the Cerrados of eastern Brazil where fire prevents certain species from colonising the savannas. We argue that the response of the trait sapling resistance to fire is different from sapling and tree death rates possibly due to the nonlinear relationship between the trait and survival from fire (Moran et al., Reference Moran, Hartig and Bell2016); this result is also broadly consistent with the effect predicted by Jensen’s inequality (Ruel and Ayres, Reference Ruel and Ayres1999).
Interestingly, we find that communities with similar environmental conditions can be very dissimilar in their community composition, not only as a result of complexity resulting in alternative stable states but also due to historical contingency even when there are no multiple stable states. In the woodland regime, the steady-state tree cover for different cases of initial conditions was identical, yet the population compositions differed for trait
$ \theta $
(Figure 3(C1)–(C3)). This difference results from the influence of initial grass cover on steady-state distribution. When grass cover (which governs the frequency of fire) is reduced, even saplings with lower
$ \theta $
find conditions favourable for them to transition to trees (Figure 4(C)–(E)); note that this interplay is possible only because of initial trait diversity. Thus, even in the absence of bistability, the community shows a strong impact of historical legacy of grasses. This impact of entirely different taxa (i.e., grasses) on tree community trait assemblage is mediated by disturbance frequency experienced by the system. In this way, historical disturbance regimes would have a signature on the current community composition.
Potential practical implications
We discuss some potential practical insights from our study. Across continents, tropical grasslands and savannas are at risk of woody encroachment (O’Connor and Crow, Reference O’Connor and Crow1999; Ratnam et al., Reference Ratnam, Bond, Fensham, Hoffmann, Archibald, Lehmann, Anderson, Higgins and Sankaran2011; Stevens et al., Reference Stevens, Lehmann, Murphy and Durigan2017; Buisson et al., Reference Buisson, Le Stradic, Silveira, Durigan, Overbeck, Fidelis, Fernandes, Bond, Hermann and Mahy2019). One method employed for grassland restoration includes prescribed fires and mechanical removal of trees (or tree thinning) (Brockway et al., Reference Brockway, Gatewood and Paris2002; Cuevas and Zalba, Reference Cuevas and Zalba2010). Tree thinning is expected to achieve the effect of reduced resistance to fire, but requires careful consideration of the dynamics of the system as the restoration efforts can also lead to an undesirable state (Smit, Reference Smit2004). In fact, our model analyses show that the removal of highly fire-resistant saplings (i.e. high
$ \theta $
values (0.65–0.95) in the steady-state population represented by Figure 3(C1)) is insufficient to revert a woodland to a grassland state. To restore a grassland, it is necessary to also remove saplings with high resistance to fires (either manually or through controlled prescribed burning). In other words, restoration efforts can be made more efficient by focusing on the removal of both tree and sapling types. Our work adds credence to the idea that understanding the dynamics of grassland systems, along the role of trait diversity, is an important research gap for achieving better grassland restoration (Török et al., Reference Török, Brudvig, Kollmann, Price and Tóthmérész2021).
While most studies look at the coarse-grained properties of the system, our trait diversity model can be used to gain insights on community composition. We find that changes in tree cover and plant diversity are decoupled in the monostable woodland regime when the resistance of the saplings to fire trait (
$ \theta $
) exhibits variability. Such decoupling was also observed by Pardo et al. (Reference Pardo, Camarero, Gutiérrez and García2013) for one of the Pyrenean treeline ecotones as even with no significant change in tree cover, vegetation composition changed significantly. In the literature on savanna forest ecosystems, remotely sensed tree cover is often used as a state variable that can be used to assess ecosystem resilience as well. Changes in community composition, which arise from historical contingency and other local factors, do not affect the tree cover, but can in turn increase the vulnerability of the overall ecosystem due to changes in trait distributions. Therefore, our study cautions that resilience metrics that use tree cover alone (or their derivative properties such as variance or correlations (Génin et al., Reference Génin, Majumder, Sankaran, Danet, Guttal, Schneider and Kéfi2018; Dakos et al., Reference Dakos, Matthews, Hendry, Levine, Loeuille, Norberg, Nosil, Scheffer and De Meester2019)), without accounting for trait diversity and community composition, could have serious limitations as indicators of ecosystem stability (Baruah et al., Reference Baruah, Clements and Ozgul2020). At a very broad level, owing to differences in the predictions of no trait variations model and the more realistic trait variation model (Figure 1(B) and Figure 2), an important take away for ecosystem modellers is to consider trait diversity, not just mean trait values (Wakeling et al., Reference Wakeling, Staver and Bond2011), in ecosystem modelling.
Future directions and concluding remarks
We have adopted a simple model in which we investigated variation in one trait at a time. Future studies can incorporate variations in multiple traits simultaneously. There could be interesting trade-offs, for example, a higher resistance to fire (
$ \theta $
) through bark/stem properties can result in a slower growth rate of saplings to trees (Corrêa Scalon et al., Reference Corrêa Scalon, Domingos, da Cruz, Júnior, Marimon and Oliveras2020). Or, there could be positive correlations between different traits (e.g. a higher resistance to fire may also lead to better survival to other background abiotic stressors). Intriguingly, even in the deterministic version of the trait variation model, we predict coexistence of a diverse trait distribution of sapling resistance at steady state. These variations can be attributed to nonlinear response of vegetation to climate (rainfall) and disturbances (fires through grass cover). Naturally, including additional mechanistic factors that act as the source of trait variation (e.g. demographic stochasticity and dispersal) as well as external stochastic factors (e.g. seasonal and interannual rainfall variability) offer some interesting avenues for future work.
On the modelling front, any system with multiple competing traits is essentially an eco-evolutionary dynamical system, for which many alternative and mathematical rigorous formulations are available. These include, but are not limited to, parsing the interactions via game theoretic pay-off matrices and constructing suitable replicator-mutator equations or an adaptive dynamics approach (Geritz et al., Reference Geritz, Metz, Kisdi and Meszéna1997; Lion, Reference Lion2018; Wickman et al., Reference Wickman, Koffel and Klausmeier2023; Bhat and Guttal, Reference Bhat and Guttal2025). In addition, mutations and spatial structure to the model of trait diversity are known to have impact on ecological dynamics especially in the context of dryland ecosystems (Hirota et al., Reference Hirota, Holmgren, Van Nes and Scheffer2011; Staver et al., Reference Staver, Archibald and Levin2011; Sankaran et al., Reference Sankaran, Majumder, Viswanathan and Guttal2019; Goel et al., Reference Goel, Guttal, Levin and Staver2020; Bennett et al., Reference Bennett, Bera, Ferré, Yizhaq, Getzin and Meron2023). Another important area in the context of systems with bistable dynamics is the ability to obtain early warning signals of potential abrupt transitions (Scheffer et al., Reference Scheffer, Carpenter, Foley, Folke and Walker2001; Génin et al., Reference Génin, Majumder, Sankaran, Danet, Guttal, Schneider and Kéfi2018). How do trait variations affect, do they enhance or reduce the intensity of signals, could be a question of important practical consideration.
In summary, we have shown some potential effects that trait variations can have on dryland ecosystems with alternative stable states. Trait diversity may promote woodlands, make regime shifts less catastrophic and lead to the establishment of a community assembly that may crucially depend on the historical composition of the community. Finally, beyond dryland ecosystems, a diverse range of ecosystems are hypothesised to exhibit alternative stable states. Our research suggests that incorporating trait diversity in mathematical models could open up new avenues to gain insights on the basic ecological understanding of ecosystems and potentially offers new perspectives for management and restoration.
Supplementary material
To view supplementary material for this article, please visit http://doi.org/10.1017/dry.2025.1.
Data availability statement
All codes necessary to reproduce the results of this manuscript are available via our Git-hub repository: https://github.com/tee-lab/trait-variation/.
Acknowledgements
We are grateful to three anonymous reviewers for many helpful constructive criticisms which substantially improved the manuscript.
Author contribution
T.R. and V.G. conceived the idea. T.R. and D.R. developed the models, conducted simulations and analysed the results. S.P. performed model analyses. T.R. and V.G. wrote the manuscript, with inputs from other coauthors.
Financial support
We acknowledge funding from the Department of Biotechnology (BT/PR29719/FCB/125/19/2018), Department of Science and Technology (NGP-GS08/IISc-B/KA/2023 and Science and Engineering Research Board (MTR/2022/000273) to V.G. and a PhD scholarship to T.R. via Ministry of Education and Prime Ministers Research Fellowship, all from the Government of India.
Competing interest
The authors declare no competing interests.
Open access
Comments
Dear Editor,
We wish to submit an original research article entitled “Higher individual variation in savanna tree species may reduce bistability in favour of woodlands” for consideration for publication in the journal Drylands.
Many ecological systems exhibit alternative stable states and numerous models have been developed to study them. However, few of these theories/models incorporate a biologically important feature: individual trait variations. Here, we chose a dryland system, savanna-woodland ecosystem, as our system of study as they are important ecosystems that support huge biodiversity and many human populations. Savannas are shown to coexist with woodlands as alternative stable states, with complex dynamics. In this study, we incorporate variations in traits of savanna trees and saplings in a model that explains savanna-woodland bistable dynamics.
To the best of our knowledge, this is one of the few models that incorporates trait variations in a bistable ecosystem and perhaps the first in the contex of a dryland ecosystem. Our model is in the form of coupled ordinary differential equations. A thorough analysis of the model via numerical simulations shows that trait variations in savanna trees reduce bistability in a savanna-woodland bistable system. Further, trait variations can change the qualitative nature of regime shifts between savanna and woodlands. We also report how demographic traits and a functional trait affect the dynamics at the ecosystem and population levels differently.
We believe that this manuscript is appropriate for publication in Drylands because our study analyses a theoretical model of important dryland ecosystems (grassland, savanna and woodland) with a crucial biological feature that has largely been ignored: individual variations. Beyond drylands, this study caters to different fields such as theoretical ecology, ecosystem dynamics, tipping points, and the broad field of dryland ecology. We make comments on how our results can provide insights on restoration and management strategies.
We have no conflicts of interest to disclose. We confirm that this work is original and has not been published elsewhere, nor is it currently under consideration for publication elsewhere.
Thank you for your consideration of this manuscript.
Sincerely,
Tanveen Kaur Randhava and Vishwesha Guttal
Centre for Ecological Sciences
Indian Institute of Science, Bengaluru, 560012, Karnataka, India.