level resistance and resilience - Wiley Online Library

Utilizing gradient simulations for quantifying community-level resistance and resilience KARL A. LAMOTHE, DONALD A. JACKSON, AND KEITH M. SOMERS Department of Ecology and Evolutionary Biology, University of Toronto, Toronto, Ontario M5S 3B2 Canada Citation: Lamothe, K. A., D. A. Jackson, and K. M. Somers. 2017. Utilizing gradient simulations for quantifying community-level resistance and resilience. Ecosphere 8(9):e01953. 10.1002/ecs2.1953

Abstract. Resilience is a complex, multidimensional property of ecosystems that describes how ecosystems respond to disturbance and likely results from the interactions of species and their environments across temporal and spatial scales. Due to the complexity in how ecosystems function and respond to disturbance, measuring resilience is a challenge. Gradient analysis provides a familiar, yet somewhat neglected framework for understanding and characterizing resilience. With simulations parameterized on existing biomonitoring data, we used distance-based measures in ordination space to characterize community-level resilience, here defined as a function of resistance and recovery. Our simulations and analyses involved five steps: (1) We generated regional species pools by simulating species distributions across environmental gradients; (2) we sampled from these regional species pools to emulate temporal changes in reference (i.e., minimally disturbed) and impacted communities responding to disturbance; (3) we performed ordinations on observations from both impacted and reference communities to summarize multivariate data; (4) we calculated distance-based measures for individual community trajectories in the ordinations to quantify their relative resistance and resilience; and (5) we compared these distance-based metrics between reference and impacted communities. We conclude with an empirical example demonstrating the lack of resistance of the Harp Lake (Ontario, Canada) zooplankton community to invasion relative to the changes observed among minimally disturbed reference communities. Overall, distance measures on ordinations provide a simple and effective visual framework to quantify the relative resistance and resilience of communities to disturbance, and our simulation approach provides a novel technique to develop and evaluate quantitative metrics related to ecosystem or community-level processes. Key words: coenoplane; gradient analysis; ordination; recovery; resilience; resistance; simulation. Received 1 August 2017; accepted 10 August 2017. Corresponding Editor: Debra P. C. Peters. Copyright: © 2017 Lamothe et al. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. E-mail: [email protected]

INTRODUCTION

quantifying resilience (Angeler and Allen 2016). Arguably, the most common definition of resilience in the ecological literature is that of Holling (1973:17), who defined resilience as “a measure of the ability of [these] systems to absorb changes of state variables, driving variables, and parameters, and still persist.” This definition, now termed “ecological resilience,” describes resilience as a property of ecosystems and assumes that ecosystems can persist in multiple states (i.e., stability domains; Holling 1996). In contrast, “engineering resilience” describes the rate of

Resilience has become an increasingly important term and concept in ecology as evidenced from the increased frequency of citations in journal articles (Hodgson et al. 2015), discussion in books (Berkes and Folke 1998, Biggs et al. 2015), and the topic of special issues of academic journals (see Ecosystems, Carpenter et al. 2005; Journal of Applied Ecology, Angeler and Allen 2016). Because of its long-term use, there are many definitions of resilience and numerous techniques for ❖ www.esajournals.org

1

September 2017

❖ Volume 8(9) ❖ Article e01953

LAMOTHE ET AL.

nature (e.g., ordination axes, phase spaces that are not directly measurable during field sampling), but are useful for explaining community patterns (Austin 1985, ter Braak and Prentice 1988). A coenocline is a type of gradient that describes how species abundances change across a single environmental parameter, for example, elevation (Whittaker 1960, 1967). Similarly, a coenoplane describes how species abundances change across two orthogonal gradients (Fig. 1). Simulations of coenoclines and coenoplanes based on patterns of species in nature have been used for testing the reliability of many quantitative techniques (Gauch and Wentworth 1976, Gauch and Whittaker 1976, Hirst and Jackson 2007). Simulations of coenoplanes coupled with a space-for-time substitution sampling approach can be used to develop data sets where communities undergo known changes in composition. Space-for-time substitution uses successive sampling from different spatial locations to represent changes in species composition that would occur over time (e.g., species succession; Pickett 1989). By placing a grid on a coenoplane, we create a spatial landscape that can be sampled using a space-for-time approach (Fig. 1). At each intersecting point along the gridded coenoplane is a snapshot of a community where abundances of individual species are known and vary across the coenoplane (Fig. 1). Successive sampling at points of intersection along the gridded coenoplane results in a sequence of community composition data whereby abundances of species may vary over the order of observations (or time). Sampling repeatedly at a given point of intersection on the

return of an ecosystem to an equilibrium state following a disturbance (Pimm 1984, Holling 1996). Conceptually similar to the engineering definition and somewhat nested within the ecological definition, ecologists have also advocated for a resistance–resilience framework, where resistance describes the ability to persist during a disturbance event and resilience describes the ability to recover following a disturbance event (Nimmo et al. 2015). Whereas ecological resilience is often viewed at the ecosystem scale, the resistance–resilience framework can be applied at multiple biological levels, in turn informing us about the overall resilience of that system (Nimmo et al. 2015). A suite of novel quantitative approaches has been developed to characterize the resilience of systems to disturbance (Table 1). The variety of approaches developed has been fueled by differences in definition (Brand and Jax 2007, MyersSmith et al. 2012), the scale at which studies take place (Angeler and Allen 2016), the specific questions being addressed, the disturbances that are impacting the system (Carpenter et al. 2001), and the data available or required to answer the questions. One set of approaches that has been relatively underutilized for quantifying resilience is gradient analyses. Ecological gradient analyses were developed by plant community ecologists to understand patterns of species distributions and abundances along environmental gradients (Whittaker 1967). Gradient is a general term that can be used to describe spatial gradients, environmental gradients, or abstract gradients that may not occur in

Table 1. Examples of techniques/surrogates used for quantifying resilience. Method

Description

Functional redundancy

Multiple species in a community perform the same ecological function providing a form of insurance when a disturbance occurs Populations recover to historical population sizes or biomass levels quickly from severe disturbance events. Lack of resilience can be predicted using early warning signals Diversity in the response of species to disturbances maintains functional patterns and processes of ecosystems

Population models Response diversity Species richness Systems models

Sample references

Species-rich communities can better buffer environmental variability and are likely to contain species showing differing responses When facing disturbance events, resilient ecosystems mimic those absent of disturbance events

❖ www.esajournals.org

2

September 2017

Angeler and Allen (2016), Laliberte et al. (2010), Angeler et al. (2013) Dakos et al. (2008), Scheffer et al. (2009), Dai et al. (2012) Elmqvist et al. (2003), Mori et al. (2013), Baskett et al. (2014) Yachi and Loreau (1999), Downing and Leibold (2010) Bennett et al. (2005), Mumby and Anthony (2015), Marzloff et al. 2015


LAMOTHE ET AL.

Fig. 1. A coenoplane showing the abundance relationships for 50 species along two orthogonal gradients (e.g., environmental gradients). Each numbered point represents an individual observation, where species within a community vary in their abundance.

trajectories of change as a framework for the space-for-time substitution sampling on simulated coenoplanes (Fig. 2g), we can mimic biomonitoring data sets where communities show differing trajectories of change in community composition and abundance (Fig. 2), allowing comparisons of community-level resilience metrics. Here, we use coenoplane simulations and the hypothetical trajectories described by Matthews et al. (2013) to demonstrate a distance-based multivariate approach to quantify the relative resilience of communities to disturbance that embodies the properties of both resistance and recovery. The objective of our simulations is to represent actual community biomonitoring data to provide a practical approach for characterizing the relative resilience of communities to disturbances. We base our coenoplane simulations on long-term crustacean zooplankton community data from lakes in Ontario, Canada, where the sampling programs used a consistent sampling protocol, ensuring reliable comparisons of zooplankton communities over time (Yan et al. 2008, Palmer and Yan 2013; Appendix S1: Table S1). We conclude with an empirical example, characterizing the resistance and resilience of the Harp Lake zooplankton community to invasion relative to minimally disturbed reference zooplankton communities.

simulated coenoplane grid will result in identical species composition measurements over time representing a constant species composition and abundance. Alternatively, successive sampling of points that are more distant from each other will result in a sequence of observations that have less in common, both in terms of species abundance and occurrence (Fig. 1). Matthews et al. (2013) described six hypothetical trajectories of temporal change in communities that can be represented spatially along a coenoplane (Fig. 2). Community trajectory patterns were described as (1) saltatory or gradual, and (2) non-directional, directional, or directional followed by a return to a historical state. Small, gradual non-directional changes in communities can reflect trajectories of natural community variation exhibiting random increments of change in magnitude and direction (i.e., loose equilibrium; Matthews and Marsh-Matthews 2016; Fig. 2a), whereas saltatorial, abrupt responses, can reflect the responses of communities to “pulse” disturbances (Fig. 2b) such as extreme hydroclimatic events or wildfires. Gradual and saltatorial directional responses can depict permanent changes in community composition (Fig. 2b, e). In some cases, recovery to a historical state can occur (Fig. 2c, f), but such recovery typically requires gradual steps (Keller et al. 2002). Using these six ❖ www.esajournals.org

3

September 2017


LAMOTHE ET AL.

Fig. 2. Conceptual reconfigurations of gradual and saltatorial community trajectories over two environmental gradients (adapted from Matthews et al. 2013). (a) Trajectory of a reference community undergoing nondirectional gradual change in composition. (b–f) Impacted communities showing gradual (b, c) and saltatorial (d–f) responses. Saltatorial responses are depicted as having greater community turnover than gradual responses. Periods of relative constancy in community composition are depicted with filled circles. (g) Space-for-time sampling approach where each trajectory represents a single community.

MATERIALS AND METHODS

a region defined by two orthogonal gradients forming a coenoplane. We then sampled the coenoplane using a space-for-time substitution approach following the hypothetical trajectories

Our simulations involved five steps. We first simulated the distributions of all species found in ❖ www.esajournals.org

4

September 2017


LAMOTHE ET AL.

described by Matthews et al. (2013; Fig. 2). We characterize community trajectories as either reference communities exhibiting natural variability (Fig. 2a) or impacted communities showing varying responses to disturbance (Fig. 2b–f). The third step ordinated community data in a lowerdimensional space. Using the ordination results, we calculated distance-based metrics to characterize the relative resilience of individual communities. Finally, to follow a common approach in biomonitoring studies approach (Bowman and Somers 2005, Palmer et al. 2013), we compared these metrics between individual impacted communities and reference communities (See Appendix S1: Fig. S1 for flow chart of steps taken).

species richness levels across the grid. Centroids of species distributions did not always fall within the 100 9 100 grid, and this represents the situation in nature where a species’ optimum along an environmental gradient may fall beyond the range of conditions sampled. Each point of intersection within the grid represents a potential sampling location and represents a snapshot of species abundances with different species composition and potentially different richness. Differences in species composition across the simulated coenoplanes (i.e., gradient lengths) were assessed among coenoplanes using detrended correspondence analysis, and it was 5.97 SD (range: 5.12– 6.93), similar to the reported average of 4.10 (range: 1.36–11.98) from terrestrial and aquatic field studies (Hirst and Jackson 2007). A measure of approximately four SD units represents complete species turnover between the two ends of the gradient (Gauch 1982).

Terminology We adopt the following terminology to describe our simulations and analyses. Here, “sampling” describes the act of locating an intersecting point along the simulated coenoplane and extracting (with replacement) the underlying community composition snapshot, or “sample.” This sample is an “observation,” or a snapshot of community composition consisting of the abundances of all simulated species present at that location. Here, a “community” refers to 30 observations collected along a single trajectory based on Matthews et al. (2013). We refer to “reference communities” as minimally impacted communities demonstrating gradual, non-directional trajectories (Fig. 2a) and “impacted communities” as those showing some other response (Fig. 2b–f). Finally, we used the term “constancy” to describe the relative consistency of community composition over time.

Step 2 sampling communities from the simulated coenoplanes We sampled the coenoplanes to resemble the hypothetical trajectories of communities described by Matthews et al. (2013; Fig. 2). Each scenario represents the trajectory of an individual community, and consists of 30 observations of species abundance data, consistent with freshwater zooplankton biomonitoring data from Ontario, Canada (i.e., 30 yr; Appendix S1: Table S1). For all six trajectories, communities were sampled to display 10 initial time steps of relative constancy in composition (i.e., minimal changes in abundance) prior to showing a response. One hundred replicates of each of the community trajectories were sampled from each of the 10 coenoplanes, thereby providing a total of 1000 possible replications of each scenario (Fig. 2a–f). Trajectories of recovery were set to track back to historical centroids while incorporating noise into the sampled coenoplane coordinates.

Step 1 simulating coenoplanes We simulated hypothetical species distributions along two orthogonal environmental gradients forming a coenoplane using the coenoflex package (Roberts 2016) in the R statistical software (R Core Team 2016). We simulated 10 environmental coenoplanes, containing 50 species each and constructed to be 100 9 100 units in length-by-width. The centroids of individual species distributions (i.e., species optima) were randomly dispersed across each gradient with approximately constant density throughout the sample space. Species abundances were simulated to have a Gaussian response along each environmental gradient. This random distribution of species forms a mosaic of ❖ www.esajournals.org

Step 3 multivariate analyses We used principal component analyses (PCAs) on Hellinger-transformed species abundance matrices to ordinate the species composition data in fewer dimensions. Hellinger transformations are a recommended distance measure for the ordination of species and were used to reduce the impact of rare species on 5

September 2017


LAMOTHE ET AL.

Step 4 distance-based metrics

ordination results (Legendre and Legendre 1998, Legendre and Gallagher 2001). For our ordinations, we consider situations that range from having many more reference communities (Fig. 2a) to impacted communities (Fig. 2b–f) to the reverse situation. For example, we performed 100 replicate ordinations where the analyses included 50 communities showing gradual, non-directional change (i.e., reference) and 10 communities showing any of the other five hypothetical trajectories (i.e., impacted). In this situation, the community data matrix contained 60 communities each consisting of 30 observations totaling 1800 observations. In the reverse situation, we performed a total of 100 replicate ordinations that each included 10 reference and 50 impacted communities. Finally, we performed 100 replicate ordinations of each of the ratios where the analyses contained 20 reference communities: 4 impacted communities (20:4), 4:20, 30:30, and 10:10. Biomonitoring programs tend to vary in the number of reference and impacted community data sets available, and therefore, understanding how the distance-based framework performs with differing numbers of communities is important. However, for the sake of brevity and clarity, we only present results for the situation where 50 communities showing gradual, non-directional change and 10 communities showing any of the other five hypothetical trajectories are considered (ideal scenario); results for differing sample size scenarios can be found in Appendix S1. Communities used in the ordination analyses were randomly selected from the communities described in Step 2. For example, in the analysis of 50 reference and 10 impacted communities, 50 reference communities were sampled without replacement from the 1000 simulated communities undergoing gradual, non-directional change (Fig. 2a) and 10 impacted communities were selected without replacement from any of the 5000 simulated communities undergoing either gradual directional change, gradual directional change with recovery, saltatory nondirectional change, saltatory directional change, and saltatory directional change with recovery (Fig. 2b–f). Two principal component axes were extracted from each ordination because the simulations were based on two orthogonal gradients. ❖ www.esajournals.org

Distances between observations in an ordination represent differences between observations; greater distances indicate greater differences in species composition. Communities experiencing little change over time will travel relatively small distances in multivariate space over time compared to communities exhibiting major departures from a historical state. Resistance and recovery were assessed by measuring the Euclidean distance of each observation relative to a baseline centroid (db; Anderson and Thompson 2004, Milner et al. 2016). We calculated distances to a baseline centroid of 10 observations for each individual community as each trajectory initially had 10-time steps of relative constancy in community composition prior to subsequent responses (Fig. 2). Low db values over time indicate a community maintaining its composition around some equilibrial state or stability domain. Additionally, distances were calculated from each observation to the centroid defined by all prior observations within a trajectory (i.e., cumulative centroid: dc; Anderson and Thompson 2004, Milner et al. 2016). Sequential distances to the cumulative centroid describe the magnitude of change in community composition in sequential time steps, with the potential for detecting pulse disturbances (Anderson and Thompson 2004, Milner et al. 2016). Short distances to a cumulative centroid over time indicate small changes in variation in community composition over time, whereas a large distance within a series of short distances can indicate the effects of a pulse disturbance on the system. To assess the ability of the ordinations to accurately capture changes in multivariate space, Euclidean distances to cumulative (dmv.c) and baseline (dmv.b) centroids were calculated from the multivariate Hellinger-transformed abundance matrices. We then calculated Pearson correlations between these full multivariate distances and the distance-based metrics calculated from the PCA site scores (db and dc) to assess any differences arising due to working in the reduced dimensionality of the PCA axes.

Step 5 comparing distance-based metrics between reference and impacted communities We compared results from distance-based metrics (db and dc) for individual impacted communities with the results for communities 6

September 2017


LAMOTHE ET AL.

undergoing gradual, non-directional change (i.e., reference communities). Ninety-five percent confidence intervals were calculated for the set of reference communities to define the boundary of distances traveled in multivariate space over time by reference communities. Distance values falling outside this reference distribution indicate changes in community composition that exceed expectations of natural variability and therefore a lack of resistance to change. To further evaluate deviations in db and dc from historical conditions among individual communities, we applied the sequential regime shift detector (SRSD; Rodionov 2004). This methodology was developed to detect nontrivial shifts in time-series data using sequential t-tests, where each additional value in a time series necessitates a new test of whether there is a significant difference between the new value and the previous l time steps. A regime shift index (RSI) is calculated as the cumulative sum of normalized differences from the hypothetical mean for a new regime (Rodionov and Overland 2005). Positive RSI values indicate shifts away from a centroid, and negative RSI values indicate shifts toward a centroid. Significant shifts to alternative states are largely determined by choices in the minimum length of l time steps to compare and by the chosen significance level (P). For the purposes of this study, we knew the length of community trajectories a priori, and chose l = 5 and P = 0.01.

(Yan and Pawson 1997, Yan et al. 2001, 2008). We combine data from Harp Lake with zooplankton monitoring data from two nearby minimally disturbed reference lakes (Red Chalk and Blue Chalk lakes) and seven minimally disturbed reference lakes from the Experimental Lakes Area (ELA now known as the IISD-ELA) in northwestern Ontario (lakes 224, 239, 240, 373, 375, 382, and 442). Reference lakes were chosen to reflect similar conditions to Harp Lake, being relatively free of impacts from anthropogenic stressors and, particularly, the presence or invasion of Bythotrephes. Methods for sampling zooplankton differed between lakes, but were generally consistent within lakes over time. Several zooplankton species or taxonomic groups were combined due to differences in identification and changes in nomenclature over time (see Palmer et al. 2013 for details). We calculated a PCA on the covariance matrix of annual-average crustacean zooplankton species abundances for all lakes and years. Prior to the ordination, Hellinger transformations were performed to reduce the impact of rare species on the ordination results (Legendre and Gallagher 2001). The number of nontrivial components to retain was chosen using a permutation approach; we permuted each column of the Hellinger-transformed zooplankton annual-average species abundance matrix and conducted a PCA 999 times. We retained a component if the proportion of variance explained in the empirical data exceeded 95% of the permuted PCAs for that component (Peres-Neto et al. 2003, 2005). We calculated and compared db and dc over time between the Harp Lake zooplankton community and the group of nine reference sites. We used 95% confidence intervals from the set of reference communities to define the boundaries of reference-site distances traveled in ordination space. Furthermore, we applied the SRSD test (Rodionov 2004) to test for nontrivial shifts in db and dc over time using l = 5 and P = 0.01. All statistical analyses were conducted in base R (version 3.2.3; R Core Team 2016) using functions from the packages “coenoflex” (Roberts 2016), “Hmisc” (Harrell 2016), “vegan” (Oksanen et al. 2016), and “ggplot2” (Wickham 2009), or Excel (SRSD add-in; Rodionov 2004). Sample code and metadata can be found in Data S1.

Characterizing the resistance and resilience of a zooplankton community to an invasion We use long-term (1980–2012) crustacean zooplankton species abundance data from Ontario, Canada, to demonstrate the utility of the distance-based approach for characterizing the resistance and resilience of communities to disturbance. Harp Lake is a relatively small lake (71.4 ha) located in southcentral Ontario (Yan et al. 2001). Crustacean zooplankton communities have been sampled from Harp Lake since the early 1980s (Yan and Pawson 1997). During this period, a nonnative zooplanktivore, Bythotrephes, appeared and has led to dramatic changes in the Harp Lake zooplankton community including a decline in species richness, reduction in body size, and decline in total zooplankton abundance


7

September 2017


LAMOTHE ET AL.

RESULTS

distance-based analyses on multivariate data (db and dmv.b range: 0.83–0.88; dc and dmv.c range: 0.73–0.85).

We performed 600 ordinations with differing numbers of reference and impacted communities (i.e., community sample sizes). Here, we present more detailed results of a single ordination containing 50 reference and 10 impacted communities and summarize the trends across community sample sizes. Detailed results on the various sample sizes can be found in Appendix S1. We report distance measures performed on ordination results because they allow us to graphically display the relationships, and because there are strong correlations with the results from parallel

Example of a single analysis containing 50 reference and 10 impacted communities Site scores for communities showing gradual, non-directional change were typically clustered, whereas impacted communities were more widely scattered across the two components (Fig. 3b). This outcome was expected and reflects the patterns of sampling (Fig. 3a): Sampling of impacted communities covered more of the coenoplane than reference communities showing

Fig. 3. (a) Sampling coordinates from the coenoplane for the 60 communities within the single ordination. Point shapes refer to inset hypothetical trajectories. Lines connecting impacted community sampling coordinates show their individual temporal trajectories. One community undergoing saltatory directional change with recovery highlighted in red (Fig. 2f). (b) principal component analysis ordination plot showing site scores for sampled communities. (c) An enlarged view of one impacted community undergoing saltatorial directional change with recovery and one reference community (Fig. 2a). Observations are connected in order of sampling. (d) Distances to baseline centroids (white shapes) over time for both the single impacted and reference community. (e) Distances to cumulative centroids (white shapes) over time for both the impacted and reference community. (f) Distance to baseline centroid over time for both communities. Dashed line indicates significant regime shift (P < 0.01). (g) Distance to cumulative centroid over time for both communities.


8

September 2017


LAMOTHE ET AL.

gradual, non-directional change (Fig. 3a). Looking more closely at individual communities, Fig. 3c displays the trajectory of a single reference community and a single impacted community undergoing a saltatorial directional response with recovery. The reference community shows variation around a mean location in multivariate space, whereas the impacted community travels around the ordination. Fig. 3d, e displays information used in db and dc calculations, respectively, for both the individual impacted and reference communities. Compared to the reference community, db of the impacted community deviated from the reference community trajectory at time-step 11 and showed a significant shift away from the initial community state (RSI = 1.414, P < 0.001; Fig. 3f). The impacted community then showed a period of relative constancy in db measurements and subsequently returned to the reference community trajectory around time-step 24 with a significant shift back toward the reference (RSI = 1.152, P < 0.001). Alternatively, the impacted community showed a greater dc over time than the reference community through the entire time series (Fig. 3g). Reference communities exhibiting a gradual, non-directional response consistently showed low db and dc over time (Fig. 4). Low db and dc values were expected as sampled reference communities lacked major change throughout the simulations. For db, deviations of impacted communities from the reference distribution occurred immediately following the initial period of relative constancy in composition (time steps 1–10), indicating a lack of resistance. Measures of db for communities showing directional trajectories (Fig. 2b, e) remained outside the range of reference communities following the initial period of relative constancy in composition (Fig. 4a), whereas communities showing resilience (Fig. 2c, f) returned to their initial state (Fig. 4a). Comparisons using dc from the PCA ordinations indicated substantial change occurring among the differing impacted communities (Fig. 4b). Impacted communities remained outside the 95% reference confidence interval for most of the time series (Fig. 4b). Distances to cumulative centroids of impacted communities were typically large during the initial 10-time steps of relative constancy in community composition, but then differed based on the community ❖ www.esajournals.org

Fig. 4. (a) Examples of the Euclidean distances to a baseline centroid from the principal component analysis over time for one ordination containing 50 reference and 10 impacted communities. Plotted are representative examples of each inset scenario. Symbols match with inset scenarios. Mean distances over time are plotted for the reference communities with 95% confidence intervals (gray ribbon). (b) Examples of the Euclidean distances to a cumulative centroid over time for the same representative examples.

trajectory. Following the initial 10 yr of constancy in community composition, communities showing gradual directional change (Fig. 2b) tracked toward the centroid and subsequently deviated away (Fig. 4b). Similarly, when communities showed patterns of recovery (Fig. 2c, f), distances to cumulative centroids were reduced and followed by subsequent deviations from the centroid (Fig. 4b). Patterns for communities undergoing saltatorial non-directional change (Fig. 2d) peaked based on the magnitude of deviation from the cumulative centroid (Fig. 4b). 9

September 2017


LAMOTHE ET AL.

Finally, following the first 10 observations, communities undergoing saltatorial directional change (Fig. 2e) tracked toward and then subsequently deviated from the cumulative centroid (Fig. 4b). Further, measures of dc were consistently lower in the second half of the time series compared to the first due to the longer period of constancy following the saltatorial change (18 times steps vs. 10).

(Fig. 4b). Similarly, the pulse disturbance of the invasion led to a substantial drop in dc scores. This drop is indicative of the trajectory in ordination space; there were an initial 13 yr of constancy in community composition followed by the invasion where the community deviated from the historical region in ordination space (i.e., upper-left quadrant) for approximately 3 yr and then remained in the new region of ordination space (i.e., bottom-right quadrant) for approximately 14 yr (Fig. 5a).

Trends across community sample sizes Despite the differences in sample size across analyses, variation in db and dc among reference communities remained low (Appendix S1: Figs. S2, S3). In contrast, variability among impacted communities increased as the number of impacted communities within the analysis decreased. However, despite the increase in variability across ordinations containing differing ratios of reference to impacted communities, trajectory patterns remained consistent for db and dc measures, allowing for the consistent characterization of resistance and resilience. In all cases, impacted communities showed a period of relative constancy in db and dc measures followed by consistent patterns based on the individual trajectories and distance metric (Appendix S1: Figs. S2, S3).

DISCUSSION Given the inherent errors in sampling environmental data and lack of knowledge about how ecological communities respond to disturbance, comparing statistical methodologies with empirical field data may display significant differences in results, but the results may remain unclear about which methods best reflect true underlying relationships. Simulations provide the opportunity to evaluate quantitative approaches and aid in the development and assessment of novel metrics. Here, we used coenoplane simulations coupled with a space-for-time substitution approach for sampling to build community-level data sets and subsequently used a distance-based framework to characterize the relative resilience of potentially impacted communities to disturbance over time. Further, we applied this approach to zooplankton monitoring data from Ontario, Canada, and demonstrated a lack of resistance to change among the Harp Lake zooplankton community resulting from the invasion of Bythotrephes (Yan and Pawson 1997, Yan et al. 2001, 2008). Trajectories of distances from an historical state (db) over time provide a relative metric of both resistance and recovery, whereas distances to a community-level centroid (dc) provide an understanding of the magnitude of change over time (Milner et al. 2016). Generally, communities showing the shortest distances traveled in ordination space, and therefore having minimal changes in species composition, were considered to represent the greatest level of resistance to change. In contrast, substantial deviations from historical ranges of natural variation, or alternatively the variation of reference communities, indicated a compositional change in the community, and therefore a lack of resistance and

Zooplankton community example The zooplankton community from Harp Lake followed a similar trajectory to simulated communities undergoing saltatory, directional change over time (Fig. 2e), beginning in the upper-left quadrant of ordination space in the early 1980s and finishing in the bottom-right quadrant of ordination space during the late 2000s (Fig. 5a). We observed a significant shift in db in 1993 (RSI = 2.70, P < 0.001) where the community deviated from the reference distribution and subsequently maintained that distance away from the historical centroid in ordination space (Fig. 5b) indicating a lack of resistance to change. This is also the year in which Bythotrephes was first observed in Harp Lake. Furthermore, measures of dc remained outside the reference distribution for most of the sampling period (Fig. 5c), with a decreasing trend over time. The decreasing trend in dc over time is consistent with the patterns observed among simulated communities undergoing saltatory, directional change ❖ www.esajournals.org

10

September 2017


LAMOTHE ET AL.

Fig. 5. (a) Principal component analysis ordination plot of the 10 zooplankton communities. Highlighted in red is the trajectory for the Harp Lake zooplankton community. Distances to (b) baseline (db) and (c) cumulative centroids (dc) over time for the Harp Lake zooplankton community (red) and reference communities (white). Mean distances over time are plotted for the reference communities with 95% confidence intervals (gray ribbon). A significant shift in db occurred in 1993 for the Harp Lake zooplankton community (P < 0.01).

resilience (Seidl et al. 2016). In the case of the Harp Lake zooplankton community, the trajectories of db over time prior to the invasion of Bythotrephes were consistent with reference communities. However, after the invasion of Bythotrephes the composition of the community changed, evident from the sudden and significant shift away from the ❖ www.esajournals.org

reference distribution in 1993, indicating a lack of resistance to change. Furthermore, we have not seen a recovery back to the historical composition post-invasion, indicating a lack of resilience to invasion. We were able to characterize the resistance and resilience among all the simulated communities 11

September 2017


LAMOTHE ET AL.

understanding of the resilience of these communities to disturbance and make inferences about how these communities may change in future. And although this approach relies on the continuation of long-term monitoring programs, which can be difficult to fund and maintain, we are optimistic that the benefits of these programs for understanding the resilience of ecosystems to change will justify their perpetuation.

despite varying the ratio of reference to impacted sites included in the analyses. This is an important finding, as empirical studies are often limited by the number of reference sites available. However, the empirical example comparing Harp Lake to the nine reference systems demonstrates that we may have underestimated the variability demonstrated by reference systems, or alternatively, that reference systems are changing and demonstrate patterns that are inconsistent with our expectations (i.e., gradual, non-directional). Regardless, the changes observed in db and dc over time for the Harp Lake zooplankton community far exceeded that of the reference communities, providing a clear visual of the impact of invasion. The use of distance-based frameworks with ordinations is certainly not a new approach in ecology (Holmes et al. 1979, Holmes and Recher 1986, Evans 1988, Hughes 1990, Boulton et al. 1992, Kilgour et al. 1998, Anderson and Thompson 2004, Sponseller et al. 2010). In fact, distancebased frameworks have been used to develop a wealth of ecological metrics that have been described as important for resilience (TimpanePadgham et al. 2017) including measures of functional diversity (Walker et al. 1999, Laliberte and Legendre 2010) and beta diversity (Anderson et al. 2006). And although we applied this approach toward understanding the resilience of biotic communities to change, these techniques can be similarly applied to replicated environmental monitoring data such as water-quality parameters or nutrient measures. To conclude, our approach of simulating coenoplanes and subsequently sampling using a space-for-time substitution approach provides a novel means to generate community-level data that mimic that of biomonitoring programs like the freshwater zooplankton monitoring program in Ontario (Yan et al. 2008) or the stream fish monitoring in northwest Arkansas and Oklahoma (Matthews et al. 2013, Matthews and Marsh-Matthews 2016). We chose to simulate reference and impacted communities because the reference condition approach represents a common approach among biomonitoring programs (Bowman and Somers 2005, Palmer et al. 2013). By comparing distance measures of impacted community trajectories over time to those of reference communities, we can gain a better ❖ www.esajournals.org

ACKNOWLEDGMENTS We thank the members of the Donald A. Jackson lab, Brian J. Shuter, Marie-Josee Fortin, W. Bill Keller, and two anonymous reviewers for comments on earlier drafts of this manuscript. We also thank W. Bill Keller, Michael J. Paterson, and James A. Rusak for providing zooplankton community data that we used to parameterize our simulations and provide an empirical example. Funding for Karl A. Lamothe and Donald A. Jackson was provided by the Natural Sciences and Engineering Research Council (NSERC), Canadian Network for Aquatic Ecosystem Services (CNAES), the University of Toronto Connaught Scholarship for International Students awarded to Karl A. Lamothe, and an NSERC Discovery grant awarded to Donald A. Jackson.

LITERATURE CITED Anderson, M. J., K. E. Ellingson, and B. H. McArdle. 2006. Multivariate dispersion as a measure of beta diversity. Ecology Letters 9:683–693. Anderson, M. J., and A. A. Thompson. 2004. Multivariate control charts for ecological and environmental monitoring. Ecological Applications 14:1921–1935. Angeler, D. G., and C. R. Allen. 2016. Quantifying resilience. Journal of Applied Ecology 53:617–624. Angeler, D. G., C. R. Allen, and R. K. Johnson. 2013. Measuring the relative resilience of subarctic lakes to global change: redundancies of functions within and across temporal scales. Journal of Applied Ecology 50:572–584. Austin, M. P. 1985. Continuum concept, ordination methods, and niche theory. Annual Review of Ecology and Systematics 16:39–61. Baskett, M. L., N. S. Fabina, and K. Gross. 2014. Response diversity can increase ecological resilience to disturbance in coral reefs. American Naturalist 184:E16–E31. Bennett, E. M., G. S. Cumming, and G. D. Peterson. 2005. A systems model approach to determining resilience surrogates for case studies. Ecosystems 8:945–957.

12

September 2017


LAMOTHE ET AL. Harrell, F. E. 2016. Hmisc: Harrell miscellaneous. R package version 3.17-4. https://CRAN.R-project. org/package=Hmisc Hirst, C. N., and D. A. Jackson. 2007. Reconstructing community relationships: the impact of sampling error, ordination approach, and gradient length. Diversity and Distributions 13:361–371. Hodgson, D., J. L. McDonald, and D. J. Hosken. 2015. What do you mean, ‘resilient’? Trends in Ecology and Evolution 30:503–506. Holling, C. S. 1973. Resilience and stability of ecological systems. Annual Review of Ecology and Systematics 4:1–23. Holling, C. S. 1996. Engineering resilience versus ecological resilience. Pages 31–44 in P. Schulze, editor. Engineering with ecological constraints. National Academies Press, Washington, D.C., USA. Holmes, R. T., R. E. Bonney, and S. W. Pacala. 1979. Guild structure of the Hubbard Brook bird community: a multivariate approach. Ecology 60:512– 520. Holmes, R. T., and H. F. Recher. 1986. Determinants of guild structure in forest bird communities: an intercontinental comparison. Condor 88:427–439. Hughes, J. M. R. 1990. Lotic vegetation dynamics following disturbance along the Swan and Apsley Rivers, Tasmania, Australia. Journal of Biogeography 17:291–306. Keller, W., N. D. Yan, K. M. Somers, and J. H. Heneberry. 2002. Crustacean zooplankton communities in lakes recovering from acidification. Canadian Journal of Fisheries and Aquatic Sciences 59:726–735. Kilgour, B. W., K. M. Somers, and D. E. Matthews. 1998. Using the normal range as a criterion for ecological significance in environmental monitoring and assessment. Ecoscience 5:542–550. Laliberte, E., and P. Legendre. 2010. A distance based framework for measuring functional diversity from multiple traits. Ecology 91:299–305. Laliberte, E., et al. 2010. Land-use intensification reduces functional redundancy and response diversity in plant communities. Ecology Letters 13:76–86. Legendre, P., and E. D. Gallagher. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129:271–280. Legendre, P., and L. Legendre. 1998. Numerical ecology. Second English edition. Elsevier, Amsterdam, The Netherlands. Marzloff, M. P., L. R. Little, and C. R. Johnson. 2015. Building resilience against climate-driven shifts in a temperate reef system: staying away from context-dependent ecological thresholds. Ecosystems 19:1–15.

Berkes, F., and C. Folke. 1998. Linking social and ecological systems: management practice and social mechanisms for building resilience. Cambridge University Press, Cambridge, UK. € ter, and M. L. Schoon. 2015. PrinciBiggs, R., M. Schlu ples for building resilience: sustaining ecosystem services in social-ecological systems. Cambridge University Press, Cambridge, UK. Boulton, A. J., C. G. Peterson, N. B. Grimm, and S. G. Fisher. 1992. Stability of an aquatic macroinvertebrate community in a multiyear hydrologic disturbance regime. Ecology 73:2192–2207. Bowman, M. F., and K. M. Somers. 2005. Considerations when using the reference condition approach for bioassessment of freshwater ecosystems. Water Quality Research Journal of Canada 40:347–360. Brand, F. S., and K. Jax. 2007. Focusing the meaning(s) of resilience: resilience as a descriptive concept and a boundary object. Ecology and Society 12:23. Carpenter, S. R., B. Walker, J. M. Anderies, and N. Abel. 2001. From metaphor to measurement: Resilience of what to what? Ecosystems 4:765–781. Carpenter, S. R., F. Westley, and M. G. Turner. 2005. Surrogates for resilience of social-ecological systems. Ecosystems 8:941–944. Dai, L., D. Vorselen, K. S. Korolev, and J. Gore. 2012. Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336:1175–1177. Dakos, V., M. Scheffer, E. H. van Nes, V. Brovkin, V. Petoukhov, and H. Held. 2008. Slowing down as an early warning signal for abrupt climate change. Proceedings of the National Academy of Sciences of the United States of America 105:14308–14312. Downing, A. L., and M. A. Leibold. 2010. Species richness facilitates ecosystem resilience in aquatic food webs. Freshwater Biology 55:2123–2137. €m, G. Peterson, Elmqvist, T., C. Folke, M. Nystro J. Bengtsson, B. Walker, and J. Norberg. 2003. Response diversity, ecosystem change, and resilience. Frontiers in Ecology and the Environment 1:488–494. Evans, E. W. 1988. Community dynamics of prairie grasshoppers subjected to periodic fire: Predictable trajectories or random walks in time? Oikos 52:283–292. Gauch, H. G. 1982. Multivariate analysis in community ecology. Cambridge University Press, Cambridge, UK. Gauch, H. G., and T. R. Wentworth. 1976. Canonical correlation analysis as an ordination technique. Vegetatio 33:17–22. Gauch, H. G., and R. H. Whittaker. 1976. Simulation of community patterns. Vegetatio 33:13–16.


13

September 2017


LAMOTHE ET AL. approaches and alternatives. Springer-Verlag, New York, New York, USA. Pimm, S. L. 1984. The complexity and stability of ecosystems. Nature 307:321–326. R Core Team. 2016. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Roberts, D. W. 2016. R: Coenoflex: gradient-based coenospace vegetation simulator. R package version 2.1-0. https://cran.r-project.org/web/packages/ coenoflex/coenoflex.pdf Rodionov, S. N. 2004. A sequential algorithm for testing climate regime shifts. Geophysical Research Letters 31:L09204. Rodionov, S., and J. E. Overland. 2005. Application of a sequential regime shift detection method to the Bering Sea ecosystem. ICES Journal of Marine Science 62:328–332. Scheffer, M., J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G. Sugihara. 2009. Early-warning signals for critical transitions. Nature 461:53–59. Seidl, R., T. A. Spies, D. L. Peterson, S. L. Stephens, and J. A. Hicke. 2016. Searching for resilience: addressing the impacts of changing disturbance regimes on forest ecosystem services. Journal of Applied Ecology 53:120–129. Sponseller, R. A., N. B. Grimm, A. J. Boulton, and J. L. Sabo. 2010. Responses of macroinvertebrate communities to long-term variability in a Sonoran Desert system. Global Change Biology 16:2891– 2900. ter Braak, C. J. F., and C. Prentice. 1988. A theory of gradient analysis. Advances in Ecological Research 34:271–317. Timpane-Padgham, B. L., T. Beechie, and T. Klinger. 2017. A systematic review of ecological attributes that confer resilience to climate change in environmental restoration. PLoS ONE 12:e0173812. Walker, B., A. Kinzig, and J. Langridge. 1999. Plant attribute diversity and ecosystem function: the nature and significance of dominant and minor species. Ecosystems 2:95–113. Whittaker, R. H. 1960. Vegetation of the Siskiyou mountains, Oregon and California. Ecological Monographs 30:279–338. Whittaker, R. H. 1967. Gradient analysis of vegetation. Biological Reviews 42:207–264. Wickham, H. 2009. ggplot2: elegant graphics for data analysis. Springer-Verlag, New York, New York, USA. Yachi, S., and M. Loreau. 1999. Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis. Proceedings of the

Matthews, W. J., and E. Marsh-Matthews. 2016. Dynamics of an upland stream fish community over 40 years: trajectories and support for the loose equilibrium concept. Ecology 97:706–719. Matthews, W. J., E. Marsh-Matthews, R. C. Cashner, and F. Gelwick. 2013. Disturbance and trajectory of change in a stream fish community over four decades. Oecologia 173:955–969. Milner, A. M., A. Woodward, J. E. Freilich, R. W. Black, and V. H. Resh. 2016. Detecting significant change in stream benthic macroinvertebrate communities in wilderness areas. Ecological Indicators 60:524–537. Mori, A. S., T. Furukawa, and T. Sasaki. 2013. Response diversity determines the resilience of ecosystems to environmental change. Biological Reviews 88:349–364. Mumby, P. J., and K. R. N. Anthony. 2015. Resilience metrics to inform ecosystem management under global climate change with application to coral reefs. Methods in Ecology and Evolution 6:1088– 1096. Myers-Smith, I. H., S. A. Trefry, and V. J. Swarbrick. 2012. Resilience: easy to use but hard to define. Ideas in Ecology and Evolution 5:44–53. Nimmo, D. G., R. Mac Nally, S. C. Cunningham, A. Haslem, and A. F. Bennett. 2015. Vive la resistance: reviving resistance for 21st century conservation. Trends in Ecology and Evolution 30: 516–523. Oksanen, J., et al. 2016. vegan: community ecology package. R package version 2.4-1. https://CRAN. R-project.org/package=vegan Palmer, M. E., W. B. Keller, and N. D. Yan. 2013. Gauging recovery of zooplankton from historical acid and metal contamination: the influence of temporal changes in restoration targets. Journal of Applied Ecology 50:107–118. Palmer, M. E., and N. D. Yan. 2013. Decadal scale regional changes in Canadian freshwater zooplankton: the likely consequence of complex interactions among multiple anthropogenic stressors. Freshwater Biology 58:1366–1378. Peres-Neto, P. R., D. A. Jackson, and K. M. Somers. 2003. Giving meaningful interpretation to ordination axes: assessing loading significance in principal component analysis. Ecology 84:2347–2363. Peres-Neto, P. R., D. A. Jackson, and K. M. Somers. 2005. How many principal components? Stopping rules for determining the number of non-trivial axes revisited. Computational Statistics & Data Analysis 49:974–997. Pickett, S. T. A. 1989. Space for time substitution as an alternative to long-term studies. Pages 110–135 in G. E. Likens, editor. Long-term studies in ecology:


14

September 2017


LAMOTHE ET AL. Canada, following invasion by Bythotrephes cederstrœmi. Freshwater Biology 37:409–425. Yan, N. D., K. M. Somers, R. E. Girard, A. M. Paterson, W. B. Keller, C. W. Ramcharan, J. A. Rusak, R. Ingram, G. E. Morgan, and J. M. Gunn. 2008. Longterm trends in zooplankton of Dorset, Ontario, lakes: the probable interact effects of changes in pH, total phosphorus, dissolved organic carbon, and predators. Canadian Journal of Fisheries and Aquatic Sciences 65:862–887.

National Academy of Sciences of the United States of America 96:1463–1468. Yan, N. D., A. Blukacz, W. G. Sprules, P. K. Kindy, D. Hackett, R. E. Girard, and B. J. Clark. 2001. Changes in zooplankton and the phenology of the spiny water flea, Bythotrephes, following its invasion of Harp Lake, Ontario, Canada. Canadian Journal of Fisheries and Aquatic Sciences 58:2341–2350. Yan, N. D., and T. Pawson. 1997. Changes in the crustacean zooplankton community of Harp Lake,

SUPPORTING INFORMATION Additional Supporting Information may be found online at: http://onlinelibrary.wiley.com/doi/10.1002/ecs2. 1953/full


15

September 2017
