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.