Received: 1 April 2016

|

Revised: 25 November 2016

|

Accepted: 3 December 2016

DOI 10.1111/geb.12570

RESEARCH PAPERS

Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem Stephanie Eby1,2 | Amit Agrawal3 | Sabiha Majumder4 | Andrew P. Dobson1,5 | Vishwesha Guttal3 1 Department of Ecology and Evolutionary Biology, Princeton University, Princeton, NJ, USA 2

Department of Marine and Environmental Sciences, Northeastern University, Boston, MA, USA 3

Centre for Ecological Sciences, Indian Institute of Science, Bengaluru, Karnataka, India 4

Department of Physics, Indian Institute of Science, Bengaluru, Karnataka, India 5

Santa Fe Institute, Santa Fe, NM, USA

Correspondence Vishwesha Guttal, Centre for Ecological Sciences, Indian Institute of Science, Bengaluru, Karnataka 560012, India. Email: [email protected]

Abstract Aim: Theory suggests that as ecological systems approach regime shifts, they become increasingly slow in recovering from perturbations. This phenomenon, known as critical slowing down [CSD], leads to spatial and temporal signatures in ecological state variables, thus potentially offering early indicators of regime shifts. Indicators using temporal dynamics have been empirically validated in laboratory microcosms and other well-mixed systems, but tests of spatial indicators of regime shifts at large spatial scales in the field are rare due to the relative absence of high-resolution data and difficulties in experimental manipulations. Here, we test theoretical predictions of CSD-based spatial indicators using large-scale field data from the Serengeti–Mara grassland–woodland system. Location: Serengeti–Mara ecosystem, Tanzania and Kenya. Time period: Year 2000 Major taxa studied: Vegetation

Editor: Margaret Mayfield Funding information Department of Biotechnology, Govt of India, DST-FIST program, Govt of India, ISRO-IISc Space Technology Cell, McDonnel Foundation.

Method: We used a space-for-time substitution method to empirically test the validity of CSDbased spatial indicators, i.e., we computed indicators along a spatial [in lieu of temporal] gradient of ecological states. First we used a model of vegetation dynamics to determine if our space-fortime substitution method was appropriate. Then we tested for CSD-based spatial indicators using high-resolution spatial vegetation [30 m] and rainfall [2.5 km] data from the Serengeti–Mara ecosystem. Results: Our model predicts that CSD-based indicators increase along a spatial gradient of alternative vegetation states. Empirical analyses suggest that grasslands and woodlands occur as alternative stable states in the Serengeti–Mara ecosystem with rainfall as one of the potential drivers of transitions between these states. We found that four indices of CSD showed the theoretically expected increasing trends along spatial gradients of grasslands to woodlands: spatial variance, spatial skewness, spatial correlation at lag-1 and spatial spectra at low frequencies. Main conclusions: Our results suggest that CSD-based spatial indicators can offer early warning signals of critical transitions in large-scale ecosystems. KEYWORDS

alternative stable states, critical slowing down, critical transitions, early warning signals, grassland– woodland transitions, savanna ecosystems, Serengeti–Mara ecosystem, spatial ecology, spatial variance

Authors Stephanie Eby and Amit Agrawal contributed equally to the manuscript.

Global Ecol Biogeogr. 2017;1–12

wileyonlinelibrary.com/journal/geb

C 2017 John Wiley & Sons Ltd V

|

1

2

|

EBY

1 | INTRODUCTION

ET AL.

along spatial gradients [rather than temporal gradients] from one state to an alternative state match the theoretically expected patterns. The

A range of complex dynamical systems, from physiological to ecologi-

basic assumption of space-for-time substitution is that spatial patterns

cal, earth and climatic, can exhibit transitions from one stable state to

along spatial gradients of a driver are equivalent to the response of the

another (Scheffer, Carpenter, Foley, Folke, & Walker, 2001). Examples

ecosystem to a temporally changing driver. Although there are limita-

of ecological transitions include shifts from an oligotrophic to a eutro-

tions to this approach (Johnson & Miyanishi, 2008), it can offer useful

phic state in lakes, collapses of patchy vegetation to a bare state in dry-

insights when experimental manipulations in the field are not feasible

land ecosystems, woodland encroachment of savannas and shifts in

or long-term data are unavailable (Blois, Williams, Fitzpatrick, Jackson,

community composition in terrestrial and marine ecosystems (Bagchi

& Ferrier, 2013; Streeter & Dugmore, 2013). To test hypotheses

et al., 2012; Browning, Franklin, Archer, Gillan, & Guertin, 2014;

related to spatial indicators of ecosystem transitions, we chose a

!fi Carpenter & Brock, 2006; Carpenter, Ludwig, & Brock, 1999; Ke

savanna ecosystem, because these ecosystems contain grass- or

et al., 2007; Litzow, Urban, & Laurel, 2008; Ratajczak, Nippert, &

woody-dominated alternative stable states, with rainfall and fire acting

Ocheltree, 2014; Rietkerk, Dekker, de Ruiter, & van de Koppel, 2004;

as potential drivers of these patterns (Favier et al., 2012; Hirota,

Scheffer et al., 2001; Steele, 1998). Such changes, also referred to as

Holmgren, Van Nes, & Scheffer, 2011; Staver, Archibald, & Levin,

regime shifts or critical transitions, not only affect the ecology of these

2011; Figure 1). Savannas also cover an eighth of the global land sur-

systems but can also lead to loss of ecosystem services, thus affecting

face (Scholes & Archer, 1997), making them globally important ecosys-

the human societies that depend on them (Scheffer et al., 2001). This

tems that provide important ecosystem services [e.g., grazing for

has therefore created an interest in the development of early warning

livestock, carbon storage and tourism].

signals of ecological transitions using the mathematical framework of complex dynamical systems (Scheffer et al., 2009). Mathematical models show that as an ecological system approaches a regime shift, it takes longer for it to recover from perturbations (van Nes & Scheffer, 2007; Wissel, 1984). This phenomenon, known as critical slowing down [CSD], leads to increasing autocorrelation, variance, skewness and reddening of the power spectrum in the temporal dynamics of state variables before the system undergoes a regime shift (Carpenter & Brock, 2006; Guttal & Jayaprakash, 2008; Held & Kleinen, 2004; van Nes & Scheffer, 2007; Wissel, 1984). Based on the generality and robustness of these model predictions, these statistical indicators have been suggested as generic early warning signals of critical transitions (Dakos et al., 2012; Lenton, 2011; Scheffer et al., 2009). Empirical studies on diverse systems ranging from microbial populations grown in laboratories to whole lake ecosystems and largescale climatic shifts have lent support to this theory (Carpenter et al., 2011; Dai, Vorselen, Korolev, & Gore, 2012; Dakos et al., 2008; Lenton, 2011; Veraart et al., 2012). Theory predicts that analogous CSD-based indicators measured from spatial data, e.g., spatial variance, spatial skewness or spatial auto-

Here, we use a simple spatial model of regime shifts to demonstrate the validity of the space-for-time substitution approach as a way of testing spatial indicators of critical slowing down [Box 1]. We then employ high-resolution spatial data of vegetation and rainfall from the Serengeti–Mara savanna ecosystem in Tanzania and Kenya, East Africa (Reed, Anderson, Dempewolf, Metzger, & Serneels, 2009) [Figure 1]. We investigate one potential driver of transitions between vegetation types by looking at the relationship between vegetation type and mean annual rainfall at the large-scale landscape level and the smaller-scale transect level. We then test the hypothesis that the following spatial indicators of transitions will increase, as per theoretical predictions, along spatial transects from one vegetation type to the other [e.g., grassland to woodland vegetation]: spatial variance and spatial skewness (Guttal & Jayaprakash, 2009), spatial correlation (Dakos et al., 2010) and the spatial discrete Fourier transform at low frequencies (Carpenter & Brock, 2010).

2 | MATERIAL AND METHODS 2.1 | Data

correlation, are likely to offer robust measurements that anticipate

2.1.1 | Study site

regime shifts (Carpenter & Brock, 2010; Dakos, Nes, Van Donangelo,

The Serengeti–Mara ecosystem covers 25,000 km2. It is located in the

!fi et al., 2014). Fort, & Scheffer, 2010; Guttal & Jayaprakash, 2009; Ke

north-western part of Tanzania and the south-western part of Kenya

Empirical tests of spatial indicators require temporal sequences of spa-

and is composed of the Serengeti National Park and the Masaai Mara

tial snapshots of ecosystems before they undergo regime shifts. Recent

Reserve. Above-ground productivity in the region is controlled by rain-

work on spatial indicators has largely been limited to microbial systems

fall, which mainly falls between November and May (Norton-Griffths,

and simple aquatic microcosms on relatively small spatial scales (Dai,

Herlocker, & Pennycuick, 1975). A rainfall gradient exists across the

Korolev, & Gore, 2013; Drake & Griffen, 2010); thus, their validity for

ecosystem, with the north-west being the wettest part [1100 mm

field ecosystems remains unclear. Due to the lack of high-resolution

year21] and the south-east the driest part [400 mm year21] of the eco-

spatio-temporal data and difficulties in experimentally manipulating

system (Norton-Griffths et al., 1975; Sinclair, 1979). Vegetation in this

ecosystems in the field, testing for CSD-based spatial indicators in the

system can be influenced by rainfall and topological features (Reed

field has rarely been done [but see Cline et al. (2014)].

et al., 2009). Tree densities in the Serengeti woodlands declined from

To circumvent this problem, we take a “space-for-time” substitu-

the 1920s to the 1980s, but since then they have increased (Sinclair

tion approach where we test whether spatial indicators of transitions

et al., 2007). In the Masaai Mara, woodlands declined from 1950 to

EBY

|

ET AL.

3

F I G U R E 1 A map of the Serengeti–Mara vegetation classified as grasslands and woodlands overlaid with arrows indicating the location and direction of the nine sampling transects [Color figure can be viewed at wileyonlinelibrary.com]

1982 (Dublin, 1991); however, their status has not been systematically

maps created by M. Coughenour. The detailed methodology for proc-

monitored since then.

essing rainfall data is presented in Appendix S2.

2.1.2 | Vegetation data

2.2 | Rainfall–vegetation analyses

We used a ground-truthed vegetation map, at a resolution of 30 m 3 30 m, for the Serengeti–Mara ecosystem, which is composed of 14 different vegetation classifications (Reed et al., 2009). Using ArcMap (Environmental Systems Research Institute Inc. [ESRI], Redlands, CA, USA, version 10) we recoded the vegetation classifications at each pixel to two classes with binary values, woodlands [0] and grasslands [1] [Figure 1 and Table S1 in Appendix S1 of the Supporting Information].

First, we conducted a landscape-wide analysis of the relationship between mean annual rainfall and vegetation. For each rainfall grid of 2.5 km 3 2.5 km in the entire landscape, we computed the proportion of pixels classified as grass [referred to as grass cover] from the 30-m resolution vegetation data. Following the methods of Hirota et al. (2011) and Staver et al. (2011), we divided the entire spatial area of the Serengeti–Mara ecosystem into rainfall ranges; specifically, we chose 500–650, 650–800, 800–950 and 950–1100 mm. For each of the

2.1.3 | Rainfall data

rainfall ranges, we plotted the frequency of occurrence of different

We used a 30-year [1977–2006] mean annual rainfall GIS layer for the

grass covers. To test whether these frequency data followed a unimo-

Serengeti–Mara ecosystem with a 2.5-km grid based on yearly rainfall

dal or a multimodal distribution, we used Hartigan’s DIP test statistic

4

|

EBY

ET AL.

Box 1. A theoretical framework for space–time substitution and CSD-based spatial indicators We use a model that accounts for the simple ecological processes of establishment, death and facilitation among plants, which are key interactions shaping spatial patterns and vegetation regime shifts in savanna and semi-arid ecosys!fi et al., 2007; Lubeck, 2006; Rietkerk et al., 2004; Scheffer et al., tems (Ke 2001) [see Model and Update Rules below]. Previous studies on indicators of regime shifts assumed that the driver of the regime shift was uniform in space but gradually changed over time, resulting in the system undergoing an abrupt !fi et al., 2014). transition (Dakos et al., 2010; Guttal & Jayaprakash, 2009; Ke Here, to represent conditions relevant to the vegetation transects of Figure 1, we assume that the driver of vegetation in the model changes over distance. To test the validity of this space-for-time substitution, we measure a state variable [vegetation cover] and the following theoretically proposed CSDbased spatial indicators over distance: spatial variance, spatial skewness, spatial !fi et al., 2014). correlation at lag-1 and DFT at low-frequencies (Ke

Model and update rules We employed a cellular automaton model of vegetation in a twodimensional discrete space containing N 3 N grid cells. Each cell can be in one of two discrete states, empty [denoted by 0] or occupied by an individual [denoted by 1]. The state of each cell is updated depending on its state and those of its neighbours. If a randomly selected cell is occupied, we select another random cell in the nearest neighbourhood of the focal cell. If the neighbouring cell is empty, the focal plant reproduces in this empty cell with a probability p or the focal plant dies with a probability 1 – p. Here, we refer to the establishment probability [p] as the driver value and the death probability [1 – p] as the stressor value. On the other hand, if the neighbouring cell is also occupied, we select another neighbour randomly, and if that cell is empty the focal plant reproduces in this empty cell with a probability q. Note that the death probability of the focal plant is reduced in the presence of a neighbouring plant to a probability of [1 – q][1 – p], as opposed to [1 – p] in the absence of a neighbouring plant. Increased likelihood of establishment [or decreased likelihood of death] for an individual due to the presence of nearby plants, referred to as local positive interactions, may arise from processes such as increased water infiltration and reduced evaporation and run!fi et al., 2007; Rietkerk et al., 2004). Each time unit consists of off (Ke executing the above rules N2 times. If the positive feedback [q] is high, the system shows an abrupt regime shift from a high cover to zero cover [representing an alternative state] as a function of an increasing stressor [1 – p] [Figure S1].

Spatial gradient and results To mimic our vegetation transects, we model the spatial gradient in the driver values by decreasing the value of p over distance. To account for the fact that the vegetation data may not be in a steady state, we update the system for only 2000 times steps, which is insufficient for the model simulations to reach a steady state. Figure 2(a) shows that the vegetation cover goes from a high value of cover to zero cover as a function of the stressor. CSD-based spatial indicators, i.e., spatial variance, spatial skew-

A spatial model of vegetation dynamics used to predict the patterns of CSD-based spatial indicators within the space-for-time substitution framework. (a) The model exhibits a transition over distance from a high-cover to a lowcover state along a driver. Parts (b)–(e) show that the model predicts increasing trends of the following CSD-based spatial indicators prior to sharp decline in vegetation cover: spatial variance, spatial skewness, spatial autocorrelation at lag-1 and spatial low-frequency spectra [Color figure can be viewed at wileyonlinelibrary.com]

FIGURE 2

ness, spatial autocorrelation at lag-1 and spatial low frequency spectra, all

justification for our application of a space-for-time substitution method,

show clearly increasing trends [Figure 2b–e] over distance, even before

wherein we hypothesize that patterns of indicators along a spatially chang-

the onset of nonlinear changes in cover. These results provide a theoretical

ing driver are equivalent to those due to a temporally changing driver.

EBY

|

ET AL.

[using the “diptest” package in R 3.1.1; R Foundation for Statistical

5

allowed us to compare the trends of spatial indicators from these con-

Computing, Vienna, Austria]. In this test, the null model is that the data

trol transects with the transects described in the section Rainfall–vege-

are adequately described by a unimodal distribution. We identified the

tation Analyses above.

modes of the frequency distribution as the peak of the estimated kernel density function [the function density in R]. If the dip test rejected

2.3.2 | Coarse-graining of null model and transect data

the unimodal hypothesis we employed the method of finite mixtures of

We divided the N 3 N matrix (of vegetation or the null model) into

regression models (using the “flexmix” package in R 3.1.4) to test for

non-overlapping submatrices of size a 3 a, where a represents the size

the presence of multimodality in the data. We employed the Bayesian

of coarse-graining and obtained grass cover in each submatrix to con-

information criterion (BIC) to identify the goodness of fit between mix-

struct a ‘coarse-grained’ matrix (Figure S3). These matrices consist of

ture models with a single cluster and with multiple clusters. If the dip

continuous variables and not binary values, and they can be used to

test rejected unimodality, and if the mixture model suggested a multi-

identify if a transect’s spatial pattern is more ordered than a random

cluster model as a good fit, we concluded that the frequency distribu-

!fi et al., 2014). We chose a 5 5 (see the secpattern (Appendix S3; Ke

tion was best explained by a multimodal distribution.

tion Sensitivity Analyses below), which lowered resolution (from 30 m

We then tested the relationship between mean annual rainfall and

3 30 m to 150 m 3 150 m) as well as the dimension of the matrix

grass cover at the smaller transect scale, using transects with a width

(from 250 3 250 to 50 3 50). Unless otherwise stated we use these

of 7.5 km and lengths ranging from 40 to 90 km. The start- and end-

coarse-grained null model data and coarse-grained vegetation data for

point of the transects were randomly chosen from areas of grassland

all analyses in the rest of the paper. Table S2 summarizes the spatial

[visually identified as nearly full grass cover from Figure 1] and areas of

scales employed in our spatial indicator and rainfall–vegetation

woodland [visually identified as woody landscape from Figure 1],

analyses.

respectively. No transects are located in the south-eastern section of the ecosystem due to the presence of the Ngorongoro Crater High-

2.3.3 | Computing spatial indicators along transects

lands whose woody vegetation consists of montane forests located at

For any given sliding coarse-grained vegetation matrix [containing

a higher elevation than the rest of the woody vegetation in the ecosys-

2500 entries of grass cover at 150-m resolution], we used built-in func-

tem. Each transect was subdivided into sliding “vegetation matrices” of

tions in MATLAB to calculate mean grass cover [“mean”], spatial var-

approximate dimensions 7.5 km 3 7.5 km, with a sliding distance of

iance [“var”] and spatial skewness [“skewness”]. We computed spatial

2.5 km. With a base map resolution of 30 m 3 30 m, a vegetation

correlation at lag-1 and spatial DFT at low frequencies using standard

matrix of 7.5 km 3 7.5 km corresponded to 250 rows 3 250 columns.

definitions [Appendix S3.4].

Therefore, for example, a transect of 90 km would contain 34 sliding vegetation matrices.

For each transect, we computed spatial indicators starting from the mode of the frequency distributions of grass cover, as found in the

We computed the correlation coefficient [using the command corr

rainfall–vegetation analysis, to a decline in grass cover by 15%. If the

in MATLAB 7.1.1, The MathWorks, Inc.] between (1) distance along

number of sliding vegetation matrices for this range of grass cover was

the transect and rainfall [to test for a spatial gradient in rainfall], (2) dis-

less than five, we considered a larger range of grass covers, from the

tance along the transect and grass cover [to test for a vegetation spa-

mode to a decline of 35%. The mode of the distribution can be consid-

tial gradient], and (3) rainfall and grass cover for all the transects [to

ered a proxy for the stable state of the system (e.g., Hirota et al.,

determine if there was a relationship between rainfall and vegetation].

2011). This criterion allowed us to check whether the indicators

We computed the 95% confidence interval range of the correlation

behaved as expected, according to the theory, as the system moves

coefficient by applying the bootstrap method 1000 times [i.e., we ran-

away from its stable state. We then investigated trends of spatial indi-

domly resampled the rainfall and vegetation datasets with replacement

cators as a function of (1) rainfall gradient and (2) distance along the

and computed the correlation coefficient]. If the correlation coefficient

transect. We quantified trends of spatial indicators by slopes, estimated

for the original dataset was not biased due to sampling it should fall

by fitting a linear regression model for both the above relationships.

within the 95% confidence interval of the correlation coefficients

Positive [negative] slopes indicate increasing [decreasing] trends of

obtained by the bootstrap method.

indicators. The magnitude of the slope suggests the strength of the trend. We also obtained 95% confidence intervals for the estimation of

2.3 | Spatial indicators 2.3.1 | Null model and control for spatial indicators

slope by the bootstrap method [see the section Rainfall–vegetation Analysis]. Along transects, intermediate regions between grasslands and

We devised a null model, equivalent to random shuffling of elements in

woodlands exhibit two types of patterns. One type, as seen for trans-

the vegetation matrix, to compare whether spatial indicators computed

ects 1–5 and 9, has intermediate savanna regions containing multiple

from the transect data were different from those computed from a ran-

patches of woodland interspersed with grassland patches. The other

dom spatial pattern [see Appendix S3] (Dale & Fortin, 2014). We

type, as seen for transects 6, 7 and 8, does not have this intermediate

extracted vegetation data for three ‘control transects’ [7.5 km wide]

region, but has a distinct edge separating grasslands and woodlands

from areas where the range of a potential driver of vegetation [rainfall]

[Figure 1]. We chose four of the transects [4, 5, 8 and 9] to demon-

was relatively constant [500–650 mm year21; Appendix S5]. This

strate the analyses and results in the main text. We also conducted

6

|

EBY

ET AL.

F I G U R E 3 The frequency of occurrence of different proportions of grass cover for different ranges of rainfall over the entire Serengeti– Mara ecosystem landscape. Smooth curves represent the estimated kernel density. These results indicate the possibility of grasslands and woodlands being alternative stable states with rainfall as one potential driver

sensitivity analyses to test the robustness of our results to the size at

key driver of transitions from grasslands to woodlands in the Seren-

which we coarse-grained the data, the width of transects and the low

geti–Mara ecosystem.

frequency that we selected for DFT [Appendix S5].

3 | RESULTS

3.2 | Spatial indicators For transects 4, 5, 8 and 9, spatial variance, spatial skewness, spatial

3.1 | Rainfall–vegetation relationship For regions of the Serengeti–Mara with lower ranges of rainfall [500– 650 mm year21] the frequency distribution of grass cover was unimodal with the mode occurring at nearly full grass cover [Figure 3a] [pvalue of DIP test 5 .68 suggesting that the unimodality hypothesis of data cannot be rejected; estimated kernel density mode at 99.7% grass cover]. For regions that experienced higher rainfall [950–1100 mm

correlation at lag-1 and spatial DFT at low frequencies increased as the system changed from a grassland state to a woodland state as a function of rainfall [Figure 5]. In transect 5, for example, as the rainfall increased from 600 to 720 mm year21, the average grass cover only decreased by about 15% before beginning a steep decline to lower values of grass cover corresponding to woodlands. However, over the same change in rainfall, spatial variance showed a larger change,

] the frequency distribution of grass cover revealed a bimodal

increasing from small values less than 0.001 to 0.1 [Figure 5b1 and b2].

distribution [rejection of unimodality by DIP test with p-value < 1026;

Changes in the other indicators were less dramatic, but still substantial:

“flexmix” suggests a mixture model] with modes occurring at 2% grass

Spatial skewness increased from 210 to values close to zero; spatial

cover [corresponding to woodlands] and at 80% grass cover

correlation at lag-1 increased from 0.2 to 0.6; and the spatial low-

[Figure 3d]. At intermediate levels of rainfall [650–800 mm and 800–

frequency spectra increased from 6 to 30 [Figure 5b3–b5]. The slopes

21

year

], the statistical test rejected unimodality [p-values of

for spatial indicators for three transects [5, 8 and 9] were positive and

DIP tests < .001], “flexmix” suggested a mixture model with kernel den-

statistically significant, indicating increasing trends. These slopes repre-

sity estimating two [10% and 95.7% grass cover] and four [4%, 27%,

sent the response of indicators to a 1-mm increase in mean annual

66% and 92% grass cover] modes, respectively.

rainfall, which is why the magnitudes of the slopes for the transects

21

950 mm year

At the scale of the transects, we found that grass cover and rainfall

were small. Nevertheless, these slopes were at least an order of magni-

changed considerably on moving along the transects [see Figure 4a–d

tude larger than those for the null data. Furthermore, the confidence

for the four representative transects and Figure S4 for other transects;

intervals for the slopes of the transect data and the null data do not

see Table S3 for details of statistical correlations]. For each of these

overlap, with the exception of spatial skewness for transect 9 [Table

transects, we found that the percentage of grass cover declined as a

1a]. We did not quantify trends for transects 1, 2, 3, 6 and 7 [Figure

function of increasing rainfall [see Figure 5a1–d1 and Figures S5 and

S6] because they had fewer than six data points before declining to

Table S3 for all transect statistics], suggesting that rainfall could be a

85% grass cover.

EBY

|

ET AL.

7

Grass cover Rainfall

Transect-4

(B)

1.0

600 1.0

0.5

550 0.5

0.0 0

35

500 0.0 70 0

Transect-5 Grass cover

900 750

40

80

(C)

Transect-8

1.0

800 1.0

0.5

700 0.5

600 0.0

0

20

40

(D)

600 0.0 0

Transect-9

800 700 25

50

600

Rainfall (mm)

Grass cover

(A)

Distance along the transect (km) F I G U R E 4 Gradient of mean annual rainfall [blue squares] and grass cover [grey circles] along transects 4, 5, 8 and 9. Graphs for the other transects are shown in Fig. S4 [Color figure can be viewed at wileyonlinelibrary.com]

We also investigated trends of spatial indicators as a function

vegetation state; additionally, the rainfall control transects also

of distance along the transect: We found that slopes of spatial indi-

lacked any clear trends for spatial indicators. For one transect, num-

cators for six transects [numbered 1, 2, 5, 6, 8 and 9] were all posi-

bered 4, slopes were not statistically significant [Figure S9, Table

tive, indicating increasing trends [Table 1b]. In contrast, we found

1B] and the other two transects [3 and 7] had too few data points

that the rainfall control transects showed no [or little] change in

to compute trends statistically.

(a1)–(d1) Grass cover changes considerably as a function of mean annual rainfall for transects 4, 5, 8 and 9, suggesting rainfall to be a driver of transition between grassland and woodland states. The remaining graphs illustrate that CSD-based spatial indicators show increasing trends along gradients of high grass cover to lower cover, consistent with the theoretical predictions of Figure 2. All indicators were computed using sliding vegetation matrices along transects [blue circles] and corresponding null models [black squares]. The error bars on our null model data are negligible on the scale of our graphs owing to the fact that we have a relatively large number of data points. [See Table 1(a) for comparison of strengths of trends of vegetation data with null data. See Figures S5 and S6 for graphs of the other transects] [Color figure can be viewed at wileyonlinelibrary.com]

FIGURE 5

8

|

EBY

ET AL.

The sensitivity analyses of the size at which data were coarse-

states and abrupt regime shifts between those states. However, the

grained [Figure S10], the width of transects [Figure S11] and the low

other five transects showed linear or gradual changes in grass cover

frequency range we selected for DFT calculations [Figure S12], all

with rainfall [Figure S5], suggesting that local-scale patterns may differ

showed that the qualitative nature of our results were robust.

from the expectations based on analysis at the larger landscape scale. Our study showed that stable states in this system occur in areas with

4 | DISCUSSION

nearly complete grass cover and areas with almost no grass cover [Figure 3a]. Therefore, areas with intermediate cover at the interface of

Our analyses of high-resolution spatial data from the Serengeti–Mara

grasslands and woodlands are likely to be near the threshold/critical

ecosystem revealed the existence of alternative stable states of grass-

point of a regime shift. The theory of critical slowing down predicts

lands and woodlands. We showed that mean annual rainfall could be

that a system close to its critical point takes longer to return to equilib-

one of the drivers of transitions between these alternative stable states

rium (van Nes & Scheffer, 2007; Wissel, 1984), i.e. they are highly sen-

along spatial gradients [Figures 4, S4 and S5]. We found that the CSD-

sitive to perturbations. Indeed, we expect that stochastic variations in

based spatial metrics of spatial variance, spatial correlation at lag-1,

drivers, or disturbances such as fire and herbivory, are likely to cause

spatial skewness and spatial DFT at low frequencies, typically increased

frequent perturbations to savanna ecosystems, shifting them away

along spatial gradients of transitions from grasslands to woodlands

from their equilibrium state. Consequently, these intermediate areas

[Figures 4, S6 and S7]. These results, when interpreted within the

may undergo regime shifts on small spatial scales. Therefore, we specu-

space-for-time substitution framework [Box 1], suggest that CSD-

late that areas with intermediate grass cover are likely to be away from

based spatial indicators can offer signatures that an ecological system

their equilibrium state, thus masking the expected nonlinear relation-

is approaching a critical transition.

ship between rainfall and grass cover. An analysis of changes in vegetation along transects over time may help test and tease apart some of

4.1 | Regime shifts and spatial scales Our analysis, at the level of the entire Serengeti–Mara landscape [approximately 25,000 km2], showed that the vegetation type changes from an unimodal distribution [of high grass cover] to a bimodal distri-

these possibilities.

4.2 | Spatial indicators of critical slowing down Spatial indicators of regime shifts can be categorized into “generic spa-

bution [of high and low grass cover] as a function of mean annual rain-

!fi et al., 2014). Previous tial indicators” and “patch-based indicators” (Ke

fall [Figure 3]. The occurrence of multiple modes is an indicator that

studies using field data from semi-arid ecosystems have focused on

the ecosystem exhibits alternative stable states (Hirota et al., 2011;

patch-based indicators, such as the probability distribution of vegeta-

Scheffer et al., 2001; Staver et al., 2011). Such systems may undergo

!fi et al., 2007; Lin, Han, Zhao, & tion patches of different sizes (Ke

abrupt shifts from one state to another, either when the system

Chang, 2010). However, early warning signals based on patch-based

approaches the threshold point (Hirota et al., 2011; Scheffer et al.,

metrics could be sensitive to the details of local interactions, such as

2001) or due to environmental stochasticity even when the system is

the scale and strength of positive feedbacks in relation to negative

away from the threshold point (Beaugrand, Edwards, Brander, Luczak,

!fi et al., 2014; Weerman et al., 2012). Here, we instead feedbacks (Ke

& Ibanez, 2008; Guttal & Jayaprakash, 2007a). Theoretical studies,

used generic spatial indicators that have been shown to typically

which incorporate spatial interactions, typically assume that the driver

depend only on the proximity to a regime shift not on the details of

or stressor of the ecosystem is spatially homogeneous and predict that

the processes or mechanisms controlling an ecological system

the spatial boundary of alternative states is unstable (van de Leemput,

!fi, Rietkerk, Van Nes, & Scheffer, (Carpenter & Brock, 2010; Dakos, Ke

van Nes, & Scheffer, 2015). However, our results showed that the

!fi et al., 2014). Three indicators, 2011; Guttal & Jayaprakash, 2009; Ke

Serengeti–Mara ecosystem exhibits a driver with a spatial gradient

spatial variance, correlation at lag-1 and DFT at low frequencies,

[rainfall] along with a gradient of vegetation types [grassland to wood-

showed increasing trends from small positive values as we moved

land]. Our spatially explicit model [Box 1] offers a framework to investi-

away from the stable full grass cover mode for most of the transects

gate dynamics of regime shifts in ecosystems with a spatial gradient of

[Figures 5, S6 and S7, Table 1]. Based on the theory of critical slowing

drivers/stressors. This model showed that alternative vegetation states

down, we hypothesize that these trends in spatial indicators are signa-

can exhibit stable boundaries due to the existence of a driver that

tures of increasing sensitivity of the system to perturbations as the sys-

exhibits a spatial gradient [Figure S1]. Therefore, it is likely that the

tem moves away from one stable state towards an alternative state.

rainfall gradient is stabilizing the spatial gradient of grassland to wood-

Spatial skewness also showed increasing trends, but from negative val-

land states at the Serengeti–Mara landscape scale. However, as we dis-

ues to zero; areas with zero skewness are likely to consist of a mosaic

cuss below, abrupt shifts from one vegetation state to the other may

of grassland and woodland states, thus announcing the approach to

occur near the boundaries of these alternative states.

the alternative state. However, transect 4 in Figure S7 showed no sta-

At the smaller transect scale, four of the nine transects [4, 5 and 8

tistical signal of indicators along the transect distance [Table 1]. Here,

shown in Figure 4 and 6 in Figure S5] exhibited a threshold-like rela-

the system is moving away from the stable state [nearly full grass

tionship between grass cover and rainfall, consistent with the above

cover] in a non-monotonic fashion along the transect distance. As the

hypothesis that this savanna ecosystem exhibits alternative stable

system moves away from [or moves back to] the stable state, spatial

0.0006

4.5 3 10

0.0044

2.42 3 10

8 (vegetation)

8 (null)

9 (vegetation)

9 (null)

25

]

, 1.0 3 10

24

]

[0.003, 0.013]

[0.004, 0.045]

[0.047, 0.08]

–

[0.01, 0.06]

–

–

0.004

0.004

0.03

–

0.01

0.004

20.0001

20.0005

2.1e-04

3

4

5

6

7

8

9

Rain control A

Rain control B

Rain control C

20.04 20.157

[–0.001, 24.1 3 1025

[–6.0 3 1025, 0.0004]

[–1.08, 0.697]

[–0.96, 1.24

[–0.46, 0.71]

[0.04, 0.4]

[0.7, 3.21]

–

[1.44, 4.47]

[0.37, 2.72]

[–0.348, 0.751]

–

[1.578, 3.77]

[0.088, 0.373]

26

, 3.5 3 10

[–5.7 3 10

24

, 3.5 3 10

[0.0005, 0.006]]

[–5.7 3 10

24

[0.004, 0.009]

[–0.002, 0.0001]

[0.001, 0.005]

95% CI

0.01

20.02

20.002

0.028

0.094

–

0.135

0.027

20.02

–

0.11

0.06

Slope

[–0.009, 0.03]

[–0.05, 0.0002]

[–0.01, 0.006]

[0.004, 0.5]

[0.076, 0.124]

–

[0.07, 0.2]

[0.008, 0.04]

[–0.06, 0.018]

–

[0.066, 0.157]

[0.035, 0.075]

95% CI

Spatial autocorrelation at lag-1

25.4 3 10

0.003

25

25

6.1 3 10

0.006

25.4 3 10

0.003

Slope

Spatial autocorrelation at lag-1

24

24

]

]

28

28

[0.003, 0.02]

[0.06, 0.42]

[–1.4 3 1027, 1.4 3 1029]

[0.187, 0.538]

[–5.8 3 1028, 24.1 3 1029]

[0.13, 0.49]

95% CI

26

21.1 3 1025

6.6 3 1025

7.7 3 10

1.94

5.18

–

11.93,

1.99

2.105

–

7.415

3.256

Slope

[–0.0001, 0.0001]

[–6.9 3 1025, 0.0001]

[–1.1 3 1025, 2.6 3 1025]

[0.56, 3.51]

[3.1, 7.21]

–

[6.82, 17.1]

[1.277, 2.79]

[–1.07, 4.32]

–

[1.69, 12.98]

[1.416, 4.712]

95% CI

Spatial DFT at low-frequencies

0.01

0.21

23.8 3 10

0.333

22.3 3 10

0.21

Slope

Spatial DFT at low frequencies

The confidence intervals (CI s) for the slopes of the transect data and the null data do not overlap, with the exception of spatial skewness for transect 9, suggesting statistical significance of trends. Rainfall control transects lacked any clear trends for spatial indicators (Figures S8 and S9). For transect 4, slopes were not statistically significant (Figure S9) and the other two transects (3 and 7) had too few data points to compute trends statistically.

0.0449

0.236

1.59

–

2.86

1.66

0.22

–

2.512

0.249

]

[–0.0002, 7.4 3 10

[0.0004, 0.008]

[0.004, 0.016]

[0.0031, 0.0044]

[–0.006, 0.009]

[0.002, 0.039]

0.021

2

[0.002, 0.01]

0.007

95% CI

Slope

0.009

0.026

0.054

[0.044, 0.22]

Spatial skewness

26

, 5.1 3 10

25

]

0.105

[0.004, 0.01]

[0.03, 0.36]

Slope

95% CI

[3.95 3 10

26

[5.4 3 1025,1.0 3 1023]

[1.3 3 10

25

0.008

0.14

Spatial variance

25

, 2.8 3 10

[0.0002, 0.001]

[1.5 3 10

25

[0.0003, 0.0008]

1

Transect

(b)

2.2 3 10

5 (null)

25

0.0004

25

Slope

Slope

95% CI

Spatial skewness

Spatial variance

5 (vegetation)

Transect

(a)

95% CI

(a) Statistics for trends for spatial indicators, quantified by slope, as a function of rainfall for transects 5, 8 and 9 [other transects had too few data points for statistical analysis; also see Fig S6]. (b) Statistics for trends for spatial indicators, quantified by slope, as a function of distance showing increasing trends for most transects [see also Fig S7]

T A BL E 1

EBY ET AL.

| 9

10

|

EBY

ET AL.

indicators are expected to show increasing [decreasing] trends. Thus,

AC KNOW LEDG EME NT S

our theoretical expectation is that indicators too show non-monotonic

V.G. acknowledges research support from a DBT Ramalingaswamy

trends in such transects.

Fellowship, DBT-IISc partnership program, ISRO-IISc Space Technol-

We focused on rainfall as one of the possible drivers of vegeta-

ogy Cell and infrastructure support from DST-FIST. V.G. also

tion states in the Serengeti–Mara ecosystem. Our analyses do not,

acknowledges facilities at the Simons Centre for the Study of Living

however, preclude the role of other drivers and more complex eco-

Machines at NCBS-TIFR, where part of this work was carried out. A.

logical mechanisms in shaping the observed patterns of vegetation

P.D. and S.E. were supported by funds from the McDonnell Founda-

and spatial metrics. Indeed, various factors such as herbivory [involv-

tion. We thank Michael Coughenour for creating the rainfall data.

ing both migratory and non-migratory organisms], fire and soil char-

We thank the anonymous referees, Krishnapriya Tamma, Hari Srid-

acteristics also shape spatial vegetation patterns in savanna

har and Sumithra Sankaran and members of the Nuffield Sponsored

ecosystems at different scales (Asner, Elmore, Olander, Martin, &

Serengeti Fire Project for substantial comments on the manuscript.

Harri, 2004; Lehmann et al., 2014; Owen-Smith, 1989; Sankaran, Ratnam, & Hanan, 2004; Sankaran et al., 2005; Staver, Bond, Stock, Van Rensburg, & Waldram, 2009; Tiessen, Feller, Sampaio, & Garin, 1998; van Langevelde et al., 2003). Although mathematical models show that generic spatial indicators are insensitive to various !fi et al., 2014), theoretical and empirical system-specific features (Ke work are required to understand the efficacy of leading indicators in different types of spatial systems such as those exhibiting Turing patterns, Maxwell points or systems driven by heterogeneous, seasonal and stochastic stressors (Guttal & Jayaprakash, 2007b; Riet!fi, 2016; van de Leemput et al., kerk et al., 2004; Schneider & Ke 2015). Additionally, studies should focus on investigating other drivers of vegetation dynamics, conducting detailed investigations of spatial heterogeneity in habitat characteristics like soil and topography across the entire landscape and the role of disturbances such as fire, herbivory and anthropogenic land use, and the effect of all of these on early warning signals. Investigations of spatial patterns of vegetation over time and manipulative experiments in the field to measure the dynamics causing “slow down” would be another

AUTHOR CONTRIBUTI ONS This project began during conversations between Stephanie Eby, Andrew P. Dobson and Vishwesha Guttal on linking the theoretical models of alternative stable states and state changes with landscape-scale data from savanna ecosystems. S.E., A.P.D. and V.G. designed the study; S.E., A.A., S.M. and V.G. provided methods; A.A., S.M., S.E. and V.G. analysed results; A.A., S.M. and S.E. created the figures; S.E. and V.G. wrote the first draft of the manuscript; all authors contributed substantially to revisions. S.E. and A.A. contributed equally to the manuscript. DATA ACCESS IBILI TY Simulation codes to reproduce model results, spatial-indicator codes and vegetation data for one transect are available at GitHub (https://github.com/tee-lab/spacetime-csd). Full vegetation data are available from Reed et al. (2009). Rainfall data may be obtained from M. Coughenour.

important direction for future research (Verbesselt et al., 2016). RE FE RE NCE S

5 | CONCLUDING REMARKS Our study has revealed that simple spatial metrics such as spatial variance, spatial skewness, spatial autocorrelation at lag-1 and spatial DFT at low frequencies show theoretically expected changes along alternative vegetation states in the Serengeti–Mara ecosystem. Therefore, we argue that CSD-based spatial metrics are valid and reliable substitutes for time-series-based metrics to detect imminent critical transitions in large-scale ecosystems, since CSD-based indicators arise due to the sensitivity of ecosystems to perturbations (Scheffer et al., 2009). Furthermore, we suggest that regions of ecosystems that show strong signals of CSD-based spatial metrics may be “sensitive regions” that are prone to abrupt regime shifts. Using CSD-based spatial metrics to identify “sensitive regions” could help the conservation and management of areas of ecosystems prone to abrupt regime shifts. Future studies should combine time series of spatial data in various ecosystems to determine the strengths and weaknesses of CSD-based spatial metrics as forecasting tools of regime shifts.

Asner, G. P., Elmore, A. J., Olander, L. P., Martin, R. E., & Harri, A. T. (2004). Grazing systems, ecosystem responses, and global change. Annual Review Environmental Resources, 29, 261–299. Bagchi, S., Briske, D. D., Wu, X. B., McClaran, M. P., Bestelmeyer, B. T., !nez, M. E. (2012). Empirical assessment of state& Fern!andez-Gime and-transition models with a long-term vegetation record from the Sonoran Desert. Ecological Applications, 22, 400–411. Blois, J. L., Williams, J. W., Fitzpatrick, M. C., Jackson, S. T., & Ferrier, S. (2013). Space can substitute for time in predicting climate-change effects on biodiversity. Proceedings of the National Academy of Sciences USA, 110, 9374–9379. Browning, D. M., Franklin, J., Archer, S. R., Gillan, J. K., & Guertin, D. P. (2014). Spatial patterns of grassland-shrubland state transitions: A 74-year record on grazed and protected areas. Ecological Applications, 24, 1421–1433. Beaugrand, G., Edwards, M., Brander, K., Luczak, C., & Ibanez, F. (2008). Causes and projections of abrupt climate-driven ecosystem shifts in the North Atlantic. Ecology Letters, 11, 1157–1168. Carpenter, S. R., & Brock, W. A. (2006). Rising variance: A leading indicator of ecological transition. Ecology Letters, 9, 311–318. Carpenter, S. R., & Brock, W. A. (2010). Early warnings of regime shifts in spatial dynamics using the discrete Fourier transform. Ecosphere, 1, 1–15.

EBY

ET AL.

Carpenter, S. R., Cole, J. J., Pace, M. L., Batt, R., Brock, W. A., Cline, T., . . . Weidel, D. (2011). Early warnings of regime shifts: A wholeecosystem experiment. Science, 332, 1079–1082. Carpenter, S. R., Ludwig, D., & Brock, W. A. (1999). Management of eurtrophication for lakes subject to potentially irreversible change. Ecological Applications, 9, 751–771. Cline, T. J., Seekell, D. A., Carpenter, S. R., Pace, M. L., Hodgson, J. R., Kitchell, J. F., & Weidel, B. C. (2014). Early warnings of regime shifts: Evaluation of spatial indicators from a whole-ecosystem experiment. Ecosphere, 5, 1–13.

|

11

!fi, S., Rietkerk, M., Alados, C. L., Pueyo, Y., Papanastasis, V. P., Elaich, Ke A., & De Ruiter, P. C. (2007). Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems. Nature, 449, 213–217. Lehmann, C. E., Anderson, T. M., Sankaran, M., Higgins, S. I., Archibald, S., Hoffmann, W. A., . . . (2014). Savanna vegetation–fire–climate relationships differ among continents. Science, 343, 548–552. Lenton, T. M. (2011). Early warning of climate tipping points. Nature Climate Change, 1, 201–209.

Dai, L., Korolev, K. S., & Gore, J. (2013). Slower recovery in space before collapse of connected populations. Nature, 496, 355–358.

Lin, Y., Han, G., Zhao, M., & Chang, S. X. (2010). Spatial vegetation patterns as early signs of desertification: A case study of a desert steppe in Inner Mongolia, China. Landscape Ecology, 25, 1519–1527.

Dai, L., Vorselen, D., Korolev, K. S., & Gore, J. (2012). Generic indicators for loss of resilience before a tipping point leading to population collapse. Science, 336, 1175–1177.

Litzow, M. A., Urban, J. D., & Laurel, B. J. (2008). Increased spatial variance accompanies reorganization of two continental shelf ecosystems. Ecological Applications, 18, 1331–1337.

Dakos, V., Carpenter, S. R., Brock, W. A., Ellison, A. M., Guttal, V., Ives, A. R., . . . Scheffer, M. (2012). Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PLoS One, 7, e41010.

Lubeck, S. (2006). Tricritical directed percolation. Journal of Statistical Physics, 123, 193–221.

!fi, S., Rietkerk, M., Van Nes, E. H., & Scheffer, M. (2011). Dakos, V., Ke Slowing down in spatially patterned ecosystems at the brink of collapse. The American Naturalist, 177, E153–E166. Dakos, V., Nes, E. H., Van Donangelo, R., Fort, H., & Scheffer, M. (2010). Spatial correlation as leading indicator of catastrophic shifts. Theoretical Ecology, 3, 163–174. Dakos, V., Scheffer, M., Nes, E. H., Van Brovkin, V., Petoukhov, V., & Held, H. (2008). Slowing down as an early warning signal for abrupt climate change. Proceedings of the National Academy of Sciences USA, 105, 14308–14312. Dale, M. R., & Fortin, M. J. (2014). Spatial analysis: A guide for ecologists (2nd edn). Cambridge: Cambridge University Press. Drake, J. M., & Griffen, B. D. (2010). Early warning signals of extinction in deteriorating environments. Nature, 467, 456–459. Dublin, H. T. (1991). Dynamics of the Serengeti–Mara woodlands: An historical perspective. Forest and Conservation History, 35, 169–178. Favier, C., Aleman, J., Bremond, L., Dubois, M. A., Freycon, V., & Yangakola, J. M. (2012). Abrupt shifts in African savanna tree cover along a climatic gradient. Global Ecology and Biogeography, 21, 787–797. Guttal, V., & Jayaprakash, C. (2007a). Impact of noise on bistable ecological systems. Ecological Modelling, 201, 420–428. Guttal, V., & Jayaprakash, C. (2007b). Self-organization and productivity in semi-arid ecosystems: Implications of seasonality in rainfall. Journal of Theoretical Biology, 248, 490–500. Guttal, V., & Jayaprakash, C. (2008). Changing skewness: An early warning signal of regime shifts in ecosystems. Ecology Letters, 11, 450–460.

Norton-Griffths, M., Herlocker, D., & Pennycuick, L. (1975). The patterns of rainfall in the Serengeti ecosystem, Tanzania. African Journal of Ecology, 13, 347–374. Owen-Smith, N. (1989). Megafaunal extinctions: The conservation message from 11,000 years B.P. Conservation Biology, 3, 405–412. Ratajczak, Z., Nippert, J. B., & Ocheltree, T. W. (2014). Abrupt transition of mesic grassland to shrubland: Evidence for thresholds, alternative attractors, and regime shifts. Ecology, 95, 2633–2645. Reed, D. N., Anderson, T. M., Dempewolf, J., Metzger, K., & Serneels, S. (2009). The spatial distribution of vegetation types in the Serengeti ecosystem: The influence of rainfall and topographic relief on vegetation patch characteristics. Journal of Biogeography, 36, 770–782. Rietkerk, M., Dekker, S. C., de Ruiter, P. C., & van de Koppel, J. (2004). Self-organized patchiness and catastrophic regime shifts in ecosystems. Science, 305, 1926–1929. Sankaran, M., Hanan, N. P., Scholes, R. J., Ratnam, J., Augustine, D. J., Cade, B. S., . . . Ardo, J. (2005). Determinants of woody cover in African savannas. Nature, 438, 846–849. Sankaran, M., Ratnam, J., & Hanan, N. P. (2004). Tree-grass coexistence in savannas revisited – insights from an examination of assumptions and mechanisms invoked in existing models. Ecology Letters, 7, 480–490. Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., . . . Sugihara, G. (2009). Early-warning signals for critical transitions. Nature, 461, 53–59. Scheffer, M., Carpenter, S., Foley, J. A., Folke, C., & Walker, B. (2001). Catastrophic shifts in ecosystems. Nature, 413, 591–596. !fi, S. (2016). Spatially heterogeneous pressure Schneider, F. D., & Ke raises risk of catastrophic shifts. Theoretical Ecology, 9, 207–217.

Guttal, V., & Jayaprakash, C. (2009). Spatial variance and spatial skewness: Leading indicators of regime shifts in spatial ecological systems. Theoretical Ecology, 2, 3–12.

Scholes, R. J., & Archer, S. R. (1997). Tree–grass interactions in savannas. Annual Review of Ecological Systems, 28, 517–544.

Held, H., & Kleinen, T. (2004). Detection of climate system bifurcations by degenerate fingerprinting. Geophysical Research Letters, 31, L23207.

Sinclair, A. R. E. (1979). The Serengeti environment. In A. R. E. Sinclair & M. Norton-Griffiths (Eds.) Serengeti: Dynamics of an ecosystem (pp. 31–45). Chicago, IL: University of Chicago Press.

Hirota, M., Holmgren, M., Van Nes, E. H., & Scheffer, M. (2011). Global resilience of tropical forest and savanna to critical transitions. Science, 334, 232–235. Johnson, E. A., & Miyanishi, K. (2008). Testing the assumptions of chronosequences in succession. Ecology Letters, 11, 419–431.

Sinclair, A. R. E., Hopcraft, J. G. C., Olff, H., Mduma, S. A. R., Galvin, K. A., & Sharam, G. J. (2007). Historical and future changes to the Serengeti ecosystem. In A. R. E. Sinclair, C. Packer, S. Mduma, & J. Fryxell (Eds.), Serengeti III: Biodiversity and biocomplexity in a human-influenced ecosystem (pp. 7–46). Chicago, IL: University of Chicago Press.

!fi, S., Guttal, V., Brock, W. A., Carpenter, S. R., Ellison, A. M., Livina, V. Ke N., . . . Dakos, V. (2014). Early warning signals of ecological transitions: Methods for spatial patterns. PLoS One, 9, e92097.

Staver, A. C., Archibald, S., & Levin, S. A. (2011). The global extent and determinants of savanna and forest as alternative biome states. Science, 334, 230–232.

12

|

Staver, A. C., Bond, W. J., Stock, W. D., Van Rensburg, S. J., & Waldram, M. S. (2009). Browsing and fire interact to suppress tree density in an African savanna. Ecological Applications, 19, 1909–1919. Steele, J. H. (1998). Regime shifts in marine ecosystems. Ecological Applications, 8, S33–S36. Streeter, R., & Dugmore, A. J. (2013). Anticipating land surface change. Proceedings of the National Academy of Sciences USA, 110, 5779–5784.

EBY

ET AL.

patch-size distribution and degradation in a spatially self-organized intertidal mudflat ecosystem. Ecology, 93, 608–618. Wissel, C. (1984). A universal law of the characteristic return time near thresholds. Oecologia, 65, 101–107.

BIOSK ET CH

Tiessen, H., Feller, C., Sampaio, E. V. S. B., & Garin, P. (1998). Carbon sequestration and turnover in semiarid savannas and dry forest. Climatic Change, 40, 105–117.

STEPHANIE EBY is a lecturer in the Department of Marine and Environ-

van de Leemput, I. A., van Nes, E. H., & Scheffer, M. (2015). Resilience of alternative states in spatially extended ecosystems. PLoS One, 10, e0116859.

cifically, she is interested in understanding how burning impacts the

van Langevelde, F., van de Vijver, C. A., Kumar, L., van de Koppel, J., de Ridder, N., van Andel, J., . . . Prins, H. H. (2003). Effects of fire and herbivory on the stability of savanna ecosystems. Ecology, 84, 337–350. van Nes, E. H., & Scheffer, M. (2007). Slow recovery from perturbations as a generic indicator of a nearby catastrophic shift. The American Naturalist, 169, 738–747. Veraart, A. J., Faassen, E. J., Dakos, V., van Nes, E. H., Lurling, M., & Scheffer, M. (2012). Recovery rates reflect distance to a tipping point in a living system. Nature, 481, 357–359.

mental Sciences at Northeastern University. Her research focuses on how species influence and are influenced by their environments. Spebehaviour and distribution of large mammals and how burning and grazing impact plant species composition.

SUPPORTING INFORMATION Additional Supporting Information may be found online in the supporting information tab for this article:

How to cite this article: Eby S, Agrawal A, Majumder S, Dobson

Verbesselt, J., Umlauf, N., Hirota, M., Holmgren, M., van Nes, E. H., Herold, M., . . . Scheffer, M. (2016). Remotely sensed resilience of tropical forests. Nature Climate Change, 6, 1028–1031.

AP, Guttal V. Alternative stable states and spatial indicators of

!fi, S., Weerman, E. J., van Belzen, J., Rietkerk, M., Temmerman, S., Ke Herman, P. M. J., & van de Koppel, J. (2012). Changes in diatom

10.1111/geb.12570

critical slowing down along a spatial gradient in a savanna ecosystem. Global Ecol Biogeogr. 2017;00:1–12. https://doi.org/

Online Supporting Information: Alternative stable states and spatial indicators of critical slowing down along a rainfall gradient in a savanna ecosystem Stephanie Ebya,b,*,+, Amit Agrawalc,+, Sabiha Majumderd, Andrew P. Dobsona,d, Vishwesha Guttalc,* a

Department of Ecology and Evolutionary Biology, Princeton University, Princeton, NJ

08544, USA b

Current address: Department of Marine and Environmental Sciences, Northeastern

University, Boston, MA 02115, USA c

Centre for Ecological Sciences, Indian Institute of Science, Bengaluru, 560012, India

d

Department of Physics, Indian Institute of Science, Bengaluru, 560012, India

d

Santa Fe Institute, Santa Fe, NM 87501, USA

+These two authors contributed equally to the manuscript.

Table of contents: Appendix S1: Spatially explicit model and steady-state phase diagram Fig. S1: Steady state phase diagram of the model in Appendix S1 Appendix S2: Vegetation and rainfall data Table S1: Vegetation classification scheme Table S2: Spatial scales of analyses and terminology used Fig. S2: Close up images of Transects five and eight Appendix S3: Null model and coarse-graining (aggregating) the spatial data: Analytical results for spatial variance and spatial skewness Fig. S3: An illustration of the aggregating procedure Appendix S4: Results for other transects Fig. S4: Grass cover and rainfall over distance for each transect Fig. S5: Grass cover as a function of rainfall for each transect Table S3: Statistics for correlation coefficients for all transects Fig. S6: Spatial indicators over rainfall for all transects Fig. S7: Spatial indicators over distance for all transects Appendix S5: Spatial indicators for rainfall control transects Fig. S8: Rainfall map and location of control transects Fig. S9: Trends of indicators for rainfall control transects Appendix S6: Sensitivity analyses Fig. S10: Sensitivity of spatial indicators to size of coarse-graining Fig. S11: Sensitivity of spatial indicators to width of transects Fig. S12: Sensitivity of low-frequency DFT results

Appendix S1: Detailed description of modelling the spatial gradient of a stressor Here, we describe how we model spatial gradient of stressor in detail (model is presented in Box 1 of the main text). We represent a transect as a two dimensional horizontal rectangular grid consisting of L number of N x N matrices, i.e. the width of the transect is N cells and the length is NL. The N x N matrix at the left end of the transect has high birth rates and thus, is in a state of high cover. After every N cells along the length of this rectangular grid, the driver (p) decreases by dp; in other words, the death rate or the stressor increases from the left end of the transect towards the right end. The N x N matrix at the right end of this rectangular space has a low birth rate where the vegetation cover is zero at steady state. We initialize the grid by randomly distributing 1’s and 0’s across the entire transect. We assume that the boundary cells of a matrix can be occupied with a probability equal to the cover in the corresponding matrix. The value of the driver at which the system undergoes an abrupt transition, called the critical point and denoted by pc, was found to be at pc = 0.25. We choose the following parameter values for our simulations: N = 128 (comparable to the size of vegetation matrices); L =181 (comparable to the lengths of vegetation transects); pl = 0.258 (a value greater than the critical point); pr = 0.24 (a value less than the critical point); q = 0.92 (representing strong positive interactions)

Figure S1: Steady state phase diagram for the model presented in Box 1of the main text. The mean cover of the system exhibits a discontinuous transition, an abrupt regime shift, from high vegetation cover to a low cover state as a function of an increasing stressor value.

Appendix S2: Vegetation and Rainfall data S2.1 Vegetation classification, spatial scales and sample transect images Table S1 Original pixel classifications from the Reed et al. (2009) vegetation map and their new classifications as woodlands (0) or grasslands (1) Original Number 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Original Classification Unclassified Bare Water Cloud Agriculture or bare Kopje Sparse to open grassland Dense to closed (shrubbed) grassland Sparse to open shrubbed or treed grassland Dense to closed treed grassland Sparse to open grassed or shrubbed woodland Dense to closed (grassed) shrubland Sparse to open shrubbed woodland or treed shrubland Dense to closed treed shrubland Dense to closed (grassed or shrubbed) forest

New Number No Data No Data No Data No Data No Data No Data 1 1 1 1 0 1 0 0 0

New Classification No Data No Data No Data No Data No Data No Data Grassland Grassland Grassland Grassland Woodland Grassland Woodland Woodland Woodland

Table S2: Spatial scales of various components of analyses Terminology Vegetation pixel Coarse-grained scale Vegetation matrix Transect Rainfall-grid

Size 30 m x 30m 150 m x 150m 7.5 km x 7.5km 7.5 km x Length 2.5 km x 2.5km

(A)

(B)

!

!

Figure S2: Close up images of transects five and eight. The white pixels represent grasslands and the dark pixels represent woodlands.

S2.2 Rainfall data processing Rainfall data were collected at 40-50 rain gauges in and around the park by the Serengeti Ecological Monitoring Programme, the Serengeti National Park ecologist, Serengeti Wildlife Research Institute, and Tanzania Wildlife Conservation Monitoring. Gauges were read each month. The first data processing step was to search for instances when rain gauges were only read after 2 or 3 months instead of monthly. A computer program was developed to distribute the cumulative amount among the months with missing readings. The program used data from the nearest gauge with data for every month and the proportional distribution of rainfall among months from that gauge was used to distribute the cumulative amount for the gauge with missing readings. A program called PPTMAP (developed by M. Coughenour) was used to interpolate the rainfall data. The approach used in PPTMAP was originally developed as part of a spatially explicit ecosystem model (Coughenour 1992) and used inverse distance weighting. Here, however, Gaussian distance weighting was used, with the alpha parameter set to 6.25, as in Thorton et al. (1997). The 10 closest stations were included in the estimate for each point. Differences in elevation between a point estimate and rain gauge location were also corrected for in PPTMAP. This was done by running a regression for each rain gauge between rainfall and elevation using data from all gauges within a distance of 75 km. The significance of the regression was assessed using the F-test, and the regression on elevation was only applied if p < 0.1. Rain gauge data were corrected for significant differences in elevation between the gauge and the point being estimated. The Gaussian weighting was then carried out using the corrected data. PPTMAP searches for outlier data points using the Z statistic where a point was considered an outlier if the Z statistic was greater than 2.0 or less than -2.0. Outliers are discarded and the regression is recalculated, after which the slopes of significant regressions for each gauge are used in the distance weighting calculation. After examining the rainfall maps created as above it became clear that there were likely a variety of errors in the data due to under and overestimation. For example, there would be unusually large disagreements between adjacent gauges, and readings at one gauge were sometimes lower than from two or more of the surrounding gauges. A computer program called “ScreenData” was developed (also by M. Coughenour) to look for bad data points, discard them, and create a “clean” set of rain gauge data. Two tests were used to detect bad data points. First, the program looked for statistical outliers in the local neighborhood. The program searched for up to the 10 closest stations within a 50 km radius and calculated their mean and standard deviation. The deviation from the mean was calculated for each data point and then the z-score was calculated as the ratio of the deviation to the standard deviation. A data point was considered bad if the z-score was larger than 2.0 or smaller than 1.5. A less conservative criterion was used for positive deviations, because errors were more likely to result from under catch than over catch. A second test was applied based upon relationships between rainfall and elevation, as an elevation gradient could produce realistic deviations and such data points, even having failed the first test above, should not be discarded. As in PPTMAP, a linear regression was performed between rainfall and elevation for the neighborhood stations selected as above. The correlation coefficient and F-statistic were calculated for the regression. A regression was considered significant at the p < 0.1 level for F. When a regression was significant, the mean residual was calculated for the neighborhood points. When the residual for the data point in question was more than 2.5 times larger than the mean residual, it was considered a bad data point. The output from ScreenData was a clean data file that was then used in PPTMAP to create the final rainfall maps. PPTMAP was used to create monthly rainfall maps and then these were summed to create annual maps. Afterwards the annual maps were averaged to create the 30 year mean annual layer. Rainfall data were then extracted for each transect. References

Coughenour, M. B. (1992) Spatial modeling and landscape characterization of an African pastoral ecosystem: a prototype model and its potential use for monitoring drought. pp. 787810 in: D.H. McKenzie , D.E. Hyatt and V.J. McDonald (eds.). Ecological Indicators, Vol. I. Elsevier Applied Science, London and New York. Thornton, P.E., Running, S. W., White, M. A. (1997) Generating surfaces of daily meteorological variables over large regions of complex terrain. Journal of Hydrology 190: 214-251

Appendix S3: Null model and coarse-graining (aggregating) the spatial data: Analytical results for spatial variance and spatial skewness. S3.1 Spatial variance and spatial skewness for the unprocessed spatial data are independent of the spatial configuration. Randomization of spatial data, equivalent to bootstrapping, is a standard method in spatial statistics for estimating the significance of a spatial quantity (Fortin and Dale, 2005). Here, we show that such a method cannot work for spatial variance and spatial skewness. We first consider a spatial dataset (represented by a matrix of size !!×!) where each entry is either 1 (representing grassland) or 0 (representing woodland). We will refer to these unprocessed data as “raw vegetation data”. For these data, p denotes the proportion of the area covered by grassland. As discussed in the main test, we first created a null model matrix by creating a matrix of the same size as the raw vegetation data where any entry in the null matrix was assigned a value of 1 at a probability of p and otherwise assigned a value of 0. This null model matrix is equivalent to a random reshuffling of the entries of the raw vegetation data. For the above definition of a null matrix, the average proportion of grassland pixels will be the same as that of the raw spatial data (which is equal to p). The spatial variance (σ!! ) of such a matrix can be computed using the following: !

!!!

!

=! !!! !!!

(![!, !] − !)! !! !!

where z[i,j] is the value of the vegetation at matrix element at (i,j), ! is the spatial mean and N2 is the total number of entries in the matrix. For the random null matrix, we can analytically obtain the variance as follows: !!! = (!"#$#"%&#'!!"!1′!)(1 − !)! + ! (Proportion!of!0’s)(0 − p)! = !(1 − !)! + 1 − ! ! ! Simplifying the above equation, we get: !!! = !(1 − !) Likewise, it can be shown that the spatial skewness, denoted by!!! , is analytically given by: !

!

!! = ! !!! !!!

(![!, !] − !)! ! ! !!!!

!! = !

1 − 2! !!(1 − !)

In the derivation of the expressions for spatial variance and spatial skewness above, only the proportion of cover (p) was a relevant parameter but not the spatial configuration of the grassland (or woodland) entries. Therefore, the value of these indicators for the null model will be identical to those of the raw matrix. In other words, if we used the raw spatial data in the form of a spatial matrix containing 0’s and 1’s (representing wood and grass, respectively in our analyses), the above randomization method cannot provide a null model for computing the spatial variance and the spatial skewness. However, given that we have randomly reshuffled the data, the spatial correlation for this null matrix is expected to be zero and therefore, can act as a null model for spatial correlation. This is evident from figure S7 (see first, second and third rows of the first column). Below, we describe a method of ‘coarse-graining’ the data (referred to as ‘aggregating’ the data in rest of this section) that provides a null model for the analysis of spatial variance and skewness.

S3.2 Procedure for coarse-graining the spatial data The coarse-grained data (for both the vegetation matrix and the null matrix) were obtained by computing the average proportion of the entries covered with grassland for an !!×! sized non-overlapping submatrices of the raw data. This produced data with a lower resolution (!!×! where ! = !/!) that was composed of relatively ‘continuous’ variables ranging from 0-1 rather than the binary values of only 0 and 1 that make up the base map (figure S3).

Figure S3 An illustration of the aggregating procedure. The raw matrix on the left is of dimensions 6×6 (i.e., N = 6). The submatrix is of size 2!×!2 (i.e., a = 2). Two sample submatrices are shown in grey. The coarse-grained matrix, on the right, is of size 3 x 3 ! (! = = 3). !

Spatial variance and spatial skewness of coarse-grained null matrices will be smaller in comparison to raw null matrices In probability theory, Bernoulli trials are events where a ‘success’ occurs with a probability p (consequently, the probability of failure is 1-p). The probability that we will obtain k successful trials out of a total of b Bernoulli trials is given by a Binomial distribution (Bodine et al, 2014). Formally, for a binomially distributed random number, the probability that the number of successful events is equal to k in a total of b trials is given by: ! ! ! ! =! ! (1 − !)!!! ! Furthermore, we also know that (Bodine et al, 2014): !"#$ ! = !!" !"# ! = !!"(1 − !) 1 − 2! !"#$%#&& ! = ! !"(1 − !)

We can apply this to our raw and coarse-grained null matrices. The entries of raw null matrices are 1 with a probability p, and therefore each entry can be thought of as a Bernoulli trial. The entries of the coarse-grained null matrix, denoted by y, are an average of a2 number of raw null matrices. If the number of 1’s in the submatrix is equal to k, then k is a binomially distributed random number with ! = !" where ! = ! ! . Therefore, we can calculate statistics for the random variable y (i.e., statistics of the coarse-grained null matrix) by: ! = !! ! !!! ! ! 1−! !!! = !"# = = ! ! !! !! ! 1 − 2! !!! !! = !!"#$%#&&! =! =! ! ! ! !(1 − !) Comparing variance and skewness of the coarse-grained null matrix to the raw null matrix, we obtain: !!! !!! = ! ! ! !!! !! = ! ! It is clear from the above expressions, that whereas the mean of the coarse-grained null matrix remains unchanged from the raw null matrix, both the spatial variance and the spatial skewness reduce as we increase the size of the submatrix (!!×!!). !"#$ ! = !!"#$!

Spatial variance and spatial skewness of the coarse-grained vegetation data are not expected to reduce much due to patchiness Consider an extreme example of a landscape with high degree of patchiness of vegetation. The presence of grass (or woody vegetation) in one pixel will likely imply the presence of grass (or woody vegetation) in neighboring pixels as well. Consequently, most of the values of the submatrix will be similar (either 1 or 0), and therefore, the coarse-grained value (obtained by averaging the values of the submatrix) will be close to 1 or 0. If this is true for many of the non-overlapping submatrices, then most of the entries in the coarse-grained matrix will also be close to 1 or 0, thus retaining the binary structure of the raw vegetation matrix. Therefore, the spatial variance and skewness of such an coarse-grained matrix will not differ much from those computed for the raw vegetation matrix. In a more general scenario, where the patchiness exists but may not be very high, we argue that reduction in spatial variance and spatial skewness will not be proportional to the coarse-graining size. S3.3 Conclusion The above analyses show us that the spatial variance and spatial skewness of both the raw vegetation matrix and the raw null matrix are the same (from S2.1). Therefore, the null matrix obtained by random filling cannot be a null model for the spatial variance and skewness. However, the spatial variance and skewness of the coarse-grained vegetation matrix and the coarse-grained null matrix are affected different ways. The spatial variance and spatial skewness of the coarse-grained null matrix are both reduced by a factor proportional to the submatrix size. However, the spatial variance of the coarse-grained vegetation matrix, due to its spatial structure, is expected not to be reduced by a factor proportional to the submatrix size. Therefore, to check if spatial variance and skewness arises because of spatial structure, we should check whether the spatial variance of the coarse-grained vegetation matrix is more than that of the coarse-grained null matrix. In other words, the method of aggregating data as

described above provides a null model for the spatial variance and spatial skewness. The raw null model, however, works for spatial correlations. S3.4 Computing Spatial autocorrelation at lag-1 and spatial DFT Spatial autocorrelation at lag-1 (C) was computed using the standard definition !=

!!

! !,!,!,!!! !

!

!, !; !, ! ! !, ! − ! (! !, ! − !) ! !,!!!(!

!, ! − !)!

where N2 is the total number of entries in the matrix, W is the total number of near neighbor pairs, w[i,j;m,n] is 1 if (i,j) and (m,n) are near neighbor entries in the vegetation matrix, z[i,j] is the value of the vegetation at matrix element at (i,j) and ! is the spatial mean (Dakos et al. 2010). To obtain DFT at low frequencies (Carpenter & Brock 2010) we first computed the full DFT using MATLAB function fft2. We then computed the average of the magnitude of DFT at low frequencies. The low frequencies were defined as those up to 10% of the maximum frequency along both dimensions of the DFT (see Appendix S5 for sensitivity analysis).

S3.5 References Fortin and Dale, 2005, Spatial Analysis: A guide for ecologists. Cambridge University Press, Cambridge, UK. Bodine, Lenhart and Gross, 2014, Mathematics for the Life Sciences, Princeton University Press, Princeton, NJ, USA

Appendix S4: Results for other transects

Figure S4 Grass cover and rainfall along distance for each transect. Grass cover (open circles) and mean annual rainfall (dots) for distance moved along a transect for transects one, two, three, four, six, seven and nine. These graphs show that there is a gradient of both grass cover and rainfall along all transects. See Table S3 for quantification and statistics of relationship.

Figure S5 Grass cover as a function of rainfall for each transects one, two, three, four, six, seven and nine. These figures broadly support the claim of the main text that rainfall could be a driver of proportion of grass cover along transects. See Table S3 for quantification and statistics of relationships.

Table S3 Statistics for correlation coefficients (Kendall’s tau) for grass cover over, mean annual rainfall over distance, and grass cover versus mean annual rainfall over distance.

Grass cover over distance

Mean annual rainfall over distance 95% CI correlation based on bootstrap

correlation coefficient

95% CI correlation based on bootstrap

Grass cover versus mean annual rainfall over distance p-value

95% CI correlation based on bootstrap

-0.86

< 0.001

[-0.94, -0.76]

-0.52 -0.75

< 0.001 < 0.001

[-0.75, -0.22] [-0.92, -0.51]

[-0.05, 0.590] [0.63, 0.96] [-0.239, 0.749]

-0.27 -0.87 -0.34

0.045 < 0.001 0.0012

[-0.59, 0.03] [-0.96, -0.75] [-0.70, 0.11]

< 0.001 < 0.001

[1, 1] [1, 1]

-0.78 -0.88

< 0.001 < 0.001

[-1, -0.35] [-0.99, -0.66]

< 0.001

[1, 1]

-0.82

< 0.001

[-0.95, -0.61]

Transect Number

correlation coefficient

p-value

1

-0.86

< 0.001

[-0.94, -0.76]

0.98

< 0.001

[0.919, 1]

2 3

-0.55 -0.74

< 0.001 < 0.001

[-0.77, -0.26] [-0.92, -0.49]

0.92 0.99

< 0.001 < 0.001

[0.823, 0.988] [0.946, 1]

4 5 6

-0.85 -0.87 -0.84

< 0.001 < 0.001 < 0.001

[-0.95, -0.71] [-0.97, -0.74] [-0.97, -0.63]

0.26 0.83 0.31

0.061 < 0.001 0.025

7 8

-0.78 -0.88

< 0.001 < 0.001

[-1, -0.35] [-0.99, -0.66]

1 1

9

-0.82

< 0.001

[-0.95, -0.61]

1

p-value

correlation coefficient

Figure S6: Spatial indicators (over rainfall gradient) for all transects except 5 and 8 (results of transects 5 and 8 were shown in the main text.) Open blue circles correspond to real data and black stars correspond to null data. See Table S4 for quantification of trends for indicators for transects 5, 8, and 9. All other transects had very few points (less than six) and hence, we did not quantify their trends.

Figure S7: Spatial indicators (over spatial gradient) for all transects. Open blue circles correspond to real data and black stars correspond to null data. See Table S5 for quantification of trends for indicators. Transects 3 and 7 had very few points (less than six); hence, we did not quantify the trends of indicators for these

Appendix S5: Rainfall control transects and trends of indicators We chose control transects using the following two criteria: (a) Driver (rainfall) values must remain roughly the same along transects. We did this because all of the nine transects shown in figure 1 are located across strong rainfall gradients. Therefore, transects located across an area with relatively constant rainfall act as a control for this driver. (b) Driver values along these transects must be in a range at which the system will likely be in the vegetation state for which we compute spatial indicators (grasslands). Based on these two criteria and the results from figure 3 of the main text, we chose to locate our three transects in the low-rainfall range of 500-650 mm/year (figure S8). We found that the vegetation state is one of high grass cover throughout these transects (top row of figure S9). We then computed spatial indicators for each of these three transects using the same methods described in the main text.

Figure S8: Rainfall map of the Serengeti-Mara ecosystem with the locations of the rainfall ‘control’ transects (rectangles). Transects were placed in areas that receive low amounts of rainfall (500 – 650 mm/year) relative to the rest of the ecosystem.

!

Figure S9 Rainfall control transects show either no or irregular patterns of spatial indicators . The graphs show grass cover, spatial variance, spatial skewness, spatial correlation, and spatial low-frequency spectra for the transects moving from south to north. The top row shows that vegetation state (grassland) does not change for each transect, which is consistent with our results in the main text that demonstrated that low rainfall supports a stable grassland state (figure 3). The trends of these indicators were quantified using slope of the linear regression (Table S5).

Appendix S5: Sensitivity analyses We conducted sensitivity analyses to test the robustness of our results to the size at which we coarse-grained the data, the width of transects, and the low frequency that we selected for DFT. We employed smaller (a = 1 and 2) and larger (a = 10) coarse-graining sizes in addition to a = 5. In addition to the 7.5 km transect width analyzed above, we analyzed transects of widths 2.5 km and 5 km, while keeping their lengths and locations the same. To compute DFT at low frequencies, in addition to the one used above (10% of maximum frequency), we also used 20% and 30% of the total frequency components of the spectrum.

Figure S10 The figure shows that coarse-graining the data is a crucial step necessary to allow for the interpretation of spatial variance and spatial skewness. The figure also shows that the results are valid for a range of coarse-graining dimensions mentioned at the top of each column. Raw vegetation data are 30 m x 30 m in resolution. Therefore, coarse-grained vegetation data of: 2 x 2 corresponds to 60 m x 60 m resolution; 5 x 5 corresponds to 150 m x 150 m resolution; and 10 x 10 corresponds to 300 m x 300 m resolution. Blue open circles represent real data whereas black stars represent the null model. The grass cover for the vegetation data and null model data are identical by design (see Methods and Appendix S3), hence these data points overlap for all coarse-graining lengths (top row). However, this is not true for CSD-based spatial indicators. The raw data column shows that spatial variance and

spatial skewness are ineffective as indicators in the raw vegetation data because values of these two indicators are exactly the same as for the null model data (see Appendix S3 above). However, after coarse-graining, the spatial indicators for the vegetation data are different from those of the null data. Furthermore, coarse-graining also improves the distinction between the vegetation data and the null data for the spatial low-frequency spectra.

Figure S11 Spatial indicators are not sensitive to choice of transect width, as shown for transect 5. Blue open circles represent vegetation data and black stars represent null model data. We did not compute Spatial Low-Frequency Spectra for the transect with width 2.49 km since the size of the spectral matrix is small (25 x 25 pixels, of which the low frequency range corresponds to a very small sized matrix of 3 x 3.

Figure S12 Results are not sensitive to the choice of low frequency cut-offs (10%, 20% or 30%) used to compute DFT at low frequencies, shown for transect 5. Blue open circles represent real data and black stars represent null model data.