Synchronous dynamics of zooplankton competitors prevail in temperate lake ecosystems

1. Introduction

The extent to which communities of interspecific competitors exhibit synchronized fluctuations over time is intrinsically linked to community stability and moreover may reflect important attributes of the functional diversity of communities [1,2] and of the processes structuring community composition [3–5]. Synchrony and its alternative pattern, compensation (or compensatory dynamics), are mutually exclusive phenomenological features of community dynamics that can arise from many underlying mechanisms [6,7]. Compensation may reflect a reciprocal negative interaction among competitors [4]; however, it also can arise when species respond oppositely to changes in an environmental factor [8], or when they respond to different environmental factors which are themselves negatively correlated through time. By contrast, species that exist in a temporally variable environment must also be qualitatively similar in their ability to tolerate the environment's extremes (e.g. phosphorus is essential for plankton but some species better tolerate phosphorus depletion), suggesting that the environment may select species with similar traits, which in turn predisposes species to exhibit synchronized dynamics [1]. The extent of synchronous and compensatory dynamics in competitive communities may therefore reflect differences in the relative strength of environmental forcing and competition [1,3].

Synchronous and compensatory dynamics are also of fundamental interest in community ecology because they regulate the relationship between population and temporal community variability. At the scale of populations, temporal variability is an indicator of stability which is often linked to extinction probability [9,10]. Similarly, temporal variability at the aggregate or community level, ‘community variability’ [7,11], is of fundamental interest as an index of community stability [12,13] and is of applied interest as a potential leading indicator of nonlinear changes in ecosystem state [14,15]. When species contribute additively to a particular community attribute (e.g. total biomass), the temporal variability of the attribute is given by the sum of all terms in the community variance–covariance matrix. All else being equal, communities with species dynamics that negatively covary over time (compensatory dynamics) will exhibit lower community variability than communities whose covariance terms are positive on average (synchronous dynamics) [7]; however, the same processes that generate negative covariation can also destabilize communities by increasing species fluctuations [16].

Many studies have measured the synchrony of competitors in experiments and nature [3,7,17–20] and some have attempted to infer the relative contributions of competition and environmental factors in these communities [1,3,20]. For example, using annual census of terrestrial plants and animals, Houlahan et al. [3] found that species are on average synchronized and suggested that environmental drivers were relatively more important than competition for shaping community dynamics. Although appealing, the argument that synchronous dynamics result when environmental drivers are relatively more important than density-dependent biotic feedbacks is challenging to verify, because the drivers leading to synchronous and compensatory dynamics need not necessarily be in conflict. Synchronous and compensatory dynamics may operate on different timescales [20,21] or during different periods [21] in which case the frequency and temporal extent of sampling largely dictates which drivers are perceived to be the most critical. For example, if competing species are commonly and uniformly affected by environmental variation, it simply generates a temporally changing arena in which the outcomes of competition are enacted. If species are sampled once annually as in [3], the signature of competition may be entirely obscured by the synchronizing effect of ‘good’ versus ‘bad’ years (e.g. drought versus wet), even if it is an important driver of species' relative fitness and persistence [22]. Furthermore, it has been shown that embedding competitors in a food web context can cause them to respond synchronously to environmental disturbance even when the underlying dynamics are compensatory [23], highlighting the challenge of identifying the mechanisms underlying dynamics. Measuring and interpreting synchrony and compensation thus requires, at a minimum, refined estimators of synchrony and compensation across different timescales [17–21].

Freshwater plankton communities have been a focus of recent work on the scale-dependence of compensatory and synchronous dynamics. In phytoplankton and zooplankton communities in Lake Constance (Bodensee), compensatory dynamics have been shown to occur regularly at scales ranging from six months to 1 year while synchronous dynamics dominate at most other scales [17,21]. Synchrony of zooplankton communities in Little Rock Lake intensified at the annual scale during experimental acidification, with no apparent alteration of covariation patterns at other scales [18,19]. Freshwater plankton are ideal for analyses of community variability due to their short generation times and importance in long-term monitoring programmes worldwide. In temperate regions, lakes experience environmental fluctuations on timescales ranging from short-term weather fluctuations (e.g. solar irradiance and wind-driven mixing) to long-term trends in nutrient loading, pH, temperature and other factors.

In addition to timescale, communities may transition between synchrony to compensation over time at fixed scales. In Little Rock Lake, experimental acidification altered the expression of synchronization at the annual scale in a zooplankton community [18] and elsewhere, natural and anthropogenic changes in long-term nutrient loads have been shown to impact community composition [24,25]. Within the growing season, the intensity of competition among zooplankton can vary greatly as nutrients become depleted and well-edible algae are replaced by colonial forms [26]. Transitions between synchrony and compensatory dynamics may be common at a variety of scales and yield important insight into the mechanisms structuring the distribution of temporal dynamics of species in competitive communities; however, detecting such transitions requires advanced statistical machinery capable of resolving patterns in both time and scale.

Here, we quantify patterns of dynamics within communities in long-term time series (mean 17 years) of crustacean zooplankton species from 58 study sites originating from 52 North American and European temperate lakes and two sub-tropical lakes (figure 1). We employ a wavelet analysis, which estimates the amount of variation in a time series attributable to a particular frequency (scale) at a particular point in time [27], to generate a scale and time-resolved binary metric of synchronous/compensatory dynamics among zooplankton assemblages within each of our study sites. Our extensive dataset allows us to then overlay this information from various lakes to determine the relative prevalence of synchronous and compensatory dynamics over a wide range of temporal scales (monthly to decadally) and across a large span of time (1971–2008) in a large sample of 58 study sites (54 lakes). The objective of this work is to (i) provide a critical assessment and frame of reference for the prevalence of synchronous and compensatory dynamics in competitive communities and (ii) to highlight potentially important scales and times at which transitions in community dynamics occur and may help to identify the mechanisms structuring competitor dynamics.

• Download figure
• Open in new tab
• Download PowerPoint

2. Material and methods

(b) Wavelet analysis

To provide a scale- and time-resolved metric of synchronization, we began by calculating the continuous wavelet transformation W_k,l of the time series x_k,l for each zooplankton species (k) in each lake (l) as the convolution of x_k,l with the conjugate of the scaled and translated Morlet wavelet ψ(·) [27,32]

2.1

In expression (2.1), t_i represents the set of Julian dates on which x_k,l was sampled and N is the total number of samples of x_k,l. The parameters n and s represent the time and scale localization, respectively, of the Morlet wavelet, which is the product of a Gaussian distribution and a complex waveform, . We set the wavenumber of the Morlet wavelet, ω₀, which controls the number of oscillations within the effective width of the Gaussian distribution, equal to 6 as in previous studies [18,27].

In order to facilitate comparison across our study sites, we chose standardized arrays of times (n) and scales (s) at which to sample the wavelet transformation and we subset these depending on the limitations of each time series. In the temporal dimension, we sampled the wavelet on every 10th Julian day provided that it fell within the ‘cone of influence’, which is defined as the range of points further than from the start and end of the time series [27]. In the scale dimension, we generated an array of scales according to

2.2

where δ = 0.1 and j is a sequence of integers. We set the minimum value of j such that the smallest scale was not less than three times the average time between successive samples in the time series. We set the maximum value of j such that the largest scale did not exceed half the total period of sampling in the time series; however, only scales smaller than this threshold included samples that were admissible based upon the ‘cone of influence’.

Most of our time series were unevenly sampled through time, potentially producing combinations time and scale which were unreliable within the boundaries defined above. For example, if sampling was discontinued during winter (which is common in temperate lakes), there may be few or no data points in the neighbourhood of the wavelet at shorter scales during winter. To deal with this, we further discarded any wavelet samples for which fewer than 10 observations occurred within of the centre of the Morlet wavelet.

(c) Measuring synchrony

We first quantified the relative amplitude of zooplankton crustacean species dynamics within each of l sites using the localized wavelet modulus ratio [18]

2.3

where denotes the complex modulus and is a Gaussian localization function in time n. The localized wavelet modulus ratio computes the ratio of the actual amplitude of community fluctuations against the cumulative amplitude of population fluctuations at each time and scale. The numerator in equation (2.3) uses the modulus of the sum of the wavelet transform; by summing wavelet coefficients prior to computing the modulus, the realized amplitude is discounted by species whose dynamics have opposite phases. The denominator, the sum of species' moduli, is the maximum possible amplitude given a perfect alignment of phases. Thus, the wavelet modulus ratio is bounded between values of zero and one. A value of zero indicates that compensatory dynamics result in a community variance equal to zero, whereas a value of one indicates perfect synchronization of population dynamics and a maximum value of the community variance. The modulus ratio is akin to a family of multivariate indices of synchronization [7] that measure the sum or mean of pairwise species covariances (or here co-amplitudes) by comparing realized variances or amplitudes of communities to the values that would be achieved under perfect synchronization, thus superseding the need to compute pairwise relationships directly using correlative or co-spectral techniques. Although uneven sampling can generate bias in the estimates of wavelet power (W²), this bias is consistent within each sampling site because sampling occurred for all zooplankton species at the same points in time, making our measure of synchronization (based on relative amplitude) unbiased in the presence of uneven sampling.

To measure synchronization, we generated 1000 null-model outcomes in each site (l), at each scale (s) and time (n) by attributing a random phase-shift on the interval (0, 2π) to each of the k species and recalculating the modulus ratio (2.3) [18]. Because this procedure retains the distribution of amplitudes of variation expressed by different species, it generates an accurate distribution for the value of the local wavelet modulus ratio under the assumption that species’ dynamics are unrelated. Although it is possible to use this distribution of outcomes to statistically differentiate synchronous and compensatory dynamics from independent fluctuations [18], the p-value adjustment required by the combined number of scales, times and lakes to which we would apply this test would require on the order of 10⁶ randomizations at each of the ca 1.1 million combinations of time, scale and sample site, an operation that is computationally infeasible. Instead, we generate a binary response metric at each combination of scale, time and lake by assigning a value of one to cases where p_l(n,s) exceeded the median value of the 1000 bootstrapped realizations and a zero otherwise. We then summarized this metric across our study regions and across the entire dataset to determine the fraction of all lakes demonstrating synchronized dynamics at any particular time and scale.

For the entire dataset, we determined whether the set of lakes contributing information to a particular time and scale were more or less synchronized than expected by chance using a binomial test. Given the large number of tests performed (all combinations of time and scale), we applied a Benjamini–Hochberg–Yekutieli correction to control the false discovery rate [33]. Similar analyses were not possible at the regional scale due to small sample sizes (number of lakes). The wavelet analysis was performed using custom code in R v. 2.6.0 (R development core team), and summary statistics and plots were generated using Mathematica v. 9.0. A summary schematic of our analysis is available (electronic supplementary material, figure S1).

3. Results

Synchronous dynamics are common in freshwater zooplankton communities. Considering all temporal scales (approximately monthly to decadal scales of variation), times (sampling dates ranging from 1971 to 2008) and study sites (58 sites in 54 lakes), we found that zooplankton communities were more synchronized than the median of our null model in 85.8% of these combinations (14.2% showed compensatory dynamics) (figure 2). The vast majority of incidences (79%) were associated with a q-value equal to 1, meaning that the observed dynamics were more synchronized than all 1000 of our null-model realizations (result not shown). Values of the wavelet modulus ratio less than approximately 0.4 were predominantly compensatory, whereas those greater than approximately 0.4 were predominantly synchronous. This implies that communities that were more synchronized than the median of null-model realizations still expressed substantial variation in their extent of synchronization. Strong synchrony (wavelet modulus ratios close to 1) mainly occurred at scales near 1 year; shorter and longer scales tended to exhibit weaker synchronization (cf. figure 2b–d).

• Download figure
• Open in new tab
• Download PowerPoint

We searched for the consistency of dynamics (synchrony or compensation) across study sites by region (using six different study regions) and across the entire 58-site dataset (figure 3). With the entire dataset, we applied a binomial test to determine whether lakes exhibited similar dynamics at each combination of time and scale. At the regional scale, smaller sample sizes precluded the use of statistical tests. Synchronous dynamics are pervasive at the annual scale in all regions and across the entire dataset (p < 0.0066); however, there are no combinations of time and scale at which a significant fraction of sites exhibit compensatory dynamics across the entire dataset (p > 0.05). Rather, times and scales at which compensatory dynamics are predominant tend to differ across regions; while both the Experimental Lakes Area (ELA) and Sudbury area lakes (figure 3a,e) show periods of compensatory zooplankton dynamics at scales greater than 1 year, these are expressed at different times in the two areas and thus are not visible in the aggregate plot (figure 3m). Similar mismatches among regions are also evident at scales smaller than 1 year.

• Download figure
• Open in new tab
• Download PowerPoint

At longer timescales (2–10 years), the entire dataset exhibits a larger fraction of lakes with synchronized dynamics after ca 1984, but not prior (p < 0.0066; figure 3m). This temporal transition is echoed at the regional scale in the Sudbury and Dorset areas (figure 3e,g) suggesting that important and widespread changes that yielded greater long-term synchronization of zooplankton dynamics within lakes occurred at this time.

4. Discussion

Characterizing patterns across all combinations of timescale, sampling date and study site, we found that synchronous dynamics were far more prevalent in lake zooplankton communities than compensatory dynamics. Although we found a range of scales over which a significant fraction of our study lakes demonstrated synchronized dynamics, we found no scales at which a significant fraction demonstrated compensatory dynamics. Furthermore, we detected a significant change over time in the fraction of study sites demonstrating synchronized dynamics at long timescales, an effect that may be caused by recovery from anthropogenic disturbance in a subset of our study lakes.

It is no surprise that we find a predominance of strong synchronization of zooplankton species at and around the scale of 1 year in our analysis (figures 2 and 3) and that this predominance is unchanged through time. Seasonal variation drives strong patterns in both the abiotic (e.g. temperature and nutrients) and biotic (e.g. resources and predators) factors that determine zooplankton abundance [26,34] and the extent to which the prevailing factors seasonally vary in temperate systems leads to synchronous variation among a large fraction of zooplankton species [26]. More interesting, is the rarity of compensatory dynamics at scales below and above 1 year. Below the 1-year scale, we expect resource competition to be ongoing and species may exhibit compensatory dynamics via their temporal partitioning of the growing season [26]; above 1 year and competitive exclusion and species replacement may be important processes driving compensatory dynamics [7].

Given the wealth of theory suggesting that compensatory dynamics should be a common facet of coexisting competitors [4,6–8,35,36], the rarity of compensatory dynamics in our analysis is surprising. A possible explanation is that freshwater zooplankton species are highly differentiated and therefore only weak interspecific competitors. However, this is unlikely given that co-occurring species of zooplankton are known to demonstrate little variation in life-history traits suggesting use of common resources [37] and negative reciprocal interactions among zooplankton taxa common to our study lakes have been previously demonstrated [38,39]. A more likely hypothesis explaining the rarity of compensatory dynamics is that the period during which intense competition occurs may be confined to a rather short portion of the seasonal cycle [21,26]; detecting compensatory dynamics would therefore then require sampling frequencies greater than those employed in most of our study lakes. In addition, our analysis of the complete assemblage of zooplankton species includes interactions other than competition. Larger bodied crustacean zooplankton are known to prey on smaller cladocerans and copepods and although these species make up a relatively small fraction of the zooplankton biomass in most lakes they may have important impacts on the expression of patterns in the dynamics.

Although our analysis revealed compensatory dynamics in 14.2% of all scales and times in our 58 zooplankton communities, we detected no times or scales at which compensatory dynamics co-occurred in the majority of these lakes and very few instances where compensatory dynamics dominated any particular region. This suggests that the mechanisms driving compensatory dynamics may be local in their extent and/or vary substantially in their strength across lakes. Because compensatory dynamics are more likely to occur during periods where resources are limiting for a subset of species [4,7,21], and such periods occur only during certain parts of the seasonal succession [26], differences in the physical and biological characteristics of nearby lakes, and differences in the seasonal forcing regime of spatially distant lakes, may ensure that compensatory dynamics lack coherence across study sites and/or regions. In addition, metacommunity models have shown that locally driven compensatory dynamics can be quickly overcome by environmental- and dispersal-driven synchronization [40] suggesting that compensatory dynamics may be unstable in the presence of other regional factors.

In addition to these broader patterns, our method also revealed scale-specific temporal changes in synchrony and compensation. Although we find no consistent temporal changes at short scales (less than 1 year), we find that zooplankton communities in the Sudbury and Dorset regions became more synchronized at longer scales (more than 1 year) from ca 1984 onwards. This shift towards long-term synchronization of the zooplankton communities in these lakes is potentially a result of biological restructuring and recovery from anthropogenic disturbances such as lake acidification and eutrophication. In the Sudbury region, it has been well established that the growth of acid-sensitive species was substantially impaired prior to 1980, but as pH increased over the mid-1980s the biomass of acid-sensitive species increased [41]. The shift towards synchronized dynamics at long timescales suggests that there was no compensating decline in the biomass of acid-tolerant species in this system, but rather that both acid-tolerant and acid-intolerant species exhibited slowly increasing biomass (albeit at different rates), during the recovery period. In the Dorset area lakes, larger zooplankton body sizes have increased biomass due to a variety of factors stemming from acidification recovery [42]. Similar results have been witnessed in response to eutrophication in Lake Constance from 1920 to 1980 (Southern Germany); during eutrophication, crustacean species that did not go extinct all increased in biomass but at relatively different rates [25]. Previous work has shown that whole-lake acidification induced synchronized dynamics at the annual scale in the zooplankton community in Little Rock Lake, because cold-tolerant species were more greatly affected, thereby selecting differentially for warm-season growers [18]. By contrast, we find that synchronous dynamics are prevalent at the annual scale and across most other scales even in the absence of anthropogenic disturbances.

Our results provide a new and detailed baseline against which the study of zooplankton community dynamics, and the dynamics of competitive communities in general, can be compared and contrasted. For example, the compensatory dynamics described in the acidified basin of Little Rock Lake [18] are an anomaly relative to not only the reference basin of the same lake, but to our entire set of study lakes, suggesting that anthropogenic disturbances may be possible to detect by identifying anomalies against (or within) our reference data. Additionally, contrasting our results to zooplankton communities from lakes without fish, with introduced top predators such as rainbow trout, or where fish are anadromous, could provide important insight into the top-down structuring forces of predators. Furthermore, coupling our results with lake hydrology, water chemistry and climate data may provide insight into the mechanisms linking the environment to the dynamics of competitors.

Although we find a predominance of synchronized dynamics, our results also demonstrate that synchronized communities can have an aggregate (e.g. total biomass) variance that is relatively low. For example, synchrony was more prevalent than compensatory dynamics when the modulus ratio varied between 0.35 and 1.0 (figure 2) but at the lower end of this range the community variance is more than four times smaller than that of a perfectly synchronized community. Thus, the widespread synchronization we detect in our analysis still represents a significant reduction in community variability and therefore may not be as detrimental to ecosystem stability and population persistence as is commonly thought [6,17,43].

Unlike previous work that used methods sensitive to only the dominant mode of variation and produced a measure of variability averaged across time, our analysis parses out the effect of variation occurring at different timescales and during different time periods. Although scale- and time-resolving methods have become commonplace in understanding patterns of species covariation [17,18,20,44], ours is the first analysis to seek patterns across broad spatial and temporal scales. We found evidence for widespread synchronous dynamics at a broad range of times and scales, while compensatory dynamics were relatively rare and incoherently distributed among different times and scales in our study lakes. Although, we found only few regionally consistent changes in zooplankton synchronization over time, this result is itself important—indicating that local drivers of zooplankton dynamics and heterogeneity among lakes are strong enough to obscure larger scale patterns. The rarity of compensatory dynamics, which are assumed to be an essential feature of competitive coexistence, brings into question our understanding of the mechanisms maintaining coexistence of competitors. Future theoretical and experimental work should aim to address how diverse species assemblages persist despite their tendency to display synchronous dynamics across long reaches of time and scale.

Data accessibility

Data can be made available upon request.

Funding statement

This paper emerged as part of the Predicting Ecological Change: Multi-Scale Analysis of Plankton Diversity and Dynamics working group funded by the Canadian Institute for Ecology and Evolution (CIEE/ICEE). This is Kellogg Biological Station contribution no. 1718. R.A. was funded by Limnotip (O1LC1207A; Biodiversa) founded under the FP7 ERA-Net Scheme and the DFG-LakeRisk Project; W.A.N. and B.E.B. acknowledge support of NSERC discovery grants; A.G. is supported by the Canada Research Chair programme; D.A.V. acknowledges support of NSF DEB grant no. 1050803.