# Temporal scale dependent interactions between multiple environmental disturbances in microcosm ecosystems

Aurélie Garnier, Frank Pennekamp, Mélissa Lemoine, Owen L Petchey

► To cite this version:

Aurélie Garnier, Frank Pennekamp, Mélissa Lemoine, Owen L Petchey. Temporal scale dependent interactions between multiple environmental disturbances in microcosm ecosystems. Global Change Biology, 2017, 23, pp.5237 - 5248. 10.1111/gcb.13786 . hal-05129819

HAL Id: hal-05129819

https://hal.science/hal-05129819v1

Submitted on 25 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Temporal scale dependent interactions between multiple environmental disturbances in microcosm ecosystems

Aurélie Garnier ID

Frank Pennekamp

Mélissa Lemoine

Owen L. Petchey

Department of Evolutionary Biology and
Environmental Studies, University of Zurich,
8057 Zurich, Switzerland

**Correspondence**

Department of Evolutionary Biology and
Environmental Studies, University of Zurich,
Winterthurerstrasse 190, 8057 Zurich,
Switzerland.

Email: aurelie.garnier.research@gmail.com

**Funding information**

Schweizerischer Nationalfonds zur
Förderung der Wissenschaftlichen
Forschung, Grant/Award Number:
31003A_159498, PMPDP3_151361,
PMPDP3_161858; University Research
Priority Program: Global Change Biology

**Abstract**

Global environmental change has negative impacts on ecological systems, impacting the stable provision of functions, goods, and services. Whereas effects of individual environmental changes (e.g. temperature change or change in resource availability) are reasonably well understood, we lack information about if and how multiple changes interact. We examined interactions among four types of environmental disturbance (temperature, nutrient ratio, carbon enrichment, and light) in a fully factorial design using a microbial aquatic ecosystem and observed responses of dissolved oxygen saturation at three temporal scales (resistance, resilience, and return time). We tested whether multiple disturbances combine in a dominant, additive, or interactive fashion, and compared the predictability of dissolved oxygen across scales. Carbon enrichment and shading reduced oxygen concentration in the short term (i.e. resistance); although no other effects or interactions were statistically significant, resistance decreased as the number of disturbances increased. In the medium term, only enrichment accelerated recovery, but none of the other effects (including interactions) were significant. In the long term, enrichment and shading lengthened return times, and we found significant two-way synergistic interactions between disturbances. The best performing model (dominant, additive, or interactive) depended on the temporal scale of response. In the short term (i.e. for resistance), the dominance model predicted resistance of dissolved oxygen best, due to a large effect of carbon enrichment, whereas none of the models could predict the medium term (i.e. resilience). The long-term response was best predicted by models including interactions among disturbances. Our results indicate the importance of accounting for the temporal scale of responses when researching the effects of environmental disturbances on ecosystems.

**KEYWORDS**

environmental changes, microbial aquatic system, multiple drivers, predictability, resilience, resistance, return time, temporal scales

## 1 | INTRODUCTION

Global environmental change is known to affect ecological systems in harmful ways and threatens the stable provisioning of functions, goods, and services that ecosystems provide (Chapin III et al., 2000;

Daily et al., 2000). Among the most important types of global environmental change are habitat loss and fragmentation, overexploitation, invasive species, and coextinctions, aptly depicted as 'the four horsemen of the ecological apocalypse' (Diamond, Ashmole, & Purves, 1989). Added to this evil "quartet" nowadays is climate

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.

© 2017 The Authors. *Global Change Biology* Published by John Wiley & Sons Ltd.
change, and to make things worse, there exists the potential for synergies between co-occurring environmental changes (Brook, Sodhi, & Bradshaw, 2008). Synergies would exacerbate pressure on natural ecosystems and, if they are difficult to predict, could lead to "ecological surprises", with potentially severe and irreversible consequences (Carpenter, Fisher, Grimm, & Kitchell, 1992; Heugens, Hendriks, Dekker, van Straalen, & Admiraal, 2001; Brook et al., 2008; Griffen & Drake, 2008; Holmstrup et al., 2010; Mantyka-pringle, Martin, & Rhodes, 2012). On the other hand, antagonistic interactions mitigate each other's effect (Folt, Chen, Moore, & Burnaford, 1999; Brook et al., 2008; Crain, Kroeker, & Halpern, 2008). A special case of such antagonistic interactions is when the combined effect of multiple environmental disturbances is equal to the largest effect of any of the disturbances when they occur in isolation (Sala et al., 2000; Brennan & Collins, 2015).

The presence and strength of interactions among multiple environmental disturbances can have large effects on predictions. For example, Sala et al. (2000) compared the future global distribution of biodiversity for scenarios with different assumptions about how multiple environmental disturbances combine. The biome in which biodiversity was most threatened depended greatly on whether one assumed additive/synergistic or dominant combining of the effects of multiple environmental disturbances. The study concluded that the most plausible scenario for the future would be between the additive and synergistic hypothesis, and highlights the importance and priority of research about how multiple environmental disturbances combine.

Although numerous conceptual frameworks for discriminating between synergistic and antagonistic effects exist (Piggott, Townsend, & Matthaei, 2015), experimental approaches that manipulate environmental disturbances in a factorial manner, which allows to rigorously test for interactive effects, are still rare (but see Doyle, Saros, & Williamson, 2005; Christensen et al., 2006; Brown et al., 2012; Griffiths, Warren, & Childs, 2015). Often these studies concern only interactions between two environmental factors, and evidence regarding the occurrence and types of interactions is mixed (Darling & Côté, 2008; Jackson, Loewen, Vinebrooke, & Chimimba, 2016). To evaluate the reliability of the scenarios for management decisions, we urgently need to understand how important interactions are and how well we can forecast with models that neglect interactions (Côté, Darling, & Brown, 2016).

Aquatic systems are particularly vulnerable to environmental changes due to their importance for and proximity to human settlements (Jenny et al., 2016). Land use changes, invasive species, climate change, nitrogen deposition, and atmospheric carbon dioxide are considered major threats for aquatic organisms (Carpenter et al., 1992; Sala et al., 2000; Stendera et al., 2012). Although many of these environmental disturbances were studied individually to understand their consequences at different levels of ecological organization, studies investigating their effects in combination are rare (Jackson et al., 2016).

A key indicator of the health of aquatic ecosystems is dissolved oxygen (DO) (Walker, 1979; Wetzel & Uchman, 2001; Hanson, Carpenter, Armstrong, Stanley, & Kratz, 2006). Dissolved oxygen is a

measure of ecosystem productivity that integrates production and respiration across trophic levels and thus estimates a whole-ecosystem response. Change in dissolved oxygen is hence a functional metric that provides the net effect of different processes. Biologically driven processes provide an integrative measure of the ecosystem functioning (Webster & Benfield, 1986) over time and across organisms at different organizational levels. Because functional metrics are independent on the identities of the species in a community, they provide a more generalizable picture than the specific structure of a given community (Denny & Benedetti-Cecchi, 2012). Nevertheless, function influences structure and vice versa, and both should be considered to assess the integrity of an ecosystem as a whole. Community structure and ecosystem functioning are strongly affected by low dissolved oxygen concentration (i.e. hypoxia), which may be insufficient to support heterotrophic organisms (>30% saturation needed; Wu, 2002). Hypoxic environments have become more common in the last three decades (Diaz & Rosenberg, 2008; Diaz & Breitburg, 2009) due to increased human pressure on freshwater ecosystems (Jenny et al., 2016). Temperature, among other factors, affects dissolved oxygen directly and indirectly by affecting its solubility (Garcia & Gordon, 1992) as well as the physiology of organisms (Brown, Gillooly, Allen, Savage, & West, 2004; Savage, Gillooly, Brown, West, & Charnov, 2004). In parallel, nutrient input can trigger bacterial growth (eutrophication), potentially leading to hypoxic condition due to excessive bacterial respiration. Moreover, the interaction of increased temperature and nutrients inputs can intensify hypoxic conditions and ultimately lead to fish extinctions (Moran et al., 2010). Hence, understanding how dissolved oxygen levels respond to (e.g. their resistance to) environmental disturbances, and their recovery (e.g. resilience and return time) from environmental disturbances is important for understanding and predicting responses of species and community composition.

Maintaining stability of ecosystems is often desired, as only stable ecosystems can provide functions and services (Isbell et al., 2015). The ability of ecosystems to buffer disturbances such as induced by global environmental change is therefore an important aspect of ecosystem functioning. Stability may be also a function of time, therefore the temporal scale of the disturbance and the response should be considered (Christensen et al., 2006; Donohue et al., 2016). We chose to apply disturbances in a press manner (rather than pulse) in which disturbances were instantaneously applied and then maintained. We considered three temporal scales of response: the short-term effect of disturbance on dissolved oxygen (i.e. resistance); in the medium term the rate of return of dissolved oxygen to control treatment levels (i.e. engineering resilience); and long-term recovery to control treatment levels ("return time") (Pimm, 1984). We use the term "scale" to describe this variation in the temporal extent over which the responses occur, and thus also the temporal scale of the processes underlying the responses.

We studied the effect of four environmental disturbances, and the direction of interactions among them, on dissolved oxygen availability. As factorial manipulations of environmental disturbances are difficult to achieve in the field, we used an aquatic experimental
system consisting of a community of algae, bacteria, ciliates, and
rotifers (Petchey, McPhearson, Casey, & Morin, 1999; Altermatt
et al., 2015). We selected temperature, nutrients, carbon enrichment,
and light availability as experimental environmental disturbances due
to their relevance for natural aquatic systems (Carpenter et al.,
1992; Piehler, Twomey, Hall, & Paerl, 2004; Llames et al., 2009;
Stanley, Powers, Lottig, Buffam, & Crawford, 2012; Yankova, Villiger,
Pernthaler, Schanz, & Posch, 2016), and manipulated these in a fac-
torial design to detect the effect of potential interactions on DO.

Dissolved oxygen concentration is determined by the action of two biological processes, namely, the respiration of all organisms and the photosynthesis of autotrophs. Effects of the four environmental disturbances on DO will therefore be indirect via effects on community respiration and photosynthesis, and one might expect different effects of each disturbance on each process. For example, carbon enrichment should increase growth, biomass, and therefore respiration at least in the short term, with little effect on photosynthesis, leading to decreased DO. Temperature has stronger effects on respiration than photosynthesis (Yvon-Durocher, Jones, Trimmer, Woodward, & Montoya, 2010) which predicts that increased temperature decreases DO. Shading should decrease DO due to reduced photosynthesis, at least in the short term. In the absence of a quantitative model of the effects of these various disturbances on photosynthesis, respiration, and DO, predictions about how they will interact are difficult to make. Hence, we tested whether multiple environmental disturbances combine according to hypotheses of additivity (combined effect equal to sum of individual effects), synergy (combined effect greater than sum of individual effects), or antagonism (combined effect less than sum of individual effects). A specific form of antagonistic interaction, that of dominance, was also tested (combined effect equal to the largest individual effect).

Dominant and additive combining of multiple types of distur-
bance represents a more predictable situation because then only
information from each individual disturbance is required for predic-
tion. In contrast, interactions between disturbances require addi-
tional, and potentially difficult to obtain, information about the sign
and strength of the interactions. For a particular model of combining
disturbances (e.g. dominant, additive, interactive), we can also ask
how predictability changes with the temporal scale of response.

We expected greater predictability at shorter time scales of
response, and lower predictability at longer time scales of response
due to greater opportunity for indirect effects at longer time scales.
Put another way, direct effects should dominate in the short term,
and direct effects should be more additive/dominant in their combi-
nations, with subsequently greater predictability. In the longer term,
indirect effects, such as those mediated via changes in environmen-
tal conditions and community composition, create greater opportuni-
ties for interactive combinations of effects of environmental
disturbances. Such contributions of indirect effects to unpredictabil-
ity can cause indeterminacy (i.e. unpredictability) of theoretical per-
turbation experiments and ecological surprises (Doak et al., 2008).

2 | MATERIALS AND METHODS

2.1 | Experimental system

Experimental microcosms were sterile 250 ml glass jars containing 100 ml of Protozoan Pellet Medium (PPM) (Lawler & Morin, 1993; Altermatt et al., 2015). Media consisted of 0.28 g of crushed Protozoan Pellets (Carolina Biological Supply Co., Burlington, NC, USA) in 1 l of Chalkley's medium. Protozoan pellets provide an organic food source (nutrient and carbon) for bacteria and protists (Kaunzinger & Morin, 1998). Two additional wheat seeds provided a slow-release nutrient source. Microcosms were placed randomly in six temperature and light-controlled incubators with a 16–8 hr light–dark cycle, at an intensity of 5000 lux during light phase.

2.2 | Microbial aquatic community

Our aim was to assemble a moderately complex microbial commu-
nity with multiple species in multiple trophic groups, so a range of
ecological processes were occurring. This was accomplished by
assembling a community initially composed of two species of bacte-
ria (Serratia fonticola and Bacillus subtilis – generally used in labora-
tory experiments with ciliates cultures (Altermatt et al., 2015)),
although bacterial composition was not subsequently controlled and
was likely higher as the sampling was conducted in nonsterile envi-
ronment, four species of algae (Chlamydomonas reinhardtii, Scenedes-
mus quadricaula, Staurastrum gracile, and Desmidium swartzii), one
species of rotifer (Rotifer sp.), and twelve species of ciliates; one was
algivorous (Nassula aurea), five were bacterivorous (Tetrahymena
thermophila, Colpidium striatum, Paramecium caudatum, Blepharisma
japonicum, and Euplotes sp.), and six were omnivorous (Euplotes
daidaleos, Frontonia sp., Paramecium bursaria, Stentor coeruleus, Dilep-
tus anser, and Actinophrys sol; the last two have a preference for cili-
ates, flagellates, amoebae, and rotifers). Based on results of previous
experiments, extinctions of some species will have happened, partic-
ularly at the higher trophic levels, leading to a community with more
species at lower trophic levels and fewer at high trophic levels (we
did not have access to species composition data when this article
was prepared).

Before the experiment, all species were cultured in monoculture
in 0.28 g L⁻¹ PPM at 20°C. At day 0, all species were combined
with different volumes according to their trophic position (10 ml for
each algae and bacteria species, 2 ml for each ciliate and rotifer spe-
cies), and topped up with 13.8 ml of 0.28 g L⁻¹ PPM and 100 µl of
each nutrient solution (NH₄Cl and KH₂PO₄ in mg L⁻¹) to a total of
100 ml per microcosm. To assure the presence of predators in the
system, five individuals of *Stentor coeruleus*, *Dileptus anser*, and *Acti-
nophrys sol* were added to all microcosms the day before the pertur-
bation treatments. Samples were taken from which we aimed to
estimate the abundance of each species; analyses concerning this
data will appear in a subsequent publication (although we do here
report some preliminary bacterial abundance data in the Figure 5).
## 2.3 | Experimental design

The experiment was four-way fully factorial with two levels of each treatment, with 6 replicates of each of the 16 treatment combinations making for a total of 96 experimental ecosystems. This constitutes a quite large and time consuming experiment, and with available resources we could not include more than two levels in each treatment.

An important aspect of this design was the choice of the two levels of each treatment, and we provide justification of these choices below. Nevertheless, it is important to note considerable variability in the predicted real change in these environmental disturbances; this variability results from uncertainty about what is likely, but also from variability through spaces (e.g. some locations likely to be warmer than others). With such variability, choosing most realistic treatment levels for any single environmental disturbance is somewhat arbitrary.

Perhaps more important than individual treatment levels are their relative levels. If we unwittingly made one of the four treatments large in magnitude and the other three small, we could accidentally favour the dominance hypothesis, for example. To avoid this problem, we tested several levels of each environmental change in preliminary experiments. Temperature was held constant at 20°C or increased to 25°C during manipulation using temperature-controlled incubators. These temperatures were chosen to align with interannual variation in summer water temperature in ponds and lakes (Moore, Folt, & Stemberger, 1996; Jankowski, Livingstone, Bührer, Forster, & Niederhauser, 2006; Yankova et al., 2016). Moreover, this increased temperature (+5°C) falls within the projection of increases in surface water (A1FI scenario, IPCC, 2007). These temperatures may or may not translate into large effects on physiological rates, depending on the temperature response curves of the diverse species in our communities. Lack of knowledge of many of these response curves limited our ability to use such information when deciding treatment levels.

If light availability is not directly a driver of global change, increased dissolved organic matter due to runoffs can result in an increase in turbidity and therefore a decrease in light availability (Anneville et al., 2002). Light availability has been shown to affect phytoplankton photosynthesis (Kirk, 1983), turbidity (Llames et al., 2009), and phytoplankton biodiversity (Stomp et al., 2004). A reduction of 25% light availability was previously found to decrease plankton abundance and to increase respiration rate (Llames et al., 2009). Therefore, we chose to reduce the light availability by 30% using shade cloth around the microcosms to assure a relative magnitude of the light disturbance in comparison to the other disturbances.

Human activities (e.g. agriculture) have resulted in increased loading of nutrients in freshwater systems that affect community structure and function (Smith, Tilman, & Nekola, 1999; Piehler et al., 2004). If nitrogen and phosphorus are limiting resources for primary production, high nutrient inputs can lead to eutrophication of natural systems (Carpenter et al., 1998). Oligotrophic lakes are characterized by high nitrogen:phosphorus (N:P) ratio, whereas eutrophic lakes have a lower N:P ratio, generally below 20:1 (Wetzel, 1983; Downing & McCauley,

1992; Stets & Cotner, 2008; Kratina, Greig, Thompson, Carvalho-Pereira, & Shurin, 2012). And even within a lake, the N:P ratio can vary seasonally from 8 to 60 (Kolzau et al., 2014). We prepared nitrogen and phosphorus solutions to a ratio of 40:1 using NH₄Cl (at 0.460 mg L⁻¹, corresponding to 1.7576 mol L⁻¹ of N) and KH₂PO₄ (at 0.010 mg L⁻¹, corresponding to 0.0439 mol L⁻¹ of P), respectively. We manipulated N:P ratios by increasing the amount of phosphorus (0.027 mg L⁻¹ KH₂PO₄, corresponding to 0.1185 mol L⁻¹ of P) with the same amount of nitrogen, resulting in N:P = 15:1.

Similar to inorganic nutrients (nitrogen and phosphorus), dissolved organic carbon has increased in aquatic ecosystems due to anthropogenic pressure (Stanley et al., 2012; Williams et al., 2016). Carbon enrichment consisted of 0.56 g.L⁻¹ of PPM, and the low concentration was 0.28 g L⁻¹ of PPM (Lawler & Morin, 1993). Importantly, this approximately factor two difference between levels of the carbon enrichment treatment is small relative to many experimental manipulations, which often cover orders of magnitude (Kaunzinger & Morin, 1998), and is small compared to the differences that can occur as a result of inputs into naturally occurring water bodies (~0.1 PgC yr⁻¹; Regnier et al., 2013). Every 3 days, we removed 5 ml of medium from each experimental unit, and replaced with 5 ml of specific medium for the treatment.

During the first week all communities experienced control treatment levels. On the eighth day, we applied a full factorial combination of four press disturbances (temperature, nutrient, carbon enrichment, and light). Responses to the perturbations were monitored until dissolved oxygen had returned to control levels in a large majority (90%) of the replicates (this was achieved by 16 days).

## 2.4 | Quantification of dissolved oxygen content

The dissolved oxygen (DO) saturation was measured daily at the end of a light period of 16 hr using a noninvasive method called chemical-optical sensor (Fibox 4 trace, PreSens, Germany; Altermatt et al., 2015). We assessed the net effect of respiration and photosynthesis on DO. Note that 100% DO saturation corresponds to an oxygen partial pressure of 21%.

## 2.5 | Responses variables

We quantified the treatment effects on DO at three time scales: a short-term response (resistance); a medium-term response (resilience); and a long-term response (return time) (Figure 1). Resistance was the effect observed within 3 days after the perturbation (Pimm, 1984). To measure resistance, we determined the maximum difference between DO in a replicate treatment and average DO across the control replicates. We chose a period of 3 days because visual inspection of the DO time series showed this was long enough to always include the minimum DO caused by the environmental change treatment.

Resilience, in this study, is considered as the rate of recovery following a perturbation (Pimm, 1984), also known as "engineering resilience" (Holling, 1996). Theoretically, the resilience is measured as
the asymptotic rate of return (Arnoldi, Loreau, & Haegeman, 2016). Empirical measures of resilience are challenging and less well defined. We estimated the resilience as the rate of change in log difference between a treatment replicate and the average of the control replicates from the day at which DO reached the maximum displacement; this excluded the possibility for system reactivity (Neubert, Caswell, & Solow, 2009) to interfere with our measure of resilience. Calculating the log difference is equivalent to calculating the rate of relative return, rather than absolute rate, rendering the resilience at least conceptually independent of resistance (Figure 1). The rate of change was estimated by fitting a polynomial of degree three (cubic regression) as this was well supported by the data. Resilience was the first derivative of this polynomial after the system started to return towards DO levels in controlled microcosm (Figure 1).

The return time was estimated as the amount of time taken for DO in a perturbed treatment to recover to the level in control treatments. In practice, this requires accounting for variability in DO among and within control replicates, accomplished by calculating a 95% confidence interval for control DO levels. We also needed to account for variability in DO levels among treatment replicates, again accomplished by calculating a 95% confidence interval around the order-3 polynomial fitted to the return dynamics (the same as used to calculate resilience). Mean return time was the time it took for the mean DO of treatments to fall within the 95% CI of the control, and lower and upper bounds on the return time were when the upper and lower bounds of DO from the treatment 95% CI first fell within the 95% CI of control (Figure 1). If the DO of treatments did not return to within control levels during the experiment, return time was right censored (i.e. the event was not observed at the end of the experiment). In the theoretical setting of exponential return, resilience (rate of exponential return) is the inverse of time to return

(Pimm, 1984). We did not observe such return dynamics, and analysed resilience and return time independently as they were not correlated (Pearson's r = 0.057; t = 0.5521, df = 93, p = .5822). Furthermore, it is important to note that responses were always relative to average control levels to account for any directional changes in control treatments.

## 2.6 Statistical analyses

In a first step, resistance, resilience, and return time were analysed separately to test for the presence and direction of interactions between environmental change treatments. Explanatory variables were the four treatments: temperature (T), nutrients (N), light (L), and carbon enrichment (C), each with two levels (control and perturbed) as well as all high-level interactions (Table 1). Resistance and resilience were analysed with a linear model using a normal error distribution with the package stats (R Core Team, 2016), and return times were modelled using survival analysis with the package survival (Therneau, 2015). The shape parameter of the survival analysis was analysed while the scale parameter was fixed at 1 to avoid lack of convergence. All models were examined visually for the homogeneity of variances and normality and found to follow model assumptions. The significance of effects was tested using two-tailed Type III F- or χ-test on the global model using maximum likelihood with the package car (Fox & Weisberg, 2011).

Then we tested the effect of the number of perturbations on each response variable using a mixed linear model with the number of perturbations as an explanatory variable, and the treatments and the replicates as random effects with the package lme4 (Bates & Mächler, 2015). To correct for the fact that a particular treatment could be involved in different combinations, the overlap between treatments was calculated according to Brennan and Collins (2015)

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Illustration of how ecological stability variables, i.e. resistance, resilience, and return time were measured. The red vertical line shows the time of the disturbance(s). The blue line and points show the dissolved oxygen levels in one treatment replicate (here carbon enrichment). The black line and white points show the six control replicates. The shaded regions show the 95% confidence interval of the control (grey) and treatment (blue) time series</figcaption></figure>

<table><caption>**TABLE 1** Models fit to the observed data, with the dataset used to parameterize the models and the numbers of parameters estimated. Notice that the full dataset corresponded to 95 microcosms (instead of 96) because one microcosm (TCL treatment) was removed due to an error in the treatment application</caption><thead><tr><th>Hypotheses/ models</th><th>Formula</th><th>Number of microcosms used to parameterize the models</th><th>Number of model parameters</th></tr></thead><tbody><tr><td>Dominant</td><td>~ Largest main effect</td><td>Subset (30)</td><td>5</td></tr><tr><td>Additive</td><td>~ T + N + C + L</td><td>Subset (30)</td><td>5</td></tr><tr><td>Interactive (two-way)</td><td>~ T + N + C + L + all two-way interactions</td><td>Subset (66)</td><td>11</td></tr><tr><td>Interactive (three-way)</td><td>~ T + N + C + L + all two and three-way interactions</td><td>Subset (89)</td><td>15</td></tr><tr><td>Interactive (four-way)</td><td>~ T + N + C + L + all 2, 3 and 4-way interactions</td><td>95</td><td>16</td></tr></tbody></table>
and used as a covariate. The shape of the relationship between the ecosystem response and the number of perturbations can inform the interplay among environmental changes. Additive effects would lead to a linear relationship between the number of perturbations and the ecosystem response while interactive (synergic or antagonist) effects would lead to a nonlinear relationship. Finally, the ecosystem response should follow a bimodal distribution when an environmental change would dominate (i.e. with and without the dominant disturbance). Therefore, we first tested the significance of the quadratic term of a polynomial regression to evaluate whether the relationship between the response and the number of perturbations was linear or quadratic. The bimodality of the distribution was investigated using the model including the number of perturbations in interaction with a categorical variable describing the presence of the dominant disturbance in the treatments. The significance of each effect was tested using a two-tailed Type I F-test with the package lmerTest (Kuznetsova, Brockhoff, & Christensen, 2016) using the Satterthwaite approximations for denominator degrees of freedom. Examination of full (linear and mixed) models and backward procedures (first removing the interactions) gave the same results.

In a second step, we examined the predictive power of three groups of hypotheses: dominant, additive, and interactive (Table 1) using a twofold cross-validation method. This involved fitting multiple linear models to the first half of the experimental data (test dataset) and then measuring how well the models predicted the second half (validation dataset). We used the adjusted R² as measure of predictive power. For resistance and resilience, the predictions were the means estimated by the linear models, whereas for return time, the predictions were the time corresponding to 50% of the survival curve. Models varied in the combinations of explanatory variables included, corresponding to (a) a nonadditive effect of treatment with only the largest main effect (dominance), (b) one model of additive effects of treatments, and (c) three models of interactive effects of treatments with up to two-, three-, or four-way interactions (the full model). We only used the minimum data required to parameterize each model: e.g. the additive model used only the experimental data for the main effects, without any interaction treatment combinations (Table 1). To examine the importance of the carbon enrichment treatment, as it appeared to be strong relative to the other treatments, we repeated the entire analysis procedure for the subset of the data that corresponds to performing the three-way factorial experiment with only the temperature, nutrient, and light treatments.

Testing models on data to which they were fitted was likely to yield overoptimistic predictive power (overfitting), therefore we fitted the models to data from three replicates of each treatment combination, and compared their predictions with the other three (or two for TCL). Notice that the full dataset corresponded to 95 microcosms (instead of 96) because one microcosm (TCL treatment) was removed due to an erroneous treatment application. To obtain confidence intervals (95%) of predictive power, we repeated the entire process 1000 times, with replicates randomly assigned to training and test datasets. All statistical analyses were performed in R (R Core Team, 2016).

## 3 RESULTS

At the start of the experiment, DO dynamics were similar across all replicates, including the control, increasing from about 60 to 100% (Figure 2, $F_{15,79} = 0.7439$, *p* = .733), and then little directional change in control replicates from when the treatments were applied. The DO sometimes exceeded 100% saturation, most likely due to production of oxygen by algal photosynthesis at a faster rate than loss by respiration. As expected, carbon enrichment and shading decreased DO, while increased temperature and changed nutrient ratios had no apparent effect (Figure 3a).

Resistance was lowest in the carbon enrichment perturbation (C), with saturation decreasing to around 17% within 2 days of the press perturbation (Figures 2, 3a, Table 2, mean effect: -82.9% oxygen saturation with 95% CI=[-99.5, -66.3]). Light availability had a relatively smaller negative short-term effects on DO (Figures 2, 3a, with mean effects of -30.9% oxygen saturation (95% CI=[-47.5, -14.2])) and nutrients with an even smaller effect size with 95% confidence interval including zero: mean = -3.8% oxygen saturation (95% CI=[-20.3, 12.7]). In contrast, temperature had a small positive short-term effect on DO (Figure 3a, 2.7% oxygen saturation (with 95% CI=[-13.8, 19.2))). There were no significant interactions between the four disturbances affecting the short-term response to DO (Table 2).

Resilience, showing the medium-term response, was mainly negative, meaning that DO was returning towards initial levels (Figures 2, 3b). Resilience did not differ among treatments except for carbon enrichment (Table 2, Figure 3b), which caused a faster rate of recovery (more negative values) (Table 2, -0.36% oxygen saturation per day with 95% CI=[-0.72, -0.01]).

Observed return time (Figure 3c), showing the long-term response, was analysed with survival analysis (survival curves shown in Suppl. Figure S1). An increase in the shape parameter corresponded to delayed recovery. Return time did not differ for disturbances applied independently (Table 2), but some positive, two-way interactions were significant (T:C and N:L) increasing recovery time.

The relationship between response of dissolved oxygen and the number of perturbations was linear for resistance (quadratic term: $F_{1,12} = 0.03$, *p* = .876; linear term: $F_{1,12} = 6.14$, *p* = .029) and return time (quadratic term: $F_{1,11.8} = 2.76$, *p* = .123; linear term: $F_{1,11.8} = 17.81$, *p* < .01) but not for resilience (Figure 3e, quadratic term: $F_{1,12} = 3.05$, *p* = .106; linear term: $F_{1,12} = 0.86$, *p* = .372). Only the presence of the dominant disturbance (i.e. carbon enrichment) in the treatment affected resilience (dominant disturbance effect: $F_{1,85} = 6.7$, *p* = .011; number of perturbations: $F_{1,85} = 0.91$, *p* = .343; interaction: $F_{1,85} = 1.68$, *p* = .198). In contrast, the number of perturbations as well as the presence of the dominant disturbance (i.e. carbon enrichment) had a significant effect on resistance of DO (Figure 3d, dominant disturbance effect: $F_{1,11} = 155.02$, *p* < .001; number of perturbations: $F_{1,10.9} = 85.23$, *p* < .001; interaction: $F_{1,11} = 0.44$, *p* = .520). The number of perturbations only affected the return time (Figure 3f, dominant disturbance effect: $F_{1,10.9} = 2.11$, *p* = .174; number of perturbations: $F_{1,10.8} = 16.66$, *p* < .005;
<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Time series of dissolved oxygen concentration for all 16 treatments and replicates within each treatment. The red line shows the time of the disturbances (temperature [T], nutrients [N], light [L], and carbon enrichment [C]). The grey dashed line represents 100% oxygen saturation. Colours indicate the number of combined perturbations (yellow = control, green = 1 perturbation, turquoise = 2, blue = 3, purple = 4)</figcaption></figure>

interaction: $F_{1,10,9} = 0.73, p = .410). Overall, increasing the number of perturbations decreased resistance of DO and increased return time linearly, whereas it did not affect resilience (Figure 3d–f). The time required to recover increased by about 2 days per additional perturbation. And for comparison, the carbon enrichment treatment decreased the amount of oxygen by an average of 82.9%, while one additional perturbation caused, on average, a decrease of about 15% DO.

We compared predictability among the three temporal scales (Figure 4). Predictability was higher for resistance (adjusted $R^2$ always above 0.5) than it was for resilience and return time (adjusted $R^2$ below 0.5). The 95% confidence intervals of predictive power overlapped for all hypotheses, suggesting that no model performed significantly better for any response. Nevertheless, variation among models was still observed. Including all the interactions among environmental changes explained almost 90% of the variation observed in resistance (median of 88%). The dominant model, despite that it uses, arguably, the lowest number of predictors, already explained 84% of variation. The 95% confidence intervals of resilience included zero for all hypotheses tested. For the return time, while its predictive power did not differ among the hypotheses, the 95% confidence interval of the dominant and additive hypotheses included zero. Including interactions to make predictions increased the predictive power up to 37%.

Analyses that excluded the carbon enrichment treatment had generally lower predictive power, although the ranking of the various

models remained similar (Figure 4). For example, the dominance and interactive hypothesis had similar accuracy for resistance, and including interactions increased predictive power for return time.

4 | DISCUSSION

There is widespread concern that negative effects of global environmental change on aquatic systems will be exacerbated by interactions among multiple environmental changes (Darling & Côté, 2008; Côté et al., 2016). We found scale-dependent importance of interactions between disturbances on dissolved oxygen dynamics. The dominance model (i.e. when the disturbance with the largest effect is used to predict the combined effects of multiple disturbances) was a more parsimonious description of short-term response (i.e. resistance) than the interactive model, and the dominance model was similarly supported in the absence of the large effect of carbon enrichment on DO availability (Figure 4). There was little apparent effect of disturbances in the medium term (i.e. for resilience), and interactions were more important in the long term (i.e. for return time). The predictability of the short-term response was almost 90%, was around 0% in the medium term, and about 40% for long-term response. Our results highlight that importance of interactions may be temporal scale dependent and that models of multiple environmental changes need to account for interactions when making
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 3** Upper panels: observations of each disturbance combination for resistance (a), resilience (b), and return time (c). The dashed lines represent the mean of the control treatment. The colours represent the number of disturbances(s) (as in Figure 2). Lower panels: relationships between the responses and the number of perturbations (d–f). The colours represent the presence (red) and absence (blue) of the dominant driver (i.e. carbon enrichment perturbation “C”). Regressions represent the best model describing the relationship (comparison between linear and quadratic)</figcaption></figure>

<table><caption>**TABLE 2** Analysis of variance (type III, for resistance and resilience) and analysis of deviance (type III, for return time) of four-way linear model on the full dataset. Bold values indicated significant effects (p < .05)</caption><thead><tr><th rowspan="2"></th><th colspan="4">Resistance</th><th colspan="4">Resilience</th><th colspan="3">Return time</th></tr><tr><th>Sum Sq</th><th>df</th><th>F value</th><th>Pr(&gt;F)</th><th>Sum Sq</th><th>df</th><th>F value</th><th>Pr(&gt;F)</th><th>LR Chisq</th><th>df</th><th>Pr(&gt;Chisq)</th></tr></thead><tbody><tr><td>(Intercept)</td><td>255.8</td><td>1</td><td>1.199</td><td>0.277</td><td>4.349</td><td>1</td><td>43.847</td><td><strong>&lt;0.001</strong></td><td></td><td></td><td></td></tr><tr><td>Temperature (T)</td><td>22.0</td><td>1</td><td>0.103</td><td>0.749</td><td>0.155</td><td>1</td><td>1.562</td><td>0.215</td><td>0.099</td><td>1</td><td>0.753</td></tr><tr><td>Nutrient (N)</td><td>43.8</td><td>1</td><td>0.205</td><td>0.652</td><td>0.049</td><td>1</td><td>0.497</td><td>0.483</td><td>0.328</td><td>1</td><td>0.567</td></tr><tr><td>Carbon enrichment (C)</td><td>20592.5</td><td>1</td><td>96.533</td><td><strong>&lt;0.001</strong></td><td>0.695</td><td>1</td><td>7.011</td><td><strong>&lt;0.01</strong></td><td>3.405</td><td>1</td><td>0.065</td></tr><tr><td>Light (L)</td><td>2873.5</td><td>1</td><td>13.470</td><td><strong>&lt;0.001</strong></td><td>0.082</td><td>1</td><td>0.823</td><td>0.365</td><td>1.188</td><td>1</td><td>0.276</td></tr><tr><td>T:N</td><td>480.1</td><td>1</td><td>2.251</td><td>0.138</td><td>0.117</td><td>1</td><td>1.177</td><td>0.281</td><td>0.461</td><td>1</td><td>0.497</td></tr><tr><td>T:C</td><td>552.8</td><td>1</td><td>2.592</td><td>0.111</td><td>0.099</td><td>1</td><td>1.001</td><td>0.320</td><td>5.498</td><td>1</td><td><strong>0.019</strong></td></tr><tr><td>N:C</td><td>286.6</td><td>1</td><td>1.343</td><td>0.250</td><td>0.006</td><td>1</td><td>0.058</td><td>0.811</td><td>0.417</td><td>1</td><td>0.518</td></tr><tr><td>T:L</td><td>22.6</td><td>1</td><td>0.106</td><td>0.746</td><td>0.062</td><td>1</td><td>0.623</td><td>0.432</td><td>0.383</td><td>1</td><td>0.536</td></tr><tr><td>N:L</td><td>26.2</td><td>1</td><td>0.123</td><td>0.727</td><td>0.000</td><td>1</td><td>0.002</td><td>0.963</td><td>11.533</td><td>1</td><td><strong>&lt;0.001</strong></td></tr><tr><td>C:L</td><td>0.9</td><td>1</td><td>0.004</td><td>0.949</td><td>0.168</td><td>1</td><td>1.689</td><td>0.198</td><td>0.876</td><td>1</td><td>0.349</td></tr><tr><td>T:N:C</td><td>758.5</td><td>1</td><td>3.556</td><td>0.063</td><td>0.016</td><td>1</td><td>0.165</td><td>0.686</td><td>0.843</td><td>1</td><td>0.359</td></tr><tr><td>T:N:L</td><td>17.2</td><td>1</td><td>0.081</td><td>0.596</td><td>0.015</td><td>1</td><td>0.152</td><td>0.698</td><td>1.480</td><td>1</td><td>0.224</td></tr><tr><td>T:C:L</td><td>156.0</td><td>1</td><td>0.731</td><td>0.794</td><td>0.015</td><td>1</td><td>0.155</td><td>0.695</td><td>0.116</td><td>1</td><td>0.734</td></tr><tr><td>N:C:L</td><td>224.1</td><td>1</td><td>1.051</td><td>0.612</td><td>0.000</td><td>1</td><td>0.004</td><td>0.948</td><td>2.599</td><td>1</td><td>0.107</td></tr><tr><td>T:N:C:L</td><td>189.8</td><td>1</td><td>0.890</td><td>0.348</td><td>0.000</td><td>1</td><td>0.004</td><td>0.952</td><td>0.561</td><td>1</td><td>0.454</td></tr><tr><td>Residuals</td><td>16852.5</td><td>79</td><td>NA</td><td>NA</td><td>7.835</td><td>79</td><td>NA</td><td>NA</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**FIGURE 4** Predictive power (adjusted R²) of different models for resistance, resilience, and return time. Median and 95% confidence intervals are shown. Circles show data for the whole experiment, triangles show data when the carbon enrichment treatment is excluded</figcaption></figure>

longer-term, but not for shorter-term predictions. This result aligns with Christensen and collaborators' experiment (2006) in which they found that interactions between three environmental changes (temperature, drought, and acidification) were stronger and synergistic at the end of their experiment due to stimulated total zooplankton biomass. Future studies should examine if these results hold for other ecosystem variables and for population dynamics and community structure.

Carbon enrichment had the greatest effect on dissolved oxygen dynamics, reducing resistance, increasing return time, while accelerating recovery from perturbations. The short-term negative effect on DO was caused by increase in bacterial per capita respiration and abundance (Figure 5), with little or no change in photosynthesis, which is often observed in natural systems (Amon & Benner, 1996; Findlay, Sinsabaugh, Sobczak, & Hoostal, 2003). Indeed, half of the variance in DO is explained by the total bacteria density, reflecting the importance of bacterial abundance and respiration for dissolved oxygen concentration (Figure 5).

It is interesting to note that this lack of resistance results from the ability of the biological community to quickly respond to the increased carbon available; such quick responses may be desirable in some situations (e.g. population recovery from low abundance) such that lack of resistance may sometime be a desirable property. Depletion of available carbon and subsequent reduction in bacteria (attributable to ciliate and rotifer consumption) could cause oxygen concentration to return to preperturbation levels. Reduced light availability had the same directional effect as carbon enrichment, although smaller in magnitude, and was likely attributable to a different underlying process—reduced light availability may have reduced photosynthesis of existing algae and also reduced algal growth and thus slowed oxygen production (Brennan & Collins, 2015). It should be noted that the time scale of the response (e.g. maximum effect within 3 days of the disturbance) might be different for other response variables; the response of longer-lived organisms than bacteria would likely take longer. Understanding how disturbances

<figure><img src="image_11.png" /><figcaption>**FIGURE 5** Relationship between dissolved oxygen measurements and total bacteria densities for the 3 days following the disturbance(s). Symbols represent microcosms with (circle) or without (square) the carbon enrichment treatment. The dark symbols highlight the bacterial abundances estimated the day of maximum amount of change (i.e. resistance)</figcaption></figure>

affected community composition and structure may be important to pinpoint the mechanisms underlying the observed responses. Our study is currently limited by lacking the population dynamic data for all species. As soon as these data get available, we will use a set of more mechanistic models to understand the dynamics of DO and its different predictabilities.

We also found that increasing the number of perturbations decreased resistance and increased recovery time, but did not affect resilience. This effect may be explained by the greater chance that the dominant disturbance would be present when the number of perturbations increases (Brennan & Collins, 2015). Interestingly, for resistance, we also observed the detrimental effect of the number of perturbations in the absence of the dominant disturbance
(Figure 3d). While two environmental changes (temperature and nutrient) did not have significant effects (Table 2), both disturbances had an effect in the presence of at least one other disturbance.

Previous conceptual work (Brook et al., 2008) and modelling studies (Sala et al., 2000) highlighted how the effects of environmental changes may be stronger than expected due to positive interactions and synergies between global change disturbances. We found limited evidence for positive interactions in our controlled experiment manipulating four common environmental changes in aquatic systems. Instead, the large effect of carbon enrichment relative to other disturbances resulted in the dominance model explaining the data best. This result conflicts with the notion of the widespread occurrences of positive interactions (Brook et al., 2008), but is in line with recent meta-analyses (Darling & Côté, 2008; Jackson et al., 2016) and experiments (Brennan & Collins, 2015). Darling and colleagues studied the presence of interactions on animal mortality in 112 factorial experiments and found truly synergistic effects in only one-third of the reviewed studies, which was consistent across different disturbance types, as well as organisms and life stages considered. Moreover, Brennan and Collins (2015) looked at the growth response of a common freshwater algae under up to eight different types of environmental changes and found that the dominance model explained the data better than any additive or multiplicative model.

Findings can be influenced by experimental design choices. We had only two levels of each treatment (e.g. lower and higher temperature), as a result of choosing a relatively large number of environmental changes. This limited the type of interaction the experiment could reveal. It was unable to detect nonlinear effects, or how such nonlinearity could be affected by other disturbances. That is, the experimental design could not evaluate if interactions among disturbances were state dependent. To do this would have required continuous manipulations of multiple disturbances to construct a disturbance-effect surface. Examining effects of continuous variation in multiple disturbances should be a priority for future research. Our findings may have also been influenced by our choice of treatment levels; for example, if we had chosen a much smaller carbon enrichment treatment, we may have found less support for the dominance hypothesis, although when we excluded this treatment from our analyses, the relative importance of the models, and their scale dependence, changed little.

To date, very few ecological studies of multiple environmental changes have attempted to predict responses across time (Petchey et al., 2015). We show that the immediate response (resistance) was very well predicted with few assumptions and data (dominance hypothesis). The additive hypothesis is commonly used to predict environmental changes (Crain et al., 2008), but here we showed that the dominance hypothesis, which estimated the same number of parameters and required a smaller experimental design, predicted the immediate effect of environmental changes very well. In contrast to resistance, resilience and return time were less predictable due to increasing variability in return time among replicates (Figure 3c) and small or no effect on the resilience (Table 2). The predictability of return time was somewhat improved by

incorporation of interaction terms. Interactions could significantly affect species, but due to species cotolerance (Vinebrooke et al., 2004) or functional redundancy (Fonseca & Ganade, 2001), ecosystem functioning may not be subjected to interactions between multiple environmental changes. Low predictability of resilience in response to the four environmental changes may have been caused by the process underlying recovery. Recovery likely resulted from arrested bacterial growth and consumption of bacteria by ciliates and rotifers. If we assume that none of these four disturbances increased the predation rate, there would have been no effect on recovery rate, although predation should at least have been higher in the increased temperature treatment level, and thereby increasing recovery rate (Pellan, Médoc, Renault, Spataro, & Piscart, 2016).

How the findings of any individual experiment performed with a specific community at a particular spatial scale apply at larger spatial scales and for different communities is an open question that will require considerable ingenuity to address. Gradual accumulation of individual experiments eventually provides opportunities for meta-analyses of such issues, but such accumulation can take a long time and is usually not part of a strategic/directed research effort. A preferable option is for multiple laboratories to coordinate to perform a carefully planned collection of individual experiments, which can then be analyses in combination. A single all-encompassing experiment would manipulate multiple environmental drivers (as we did) and include manipulations of spatial scale and community complexity. As mentioned above, such an experiment would also involve gradients (rather than discrete levels) of environmental disturbance. This experiment would require unusually large amounts of resources (time, space funding, and personnel), although is not impossible to envisage. Finally, rigorous combining of findings from experiments and observational studies is a promising approach, and may be facilitated using process-based models and appropriate statistical methods of parameter inference (e.g. Bayesian methods).

The consequences of global environmental change on ecosystem stability are difficult to foresee, despite the urgent need for accurate predictions and recommendations to policy makers. Positive interactions have the potential to hamper such predictions, however, they may be less widespread than suspected. Our results hence support the statement of Darling and Côté (2008) that the "prevailing ecological paradigm of synergies are rampant" may be overstated. Instead, we documented that the most parsimonious model for a microbial aquatic experiment showed scale dependence. Understanding what can be predicted and what cannot, and how this depends on temporal scale, is a challenge for future studies to provide accurate tools for ecosystem management.

## ACKNOWLEDGEMENTS

We thank Y. Choffat for his help with the laboratory, C. Ward for her contribution and correcting the English, and the two anonymous reviewers for their comments that improved this study. This work was financially supported by the University Research Priority
Program: Global Change Biodiversity (to AG and OLP), the Swiss National Science Foundation grants (31003A_159498 to OLP and FP), and from the Marie Heim-Vögtlin Foundation (PMPDP3_151361 and PMPDP3_161858 to ML).

REFERENCES

Alternatt, F., Fronhofer, E. A., Garnier, A., Giometto, A., Hammes, F., Klecka, J., ... Petchey, O. L. (2015). Big answers from small worlds: A user's guide for protist microcosms as a model system in ecology and evolution (ed Spencer M). *Methods in Ecology and Evolution*, **6**, 218–231.

Amon, R. M. W., & Benner, R. (1996). Bacterial utilization of different size classes of dissolved organic matter. *Limnology and Oceanography*, **41**, 41–51.

Anneville, O., Souissi, S., Ibanez, F., Ginot, V., Druart, J. C., & Angeli, N. (2002). Temporal mapping of phytoplankton assemblages in Lake Geneva: Annual and interannual changes in their patterns of succession. *Limnology and Oceanography*, **47**, 1355–1366.

Arnoldi, J.-F., Loreau, M., & Haegeman, B. (2016). Resilience, reactivity and variability: A mathematical comparison of ecological stability measures. *Journal of Theoretical Biology*, **389**, 47–59.

Bates, D., & Mächler, M. (2015). Fitting linear mixed-effects models using lme4. *Journal of Statistical Software*, **67**, 1–48.

Brennan, G., & Collins, S. (2015). Growth responses of a green alga to multiple environmental drivers. *Nature Climate Change*, **5**, 892–897.

Brook, B., Sodhi, N., & Bradshaw, C. (2008). Synergies among extinction drivers under global change. *Trends in Ecology & Evolution*, **23**, 453–460.

Brown, J. R., Blankinship, J. C., Niboyet, A., van Greonigen, K. J., Dijkstra, P., Le Roux, X., ... Hungate, B. A. (2012). Effects of multiple global change treatments on soil N2O fluxes. *Biogeochemistry*, **109**, 85–100.

Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M., & West, G. B. (2004). Toward a metabolic theory of ecology. *Ecology*, **85**, 1771–1789.

Carpenter, S. R., Caraco, N. F., Correll, D. L., Howarth, R. W., Sharpley, A. N., & Smith, V. H. (1998). Nonpoint pollution of surface waters with phosphorus and nitrogen. *Ecological Applications*, **8**, 559.

Carpenter, S. R., Fisher, S. G., Grimm, N. B., & Kitchell, J. F. (1992). Global change and freshwater ecosystems. *Annual Review of Ecology and Systematics*, **23**, 119–139.

Chapin III, F. S., Zavaleta, E. S., Eviner, V. T., Naylor, R. L., Vitousek, P. M., Reynolds, H. L., ... Diaz, S. (2000). Consequences of changing biodiversity. *Nature*, **405**, 234–242.

Christensen, M. R., Graham, M. D., Vinebrooke, R. D., Findlay, D. L., Paterson, M. J., & Turner, M. A. (2006). Multiple anthropogenic stressors cause ecological surprises in boreal lakes. *Global Change Biology*, **12**, 2316–2322.

Côté, I. M., Darling, E. S., & Brown, C. J. (2016). Interactions among ecosystem stressors and their importance in conservation. *Proceedings of the Royal Society B: Biological Sciences*, **283**, 20152592.

Crain, C. M., Kroeker, K., & Halpern, B. S. (2008). Interactive and cumulative effects of multiple human stressors in marine systems. *Ecology Letters*, **11**, 1304–1315.

Daily, G. C., Söderqvist, T., Aniyar, S., Arrow, K., Dasgupta, P., Ehrlich, P. R., ... Walker, B. (2000). The value of nature and the nature of value. *Science*, **289**, 395–396.

Darling, E. S., & Côté, I. M. (2008). Quantifying the evidence for ecological synergies. *Ecology Letters*, **11**, 1278–1286.

Denny, M., & Benedetti-Cecchi, L. (2012). Scaling up in ecology: Mechanical approaches. *Annual Review of Ecology, Evolution, and Systematics*, **43**, 1–22.

Diamond, J. M., Ashmole, N. P., & Purves, P. E. (1989). The present, past and future of human-caused extinctions [and discussion]. *Philosophical Transactions of the Royal Society B: Biological Sciences*, **325**, 469–477.

Diaz, R., & Breitburg, D. (2009). Chapter 1: The hypoxic environment. In: *Fish physiology*, vol. 27 (eds Richards J, Farrell A, Braumer C), pp. 1–23. Academic Press, New York, NY.

Diaz, R. J., & Rosenberg, R. (2008). Spreading dead zones and consequences for marine ecosystems. *Science*, **321**, 926–929.

Doak, D. F., Estes, J. A., Halpern, B. S., Jacob, U., Lindberg, D. R., Lovvorn, J., ... Novak, M. (2008). Understanding and predicting ecological dynamics: Are major surprises inevitable. *Ecology*, **89**, 952–961.

Donohue, I., Hillebrand, H., Montoya, J. M., Petchey, O. L., Pimm, S. L., Fowler, M. S., ... Yang, Q. (2016). Navigating the complexity of ecological stability. *Ecology Letters*, **19**, 1172–1185.

Downing, J. A., & McCauley, E. (1992). The nitrogen: Phosphorus relationship in lakes. *Limnology and Oceanography*, **37**, 936–945.

Doyle, S. A., Saros, J. E., & Williamson, C. E. (2005). Interactive effects of temperature and nutrient limitation on the response of alpine phytoplankton growth to ultraviolet radiation. *Limnology and Oceanography*, **50**, 1362–1367.

Findlay, S. E. G., Sinsabaugh, R. L., Sobczak, W. V., & Hoostal, M. (2003). Metabolic and structural response of hyporheic microbial communities to variations in supply of dissolved organic matter. *Limnology and Oceanography*, **48**, 1608–1617.

Folt, C. L., Chen, C. Y., Moore, M. V., & Burnaford, J. (1999). Synergism and antagonism among multiple stressors. *Limnology and Oceanography*, **44**, 864–877.

Fonseca, C. R., & Ganade, G. (2001). Species functional redundancy, random extinctions and the stability of ecosystems. *Journal of Ecology*, **89**, 118–125.

Fox, J., & Weisberg, S. (2011). *An R companion to applied regression*, 2nd ed. Thousand Oaks (CA): Sage.

Garcia, H. E., & Gordon, L. I. (1992). Oxygen solubility in seawater: Better fitting equations. *Limnology and oceanography*, **37**, 1307–1312.

Griffen, B. D., & Drake, J. M. (2008). A review of extinction in experimental populations. *Journal of Animal Ecology*, **77**, 1274–1287.

Griffiths, J. I., Warren, P. H., & Childs, D. Z. (2015). Multiple environmental changes interact to modify species dynamics and invasion rates. *Oikos*, **124**, 458–468.

Hanson, P. C., Carpenter, S. R., Armstrong, D. E., Stanley, E. H., & Kratz, T. K. (2006). Lake dissolved inorganic carbon and dissolved oxygen: Changing drivers from days to decades. *Ecological Monographs*, **76**, 343–363.

Heugens, E. H., Hendriks, A. J., Dekker, T., van Straalen, N. M., & Admiraal, W. (2001). A review of the effects of multiple stressors on aquatic organisms and analysis of uncertainty factors for use in risk assessment. *Critical reviews in toxicology*, **31**, 247–284.

Holling, C. S. (1996). Surprise for Science, Resilience for Ecosystems, and Incentives for People. *Ecological Applications*, **6**, 733–735.

Holmstrup, M., Bindesbøl, A.-M., Oostingh, G. J., Duschl, A., Scheil, V., Köhler, H.-R., ... Spurgeon, D. J. (2010). Interactions between effects of environmental chemicals and natural stressors: A review. *Science of The Total Environment*, **408**, 3746–3762.

IPCC (2007). *Climate change 2007: The physical science basis: Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change*. Cambridge; New York: Cambridge University Press, (p. 996)

Isbell, F., Craven, D., Connolly, J., Loreau, M., Schmid, B., Beierkuhnlein, C., ... Eisenhauer, N. (2015). Biodiversity increases the resistance of ecosystem productivity to climate extremes. *Nature*, **526**, 574–577.

Jackson, M. C., Loewen, C. J. G., Vinebrooke, R. D., & Chimimba, C. T. (2016). Net effects of multiple stressors in freshwater ecosystems: A meta-analysis. *Global Change Biology*, **22**, 180–189.

Jankowski, T., Livingstone, D. M., Bührer, H., Forster, R., & Niederhauser, P. (2006). Consequences of the 2003 European heat wave for lake temperature profiles, thermal stability, and hypolimnetic oxygen depletion: Implications for a warmer world. *Limnology and Oceanography*, **51**, 815–819.
Jenny, J.-P., Francus, P., Normandeau, A., Lapointe, F., Perga, M.-E., Ojala, A., ... Zolitschka, B. (2016). Global spread of hypoxia in freshwater ecosystems during the last three centuries is caused by rising local human pressure. *Global Change Biology*, 22, 1481–1489.

Kaunzinger, C. M. K., & Morin, P. J. (1998). Productivity controls food-chain properties in microbial communities. *Nature*, 395, 495–497.

Kirk, J. T. O. (1983). *Light and photosynthesis in aquatic ecosystems*. Cambridge University Press, Cambridge.

Kolzau, S., Wiedner, C., Rücker, J., Köhler, J., Köhler, A., & Dolman, A. M. (2014). Seasonal patterns of nitrogen and phosphorus limitation in four german lakes and the predictability of limitation status from ambient nutrient concentrations. *PLoS ONE*, 9, e96065.

Kratina, P., Greig, H. S., Thompson, P. L., Carvalho-Pereira, T. S. A., & Shurin, J. B. (2012). Warming modifies trophic cascades and eutrophication in experimental freshwater communities. *Ecology*, 93, 1421–1430.

Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2016). ImerTest: tests in linear mixed effects models.

Lawler, S. P., & Morin, P. J. (1993). Food web architecture and population dynamics in laboratory microcosms of protists. *The American naturalist*, 141, 675–686.

Llames, M. E., Lagomarsino, L., Diovisalvi, N., Fermani, P., Torremoreli, A. M., Perez, G., ... Zagarese, H. E. (2009). The effects of light availability in shallow, turbid waters: A mesocosm study. *Journal of Plankton Research*, 31, 1517–1529.

Mantyka-pringle, C. S., Martin, T. G., & Rhodes, J. R. (2012). Interactions between climate and habitat loss effects on biodiversity: A systematic review and meta-analysis. *Global Change Biology*, 18, 1239–1252.

Moore, M. V., Folt, C. L., & Stemberger, R. S. (1996). Consequences of elevated temperatures for zooplankton assemblages in temperate lakes. *Archiv für Hydrobiologie*, 135, 289–319.

Moran, R., Harvey, I., Moss, B., Feuchtmayr, H., Hatton, K., Heyes, T., & Atkinson, D. (2010). Influence of simulated climate change and eutrophication on three-spined stickleback populations: A large scale mesocosm experiment. *Freshwater Biology*, 55, 315–325.

Neubert, M. G., Caswell, H., & Solow, A. R. (2009). Detecting reactivity. *Ecology*, 90, 2683–2688.

Pellan, L., Médoc, V., Renault, D., Spataro, T., & Piscart, C. (2016). Feeding choice and predation pressure of two invasive gammarids, Gammarus tigrinus and Dikerogammarus villosus, under increasing temperature. *Hydrobiologia*, 781: 43. https://doi.org/10.1007/s10750-015-2312-3

Petchey, O., McPhearson, P., Casey, T., & Morin, P. (1999). Environmental warming alters food-web structure and ecosystem function. *Nature*, 402, 69–72.

Petchey, O. L., Pontarp, M., Massie, T. M., Kéfi, S., Ozgul, A., Weilenmann, M., ... Pearse, I. S. (2015). The ecological forecast horizon, and examples of its uses and determinants (ed Vasseur D). *Ecology Letters*, 18, 597–611.

Piehler, M. F., Twomey, L. J., Hall, N. S., & Paerl, H. W. (2004). Impacts of inorganic nutrient enrichment on phytoplankton community structure and function in Pamlico Sound, NC, USA. *Estuarine, Coastal and Shelf Science*, 61, 197–209.

Piggott, J. J., Townsend, C. R., & Matthaei, C. D. (2015). Reconceptualizing synergism and antagonism among multiple stressors. *Ecology and Evolution*, 5, 1538–1547.

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. Vienna, Austria: R Foundation for Statistical Computing, URL http://www.R-project.org/.

Regnier, P., Friedlingstein, P., Ciais, P., Mackenzie, F. T., Gruber, N., Janssens, I. A., ... Thullner, M. (2013). Anthropogenic perturbation of the carbon fluxes from land to ocean. *Nature Geoscience*, 6, 597–607.

Sala, O. E., Chapin III, F. S., Armesto, J. J., Berlow, E., Bloomfield, J., Dirzo, R., ... Wall, D. H. (2000). Global biodiversity scenarios for the year 2100. *Science*, 287, 1770–1774.

Savage, V. M., Gillooly, J. F., Brown, J. H., West, G. B., & Charnov, E. L. (2004). Effects of body size and temperature on population growth. *The American Naturalist*, 163, 429–441.

Smith, V. H., Tilman, G. D., & Nekola, J. C. (1999). Eutrophication: Impacts of excess nutrient inputs on freshwater, marine, and terrestrial ecosystems. *Environmental Pollution*, 100, 179–196.

Stanley, E. H., Powers, S. M., Lottig, N. R., Buffam, I., & Crawford, J. T. (2012). Contemporary changes in dissolved organic carbon (DOC) in human-dominated rivers: Is there a role for DOC management? *Freshwater Biology*, 57, 26–42.

Stendera, S., Adrian, R., Bonada, N., Cañedo-Argüelles, M., Hugueny, B., Januschke, K., ... Hering, D. (2012). Drivers and stressors of freshwater biodiversity patterns across different ecosystems and scales: A review. *Hydrobiologia*, 696, 1–28.

Stets, E. G., & Cotner, J. B. (2008). The influence of dissolved organic carbon on bacterial phosphorus uptake and bacteria-phytoplankton dynamics in two Minnesota lakes. *Limnology and Oceanography*, 53, 137–147.

Stomp, M., Huisman, J., de Jongh, F., Veraart, A. J., Gerla, D., Rijkeboer, M., ... Stal, L. J. (2004). Adaptive divergence in pigment composition promotes phytoplankton biodiversity. *Nature*, 432, 104–107.

Therneau, T. M. (2015). A package for survival analysis in S.

Vinebrooke, R. L., Cottingham, K., Norberg, M. S., Jon, I., Dodson, S. C., Maberly, S., & Sommer, U. (2004). Impacts of multiple stressors on biodiversity and ecosystem functioning: The role of species co-tolerance. *Oikos*, 104, 451–457.

Walker, W. (1979). Use of hypolimnetic oxygen depletion rate as a trophic state index for lakes. *Water Resources Research*, 15, 1463–1470.

Webster, J. R., & Benfield, E. F. (1986). Vascular plant breakdown in freshwater ecosystems. *Annual Review of Ecology and Systematics*, 17, 567–594.

Wetzel, R. G. (1983). *Limnology* (2nd ed.). Philadelphia, PA: Saunders College Publishing, 876pp.

Wetzel, A., & Uchman, A. (2001). Sequential colonization of muddy turbidites in the Eocene Beloveža Formation, Carpathians, Poland. *Palaeogeography, Palaeoclimatology, Palaeoecology*, 168, 171–186.

Williams, C. J., Frost, P. C., Morales-Williams, A. M., Larson, J. H., Richardson, W. B., Chiandet, A. S., & Xenopoulos, M. A. (2016). Human activities cause distinct dissolved organic matter composition across freshwater ecosystems. *Global Change Biology*, 22, 613–626.

Wu, R. S. S. (2002). Hypoxia: From molecular responses to ecosystem responses. *Marine Pollution Bulletin*, 45, 35–45.

Yankova, Y., Villiger, J., Pernthaler, J., Schanz, F., & Posch, T. (2016). Promotion, deepening and warming of the metalinmion change habitat conditions of the harmful filamentous cyanobacterium Planktothrix rubescens in a prealpine lake. *Hydrobiologia*, 776, 125–138.

Yvon-Durocher, G., Jones, J. I., Trimmer, M., Woodward, G., & Montoya, J. M. (2010). Warming alters the metabolic balance of ecosystems. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 365, 2117–2126.

## SUPPORTING INFORMATION

Additional Supporting Information may be found online in the supporting information tab for this article.

**How to cite this article:** Garnier A, Pennekamp F, Lemoine M, Petchey OL. Temporal-scale-dependent interactions between multiple environmental disturbances in microcosm ecosystems. *Glob Change Biol.* 2017;00:1–12. https://doi.org/10.1111/gcb.13786