A new variance ratio metric to detect the timescale of compensatory dynamics

Understanding the mechanisms governing ecological stability – why a property such as primary productivity is stable in some communities and variable in others – has long been a focus of ecology. Compensatory dynamics, in which anti-synchronous fluctuations between populations buffer against fluctuations at the community level, is a key theoretical mechanism of stability. Classically, compensatory dynamics have been quantified using a “variance ratio” approach that compares the ratio between community variance and aggregate population variance, such that a lower ratio indicates compensation and a higher ratio indicates synchrony among species fluctuations. However, population dynamics may be influenced by different drivers that operate on different timescales, and evidence from aquatic systems indicates that communities can be compensatory on some timescales and synchronous on others. The variance ratio and related metrics cannot reflect this timescale-specificity, yet have remained popular, especially in terrestrial systems. Here, we develop a timescale-specific variance ratio approach that formally decomposes the classical variance ratio according to the timescales of distinct contributions. The approach is implemented in a new R package, called tsvr, that accompanies this paper. We apply our approach to a long-term, multi-site grassland community dataset. Our approach demonstrates that the degree of compensation versus synchrony in community dynamics can vary by timescale. Across sites, population variability was typically greater over longer compared to shorter timescales. At some sites, minimal timescale-specificity in compensatory dynamics translated this pattern of population variability into a similar pattern of greater community variability on longer compared to shorter timescales. But at other sites, differentially stronger compensatory dynamics at longer compared to shorter timescales produced lower-than-expected community variability on longer timescales. Within every site there were plots that exhibited shifts in the strength of compensation between timescales. Our results highlight that compensatory versus synchronous dynamics are intrinsically timescale-dependent concepts, and our timescale-specific variance ratio provides a metric to quantify timescale-specificity and relate it back to the classic variance ratio.


52
The stability of ecosystem functions is central to the reliable provisioning of ecosystem services 53 (Oliver et al., 2015), and understanding mechanisms underlying ecological stability is a 54 fundamental question in ecology (MacArthur, 1955). A key insight into ecosystem dynamics is that 55 stable aggregate ecosystem functions such as total productivity can be composed of highly variable 56 component parts (Gonzalez and Loreau, 2009). For example, compensatory dynamics stabilize 57 productivity when different populations have offsetting fluctuations, such that increases in the 58 abundance or biomass of one or more species are accompanied by corresponding decreases in 59 others (e.g., Schindler, 1990, Frost et al., 1995, Bai et al., 2004, Hallett et al., 2014. Conversely, 60 when species increase or decrease together (e.g., when species share responses to an external driver), 61 the resulting synchrony increases aggregate community variability and has been used as a proxy for 62 instability (e.g., Houlahan et al., 2007, Keitt, 2008, Feio et al., 2015, Ma et al., 2017. 63 Population fluctuations are ubiquitous. Consequently, characterizing patterns of species 64 fluctuations over time and in relation to each other is essential to understand stability. Importantly, 65 population fluctuations can be shaped by a variety of drivers that operate on different timescales. 66 For example, acute or pulsed disturbances are a strong predictor of population dynamics in some 67 systems (e.g., heavy rainfall events or cold winter events regulating insect dynamics; Mutshinda et 68 al., 2011). In other systems, population dynamics are strongly driven by long-term climate cycles 69 (Mantua et al., 1997. As a result, communities may exhibit a timescale 70 specificity such that synchrony (or compensation) among species can occur at some timescales but 71 not others (Keitt, 2008, Downing et al., 2008. Moreover, the co-occurrence of short-and long-72 4 timescale drivers in the same system may result in communities that are synchronous over some 73 timescales but compensatory over others (as idealized in Fig. 1). 74 Evaluating timescale specificity in community variability has a tradition in aquatic ecology 75 (Vasseur et al., 2005, Keitt and Fischer, 2006, Vasseur and Gaedke, 2007, Downing et al., 2008 Keitt, 2008, Vasseur et al., 2014, Brown et al., 2016, but methods used in that literature have 77 typically relied on data-hungry techniques that may not be appropriate for shorter time series. For 78 example, Vasseur et al. (2005) found that phytoplankton in Lake Constance showed compensatory 79 dynamics at a sub-annual timescales, driven by grazing pressure and competition for nutrients 80 during summer and fall, but synchronous dynamics at most other timescales. Downing et al. (2008) 81 described how zooplankton in experimental ponds had synchronous dynamics at quite short 82 timescales (~10-day periods), but had compensatory dynamics at longer timescales (~80-day 83 periods). Overall, synchrony has tended to be more common than compensatory dynamics in 84 freshwater plankton communities, and the timescales at which compensatory dynamics occur 85 appear to be system-specific (Vasseur et al., 2014, Brown et al., 2016 (Schluter, 1984). The variance ratio compares the variance of the aggregate community to the 96 5 expected variance under the assumption of independent population fluctuations across the entire 97 time series. As such, a variance ratio greater than one indicates that populations are generally 98 synchronous, whereas a variance ratio less than one indicates compensation (Peterson, 1975, 99 Schluter, 1984. This simplicity and ease of interpretation has held wide appeal across diverse 100 scientific fields (e.g., Gotelli, 2000, Klug et al., 2000, Houlahan et al., 2007 2009). However, the variance ratio cannot disentangle the timescales at which dynamics occur ( Fig.   102 1), ultimately hindering possibilities for considering mechanisms of community dynamics.

103
Here we develop and apply a new timescale-specific variance ratio appropriate for terrestrial 104 grasslands and other systems with shorter, regularly-spaced time series. In contrast to some of the 105 timescale-specific approaches previously used on plankton data, our techniques provide a formal 106 decomposition by timescale of the classic variance ratio approach, so that appropriately 107 averaging/summing the components of the new approach across timescales recovers the classic non-108 timescale-specific results; one can then quantify the contributions of timescale bands as well as 109 individual timescales. We first develop the theory that underlies our approach and provides 110 timescale-specific measures of community and population variability. Second, we apply this theory 111 to long-term (11-30 years) records of plant community composition at six grasslands across the 112 United States (Hallett et al., 2014). We address the fundamental questions of 1) whether 113 synchrony/compensatory dynamics are timescale-dependent phenomena, and 2) whether 114 compensatory dynamics are rare, compared to synchrony, in grassland systems, as appears to be the 115 case for plankton systems (Vasseur et al., 2005, Keitt and Fischer, 2006, Vasseur and Gaedke, 2007, 116 Downing et al., 2008, Keitt, 2008, Vasseur et al., 2014, Brown et al., 2016. We aim to demonstrate 117 our timescale-specific approach in a way that deepens our understanding of grassland community 118 dynamics as developed using the traditional variance ratio metric (Hallett et al., 2014). Finally, we 119 6 provide a software package for the R programming language to facilitate wider adoption of our 120 techniques. (theorem 2 of Appendix S2; Peterson, 1975, Schluter, 1984. Thus values (respectively, < 1) 141 7 correspond to greater (resp., lesser) community variability than would be expected if dynamics of 142 different taxa were independent, reflecting synchronous (resp., compensatory) dynamics. To illustrate, we apply the theory to the artificial example of Fig. 1. The timescale-specific 187 variance ratio was greater than one for , and less than one for , capturing the 188 deliberately constructed synchronous ( ) and compensatory ( ) processes at these 189 timescales (Fig. 3a). The non-timescale-specific variance ratio , being a weighted

201
We applied our timescale-specific variance ratio to six long-term grassland plant community 202 composition datasets from sites distributed throughout the US (Table S1;  we will refer to these as "nominal" p-values. 246 We used a similar approach to test whether short or long timescales tended to contribute  plots. Using the classic approach, dynamics at this plot would be considered compensatory, as the 280 classic variance ratio was 0.457 (which is less than 1), and correspondingly CV 2 com (value 0.047) 281 13 was less than CV 2 com_ip (value 0.103). However, when we decomposed variability by timescale 282 bands, the short and long timescale bands showed opposite patterns of synchrony and compensation 283 (Fig. 4). The weighted average of φ ts (σ) across short timescales was greater than 1, indicating 284 synchronous dynamics on those timescales, but the weighted average of φ ts (σ) across long 285 timescales was less than 1, indicating compensatory dynamics (Fig. 4a). Correspondingly, the 286 average of CV 2 com (σ) was slightly higher than the average of CV 2 com_ip (σ) when these averages 287 were computed across short timescales, but was substantially less when the averages were 288 computed across long timescales (Fig. 4b-d). 289 We applied the same approach to all 150 plots in the six sites (Fig. 5), providing a 290 comprehensive overview of how a timescale approach alters our understanding of compensatory 291 dynamics in these systems. Our data showed average (across plots) greater aggregate population 292 variability at long timescales than at short timescales at all sites ( Fig. 5g-l), and likewise for 293 community variability (Fig. 5m-r), except for the Kellogg Biological Station LTER (KBS), where 294 the average CV 2 com (σ) did not differ between short and long timescales (Fig. 4n). Thus, long 295 timescales were the primary driver of population variability at all sites and of community variability 296 in five of our six sites; this was not surprising given much prior literature on the commonness of 297 temporal autocorrelation in population dynamics (see Discussion).

298
The degree of synchrony and compensation among species differed substantially by 299 timescale at some but not all sites. Short timescales had larger weighted-average values of φ ts (σ) 300 than long timescales for Jasper Ridge (JRG) and Kellogg Biological Station (KBS; Fig. 5a patterns of compensatory dynamics, timescale-specific patterns of community variability differed 308 between these sites for reasons that our approach reveals. At JRG, population variability was much 309 greater on long than on short timescales (Fig. 5g). Community variability was also higher on long 310 timescales (Fig. 5m), but to a lesser degree than population variability, because species dynamics 311 were more compensatory on long than on short timescales (Fig. 5a). At KBS, population variability 312 was also higher on long than on short timescales, but less markedly so than for JRG (Fig. 5h).

313
Because dynamics were again more compensatory on long than on short timescales for KBS ( Fig.   314 5b), community variability was similar on short and long timescales (Fig. 5n). Thus differences variability, and to about the same extent (compare Fig. 5i,o for HAY, Fig. 5j,p for JRN, Fig. 5k,q   323 for KNZ, and Fig. 5i,r for SEV). There was substantial variation in weighted-average φ ts (σ) among 324 plots at these sites for both long and short timescales. As a result, there were plots within every site 325 that exhibited a greater tendency toward compensation at long compared to short timescales, and  In contrast to the methods which have been used to study plankton communities, our 351 methods offer a formal decomposition of the classical variance ratio approach. The classical 352 16 variance ratio and the community and population variability statistics CV 2 com and CV 2 com_ip can be 353 recovered exactly from our timescale-specific quantities by summing or averaging appropriately 354 over timescales. Thus our approach formally elaborates the classical approach. In developing our 355 theory, we focused on the classical variance ratio (Peterson, 1975, Schluter, 1984. However, a responses to drivers (Ives, 1995