1. Ecosystems are connected by flows of nutrients and organisms. Changes to connectivity and nutrient enrichment may destabilise ecosystem dynamics far from the nutrient source.
2. We used gradostats to examine the effects of trophic connectivity (movement of consumers and producers) versus nutrient-only connectivity on the dynamics of Daphnia pulex (consumers) and algae (resources) in two metaecosystem configurations (linear vs. dendritic).
3. We found that Daphnia peak population size and instability (coefficient of variation; CV) increased as distance from the nutrient input increased, but these effects were lower in metaecosystems connected by all trophic levels compared to nutrient-only connected systems and/or in dendritic compared to linear systems.
4. We examined the effects of trophic connectivity (i.e. both trophic levels move rather than one or the other) using a generic model to qualitatively assess whether the expectations align with the ecosystem dynamics we observed.
5. Analysis of our model shows that increased Daphnia population sizes and fluctuations in consumer-resource dynamics are expected with nutrient connectivity, with this pattern being more pronounced in linear rather than dendritic systems.
6. These results confirm that connectivity may propagate and even amplify instability over a metaecosystem to communities distant from the source disturbance, and suggest a direction for future experiments, that recreate conditions closer to those found in natural systems.

Our gradostat flasks contained simple communities of the water flea Daphnia pulex consuming a mix of three algal species (Pseudokirchneriella subcapitata, Scenedesmus quadricauda, Ankistrodesmus falcatus). This experiment employed a 2x2x2 factorial design to test the importance of ecosystem trophic connectivity (a treatment considering movement of medium only vs. movement of media, phytoplankton and Daphnia between flasks) and metaecosystem configuration (linear or dendritic) on the stability of Daphnia populations and algal communities with two levels of enriched medium input (regular and phosphorus-enriched). Four replicates of this whole design were established, for a total of 32 metaecosystems, run in 9 blocks due to time and space constraints. Each metaecosystem consisted of four “nodes” of 500 mL Erlenmeyer flasks with a foam stopper to allow for gas exchange (128 flasks total), seeded initially with 100 mL algal mix (total average algal density of 2.22 x106 +/- 1.3x104 cells/mL) to which 50 adult Daphnia with eggs (which produce broods of about 15 individuals each week in good conditions (Schwartz 1984) were added before topping off the flask to 500 mL with FLAMES media (Celis-Salgado et al. 2008).

Configuration was controlled by unidirectionally connecting flasks in either a linear configuration (in →1→2→3→4→out) or a dendritic configuration (in →1, in→2, 1→3, 2→3, 3→4→out). We chose this as the simplest possible design in which a linear network could be compared to a branched network, with four nodes being the smallest possible number of nodes to create a dendritic configuration, and the two nodes branching into a third, similar to headwater in a river. Flasks were then connected by Tygon tubing and from an inflow reservoir of FLAMES medium (10 μgP/L) or enriched P (70 μgP /L) medium which was pumped through the array of flasks using peristaltic pumps (Watson-Marlow 503S/RL and Rainin Dynamix RP-1). Pumps were set on automatic timers to run for one hour each day at a speed adjusted to move a specific volume of media over that hour. The dilution rate was 10% of the total volume per for all flasks in the linear configurations and the “hub” (3) and “terminal” (4) nodes of the dendritic configurations (50 mL), and 5% per day (25 mL) for the “upstream” nodes in the dendritic configurations (Figure 1).

We also controlled functional connectivity, contrasting metaecosystem dynamics when only nutrients moved versus the case when nutrients, resources and consumers moved. To block the flow of organisms in the nutrient-only connectivity treatment, outflow tubing was placed inside an 80-µm nylon mesh held in place with the stopper. Due to colony formation of the phytoplankton and clogging of the mesh, this proved to be an effective retention mechanism also for the algal resources, thus we believe flow of algae was significantly reduced in these treatments compared the trophic connectivity treatments. Though it is possible a small portion of single cells were able to pass through, Scenedesmus is known to form four-cell colonies in the presence of consumers (which we also observed in our algal counts), which are too large to pass through the mesh. As D. pulex were unable to fit through the tubing or survive moving through the peristaltic pumps, in the trophic connectivity treatment, D. pulex were manually moved using a 2mL transfer pipette at a rate of 10% of the population per day (20% were moved after each sampling count as sampling was only done every two days) in all linear nodes and the hub and terminal dendritic nodes, and 5% per day (10% moved after sampling) in the upstream dendritic nodes, in the same downstream direction as media. This type of passive movement at the flow rate of the system would be typical of planktonic animals in rivers that cannot swim upstream. Inflow stock solutions were prepared using FLAMES media (10 μgP/L). 

Finally, we modified our inflow reservoirs to contain either additionally P-enriched (high P) or regular (low P) FLAMES media. To increase P in the additionally nutrient-enriched treatment without changing pH, 132 μg/L of H₂KPO₄ and 168 μg/L H₂KPO₄ were added to our increased P treatment inflow stock solution. For the less-phosphorus treatment, no additional phosphorus was added, but 218 μg/L KCl were added to control for the K added to the high-P medium. See Figure S1 for a photograph of the experimental setup.

Experimental Sampling

The gradostats were sampled every other day for 30 days. In each node, the concentration of each algal species was measured using a haemocytometer. To estimate Daphnia population size, a 2mL plastic transfer pipette was used to gently agitate, and then sample each node. The number of individuals and two age classes (adult or juvenile) in the pipette were determined and then replaced to the experimental flask. This process was repeated five times, and the average D. pulex count of the five samples was used to estimate Daphnia density/2mL (total number estimated per flask = sampled count average *250). A pilot experiment testing this method proved it had an average error of 17.41 %, equating to 2.5 Daphnia more or less than the expected count at known densities; there is no reason to believe this error was systematic in one direction or the other, or to be systematically biased among our treatments. On Day 30 of the experiment, 40mL samples were taken from each flask to be analysed for total phosphorus concentration (TP). Phosphorus samples were analysed using a standard protocol (Wetzel and Likens 2013) at the GRIL-Université du Québec à Montréal analytical laboratory.

Statistical Analysis

To quantify the instability of Daphnia populations in experimental gradostats, we determined the peak total Dapnhia population size (as estimated by our density samples) and the coefficient of variation (CV) of Daphnia population size over the course of the experiment. These variables were calculated for each node within each gradostat, as well as in aggregate summed across all nodes for additive Daphnia metapopulation peak and CV. Similarly, population CV and peak density were calculated for each species of alga but we analyse here values based on total algal community density (sum of all species present), as Pseudorkirchinella and Ankistrodesmus were undetectable in most flasks for most of the experiment. Scenedesmus was mostly observed in 4-cell colonies, which is common in the presence of consumers, but we counted the total number of cells, not colonies.

All analyses of experimental gradostat data were conducted in R version 4 (Team 2020). Statistical tests of the hypothesis were two-sided and with a level of significance of α=0.05.

To determine whether metaecosystem connectivity, configuration and nutrient enrichment, as well as node position (1 upstream to 4 terminal), influenced node Daphnia population instability downstream of the nutrient enrichment source, we analysed the effects of these factors on mean Daphnia population and algal community peak values, on mean Daphnia population and algal community CV log-transformed (natural logarithm) values, and on mean final TP concentrations values, using linear mixed-effects models with the four factors as fixed effects. The mixed model included a random effect for ‘system’ which allowed us to account for a possible clustering in the response variables since the four nodes were connected as metaecosystems. For each of these models, pairwise interactions between factors were tested and terms for non-significant interactions were removed from the final models we report. Assumptions on the model errors (randomness, normality, and homoscedasticity) and the presence of possible influential observations or outliers were assessed with diagnostic plots of the model residuals. Robust standard errors (Huang and Li 2022) were used to adjust for heteroscedasticity.

We also measured Daphnia metapopulation and algal metacommunity instability at the scale of the entire metaecosystem. To determine whether metaecosystem connectivity, configuration and nutrient enrichment influenced Daphnia metapopulation and algal metacommunity instability, we analysed the effects of these factors on mean Daphnia metapopulation and algal metacommunity peak values, and on mean CV values, using linear mixed-effects models with the three factors as fixed effects, using the block in which a metaecosystem was run as a random factor. For each of these models, pairwise interactions between factors were tested and terms for non-significant interactions were removed from the final models we present. Assumptions on the model errors (randomness, normality, and homoscedasticity) and the presence of possible influential observations or outliers were assessed with diagnostic plots of the model residuals. Robust standard errors (Huang and Li 2022) were used to adjust for heteroscedasticity.