Materials and methods

Study design

In 20 month long experiment, we manipulated the presence and absence of two keystone species: Myriophyllum spicatum (Fig. 1B;hereafter Myriophyllum ) and Dreissena polymorpha (Fig. 1C; hereafter Dreissena ) in artificial pond ecosystems (15 000L). We used a fully factorial design with either both keystone species absent as a control (C), Myriophyllum macrophytes (M),Dreissena mussels (D) or Myriophyllum and Dreissenatogether (MD). Each factorial treatment combination of eutrophied ponds was replicated four times (4 x 4 = 16 ponds total). The ponds we used were made of fiberglass with a smooth surface (Fig. 1D), had a rounded shape with approximately four meter diameter and a shallow (0.5 m) and a deep (1.5 m) end. In the first period of perturbation, we progressively increased the pulse perturbation of inorganic nutrients to all ponds, and measured the effect of presence or absence of both keystone species on several ecosystem parameters in high frequency using automated multiparameter sondes. In the second period of perturbation, one year later, we perturbed all the ponds again with our highest pulse of nutrients, to test how the responses changed between the perturbation periods, and whether treatment effects were repeatable over two years.

Experimental procedure

The ponds were initially set up on May 6th 2016 by adding a 5 cm thick layer of gravel (2-4 mm) and filling them with a tap-water. Afterwards the ponds were inoculated with a natural phytoplankton population (20 L, 30 µm filtered from Lake Greifensee). The treatments were established on May 31st by distributing 100 shoots of Myriophyllum, each attached with a cable-tie to a small rock, among the shallow and deep ends of each pond designated to the M and MD treatment. Each pond that was designated to the D and MD treatment received 25 Dreissena , distributed among the shallow and deep end. We ensured prior to the distribution of plant shoots and mussels that their size distributions were similar across all ponds of the respective treatment. To ensure that all ponds started with a similar overall amount of total biomass, we added autoclaved mussels to the M ponds, autoclavedMyriophyllum shoots to the D ponds, and both autoclaved mussels and Myriophyllum shoots to the C ponds. In May 2017, we re-established our macrophyte treatment after the winter by adding the same amount of fresh and autoclaved Myriophyllum shoots to the respective ponds to ensure effective treatment contrasts.
The first nutrient perturbation regime began on August 12th, 2016. We progressively increased phosphate and nitrate additions of 10, 20, 30, 40 and 50 µg/L of P (with a double Redfield ratio, N:P = 30) over eight weeks until October 10th 2016, with two week intervals between additions. Our second nutrient perturbation regime was single addition of 50 µg/L of P on October 10th, 2017. Using multiparameter sondes (EXO2, Xylem), installed in each pond, we tracked the following four ecosystem parameters with high frequency (15 min intervals): chlorophyll-a fluorescence (hereafter chlorophyll) and phycocyanin fluorescence, DOM fluorescence (hereafter fDOM), and dissolved oxygen. Additionally, we used the dissolved oxygen data, as well as water temperature (also measured with the sondes), light and wind data, to calculate rates ecosystem metabolism (gross primary productivity, net primary productivity, and respiration; see below). Details on sonde calibration and maintenance can be found in the Supplement.
Over the first winter period (December 1st 2016 - February 28th 2017), we could not monitor ecosystem dynamics due to ice cover in the ponds. To maintain and recalibrate the sensors, we stopped measurement from March 1st to 23rd to maintain and recalibrate all sondes (see supplement for details). A second sonde maintenance period was implemented in the fall of 2017 (September 14th - October 3rd 2017). Following this structure, we consider three phases of the experiment: Phase 1 with the first five nutrient pulses (June - December 2016), Phase 2 without nutrient pulses (March - October 2017), and Phase 3 with the final nutrient pulse (October 2017 - February 2018).

Data treatment and analysis

Data treatment - We first performed an outlier analysis by excluding values higher than 3 times the median absolute deviation of all values in a sliding window (Leys et al. 2013) of one day window size (15 min interval = 96 data points). After aggregating four measurement points to one per hour (from 96 to 24 data points per day), we used sliding windows with a one-week window size (7 * 24 = 168 data points) to calculate mean and coefficient of variation (hereafter CV) of all ecosystem parameters. We chose a seven day window size to have robust estimates of the different metrics that would not be affected by diurnal variability. Moreover, we calculated autocorrelation (hereafter AC), tailedness of the generalized extreme value distribution (hereafter GEV) and skewness, which can be used to quantify the characteristics high frequency dynamics of disturbed ecosystems are (Batt et al.2013, 2017; Gsell et al. 2016). For example, as ecosystems are disturbed they tend to become more similar to their own past, resulting in an increase in AC (Ives 1995), whereas GEV, especially of biological variables, is expected to yield higher values (i.e. “fatter tails”) in response to perturbation (Katz et al.2005; Batt et al. 2017). We also focused on these metrics because they are important early warning indicators for critical transitions in shallow lake ecosystems (Carpenter et al. 2011; Gsell et al. 2016), which, however, have been rarely investigated in factorial manipulation of foundation species.
Effects of foundation species on ecosystem parameters - Using the data derived from the sliding windows, we tested for differences between treatments using the factorial design (n=4 per treatment level, aD = main effect of Dreissena, bM = main effect of Myriophyllum, C(DxM) = interactive effect):
\begin{equation} \mathbf{y\ =\ }\mathbf{a}_{\mathbf{D}}\mathbf{+\ }\mathbf{b}_{\mathbf{M}}\mathbf{+\ }\mathbf{C}_{\mathbf{(D\times M)}}\mathbf{+error}\nonumber \\ \end{equation}
We used one linear model with Type III sum of squares per hour (24 models per day) to test for differences between treatments in mean and CV, AC, GEV, and skewness of each measured parameter. We report P values from linear models for mean and CV directly in Figure 2 and 3, where points below the time series colour coded by treatment indicate a significant difference of the respective treatment from the control. Because there were no systematic differences between treatments for AC, GEV, and skewness, we report results for these metrics in Supplementary Figures S1-S3. For better visual inference in all the presented figures, we further aggregated data from the sliding windows from 24 to one data point per day. In addition, we calculated the predicted additive response of Myriophyllum and Dreissena for each data-point by subtracting the control from the summed single species treatments ((Dreissena + Myriophyllum ) - Control). The interaction between the presence of Myriophyllum and Dreissena was considered non-additive, when the confidence interval of the MD-treatment did not overlap with the predicted additive response.
Ecosystem metabolism - We calculated gross primary productivity, net primary productivity, and respiration (hereafter GPP, NEP and R, respectively) of each pond ecosystem using the equations in Staehr et al. (2010) on time series of dissolved oxygen and temperature collected by the sondes, as well as wind speed at 10 m from a nearby weather station operated by Meteo Swiss (Dübendorf, Giessen). Because the ponds were oversaturated with respect to dissolved oxygen, we included rates of change in dissolved oxygen in the formulas as the coefficient of a linear model of hourly averages of dissolved oxygen concentrations between 13:00 and 17:00 for the day and 1:00 and 5:00 for the night, where gas exchange dynamics in the ponds were considered to have equilibrated. This assumption was tested by visually inspecting the slopes for oxygen increase or decrease, which were found to be linear within these times and across all seasons. Using the metabolism data we calculated mean and CV of all three metabolism parameters by applying a sliding window with the size of 7 days. We then tested for differences between treatments with single species (M and D - main effect) and multiple species (MD - interactive effect) and control (C) using one linear model per day. We report the results from the linear models directly in Figure 4 and 5 as colour coded points that indicate significant difference in metabolic rates of M, D or MD from C. Furthermore, we calculated the predicted additive effect in the same fashion as for the other ecosystem parameters.