An Extended Ecosystem Model for Understanding EE2 Indirect Effects on a Freshwater Food Web and its Ecosystem Function Resilience

Abstract

Freshwater species are highly impacted by human activities and the consequences on ecosystem functioning are still not well understood. In the literature, a multitrophic perspective appears to be key to advance future biodiversity and ecosystem functioning (BEF) research. This paper aims at studying indirect effects of the synthetic hormone 17α-ethinylestradiol (EE2) on a freshwater food web by creating BEF links, through the interpretation of seasonal cycles and multitrophic interactions. An ecosystem model previously developed using experimental data from a unique whole-ecosystem study on EE2 was extended with the addition of Chaoborus, an omnivorous insect. During the experimental study, a collapse of fathead minnow was measured after one year of exposure. The simulation results showed that EE2 indirect effects on other fishes (horizontal diversity) and lower trophic levels (vertical diversity) were connected to multitrophic interactions with a top-down cascade effect. The results also demonstrated that adding an omnivorous, mid-trophic level group such as Chaoborus enhances resilience. Conversely, missing such a species means that the actual resilience of an ecosystem and its functioning cannot be properly simulated. Thus, the extended ecosystem model offers a tool that can help better understand what is happening after environmental perturbations, such as with EE2.

1. Introduction

Pharmaceuticals and hormones are continuously entering the environment through release from wastewater treatment plants (WWTPs). In the last few decades, the occurrence of those emerging compounds has increasingly received attention due to their harm to the ecosystem. A review on worldwide environmental monitoring of water organic pollutants identified by European guidelines [1,2] demonstrated that pharmaceuticals represent the most studied class of Watch List compounds (42.9%), the natural estrogens E1 and E2 and the synthetic 17α-ethinylestradiol (EE2) being the second most studied (35.3%). Both the natural and synthetic hormones were previously identified as the main substances responsible for estrogenic activities in WWTPs, with concentrations as low as ng/L, sometimes exceeding their predicted no-effect concentrations (PNECs) for ecological toxicity [3]. Species responses to these environmental estrogens vary considerably, with fishes having much higher sensitivities than invertebrates [4,5]. Indirect effects on tolerant species, mediated by ecological mechanisms, may also appear in the environment but cannot be measured by single species laboratory-based toxicity tests; this is why population, community or ecosystem level studies are required [6].

EE2, one of the most potent synthetic estrogens, interferes in multiple ways with the activity of different physiological endpoints of aquatic organisms, including the endocrine system [7,8,9,10]. In a multi-year whole-ecosystem study performed at an experimental lake, [11] exposed well-defined fish and lower-trophic level populations to environmentally relevant concentrations of EE2. During this experimental study, a collapse of fathead minnow was measured after one year of exposure to EE2. A mix of reduced gamete production, increased gamete mortality and an increase of both adults and juvenile mortality was hypothesized to explain the collapse of both adult and juvenile fathead minnow due to EE2 addition in reference [12].

Since several small fishes declined in the experimental lake after EE2 additions, the authors sought evidence for a trophic cascade resulting from the associated reductions in planktivory and benthivory fish [13]. Trophic cascades (indirect effects mediated through consumer–resource interactions) are a well-studied type of indirect effect, and are generally considered in terms of ‘top–down’ (predator influence on lower trophic levels) or ‘bottom–up’ (nutrient/food/prey influence on higher trophic levels) mechanism [6]. Unlike the many studies that have shown direct effects of estrogens and their mimics on aquatic species, especially fishes, Kidd et al. [13] advance the understanding of how estrogen-induced changes in fish abundance can lead to indirect effects on freshwater food webs. This paper aims at taking further their interpretation of the direct and indirect effects of EE2 by modeling links between aquatic biodiversity and ecological functioning.

In terms of extinction risk, freshwater species are among the world’s most threatened [14]. Thus, one important question to be answered is: “To what extent does species loss from freshwaters affect ecosystem functioning and their ability to provide ecosystem goods and services for people?”. “Ecosystem function” (EF) is a general term that includes stocks of materials (e.g., carbon, water and mineral nutrients) and rates of processes involving fluxes of energy and matter between trophic levels and the environment. EFs are linked to biodiversity, and can thus be broadly defined as the biological, geochemical and physical processes that take place or occur within an ecosystem. Due to such complexity, the underlying role of biodiversity in ecosystem functioning, its relevance for ecosystem service provision in general, as well as the consequences of its decline, remain poorly understood [15]. A multitrophic perspective is key to advancing future biodiversity and ecosystem functioning (BEF) research and to address some of its most important remaining challenges [16].

Research on biodiversity and ecosystem functioning is one of the most prominent topics investigated by ecologists in the last three decades, with two new domains (Aquatic Food Webs and Agricultural Landscapes) arising in the last decade [17]. Besides, the “ecosystem service” has become the most prominent science policy term used in the last 8 years (2011–2018). Ecosystem-based management (EBM) is a collaborative management approach used with the intention to restore, enhance and protect the resilience of an ecosystem so as to sustain or improve ecosystem services (ES) and protect biodiversity, while considering nature and society [18]. In their review, Eisenhauer et al. [16] concluded that understanding why and how the strength of biodiversity effects varies with environmental conditions and at which spatial scales different mechanisms operate, will be key to operationalize BEF insights for ecosystem management, society and decision making.

Duffy et al. [19] pointed out that understanding how biodiversity affects functioning of complex ecosystems will benefit from integrating theory and experiments with simulations and network-based approaches. Indeed, integrated models could highlight priorities for the collection of new empirical data, identify gaps in our existing theories of how ecosystems work, help develop new concepts for how biodiversity composition and ecosystem function interact, and allow predicting BEF relations and its drivers at larger scales [20]. Besides, Queirós et al. [21] identified a need to translate observed mechanisms (short-term and single species) into conceptual models that include combinations of species and environmental gradients not yet observed (multiple stressors, long-term and interacting species).

According to Barnes et al. [22], integrating trophic complexity is key to understanding how biodiversity affects whole-ecosystem functioning. Therefore, this study takes up the challenge to improve an ecosystem model that was previously developed and successfully calibrated [12], with the addition of Chaoborus, an omnivorous insect that is an important food source for higher trophic levels and preys on lower trophic levels. This extended ecosystem model will be used to confirm and take further the hypotheses on EE2 indirect effects on a freshwater food web [13]. Besides, indirect EE2 effects will be compared with and without the presence of Chaoborus in the food web. Finally, how the presence of Chaoborus enhances the resilience of ecosystem functioning will be discussed.

Consultation with the ecologists observing the ecosystem was made to verify that major ecosystem mechanisms were correctly represented in the model. This was necessary as a first-order assessment of model reliability when extrapolated to conditions beyond the observed state [21]. The results presented in this paper will support the transition of BEF from a description of patterns to a predictive science [16], which can help better understand the link between biodiversity and ecosystem functioning in the context of EBM.

2. Materials and Methods

2.1. Step 1: Collecting Experimental Data

The experimental data used for the model were collected during a multi-year whole-lake study on Lake 260 in the Experimental Lakes Area (ELA) in northwestern Ontario, Canada (all sampling details can be found in [13]). This experimental lake is oligotrophic (high oxygen and low nutrient concentrations) and typical of boreal shield lakes. Lake 260 has an area of 34 ha and a maximum depth of 14.4 m. Six other experimental lakes in the same area were also studied as reference systems.

The study started with two years of baseline data, followed by three years of data collected during the addition of EE2, and continued for seven years after the addition was stopped to assess ecosystem stability and recovery after stressor removal. EE2 was added to the epilimnion for 20–21 weeks during lake stratification. Seasonal mean concentrations for the summer were between 4.8 and 6 ng/l. Lower concentrations were measured under the ice during winter.

Physico-chemical data were collected monthly during the open-water season, after lake stratification in the spring until shortly before turnover in the autumn. The epilimnion is defined as the surface layer of water with uniform temperature (ignoring any shallow temporary stratification phenomena), while the hypolimnion is defined as the bottom layer of water, also with uniform temperature. More information about the physico-chemical data used for the model, such as temperature, oxygen, light, organic matter and nutrients, can be found in [12].

A simplified food web was selected with the most relevant populations of plankton and fish naturally present in the experimental Lake 260 (Figure 1). Phytoplankton (Group 1: chlorophyte, dinoflagellates and Cyanophyta; Group 2: Chrysophyta and Cryptophyta and Group 3: diatoms) and zooplankton (cladocerans, copepods and rotifers) samples were collected biweekly. Fish abundance data were based on catch-and-release methods using trap nets (spring and autumn, all species) and short (30 min) evening gill net sets on spawning shoals for lake trout (autumn). Biomass of the minnow species (fathead minnow and pearl dace) was estimated as the product of abundance and mean size from minnow trap captures, standardized by lake area. Mark-and-recapture techniques were used to estimate the abundance of lake trout (autumn data) and white sucker (spring data). Biomass was estimated as the product of abundance estimates and mean size, standardized by the lake area.

In the extended version of this model presented here, Chaoborus was added to the studied ecosystem, a dipteran zooplankton predator and an important food item for many fish species (Figure 1). Chaoborus were sampled in the experimental lake every four weeks during the ice-free season at least 1 h after sunset at a station located over the deepest part of the lake using vertical hauls from the entire water column [13]. Instars were separated using measures of head-capsule length. The emerging adult insects were collected weekly with pyramidal emergence traps.

Discussions with biologists and ecologists were necessary to identify the seasonal cycles of the species selected for the model (Figure 2). The objective was to highlight the potential role of seasonal cycle alterations on the food web after EE2 addition, as an aspect of the relation between biodiversity and ecosystem functions.

2.2. Step 2: Building and Calibrating an Ecosystem Model

A dynamic ecosystem model was previously built and calibrated in an object-oriented framework using simplified AQUATOX equations and the software package WEST (DHI Water Environment Health, Hörsholm, Denmark) [12,23]. The model consists of a set of objects, each describing the growth of a model population in terms of its biomass concentration using differential equations including biological processes such as assimilation, photosynthesis, respiration, consumption or mortality, and additional processes such as migration, diffusion or loading. By connecting different objects and defining feeding relationships between them, a customized food web can be designed (Figure 3). In the extended ecosystem model, multitrophic interactions were complexified, with the addition of a new box for Chaoborus.

The life cycle of Chaoborus is split between aquatic ecosystems, where it spends its larvae phase (Instars 1–4), and terrestrial ecosystems, for its adult phase as an insect (Figure 4). The eggs are then released in the aquatic environment and some of them turn into instars 1–2, which is called “recruitment”. The larvae become bigger and are “promoted” to instars 3–4. After an initial phase in water, adult insects “emerge” to the terrestrial environment.

The mass balance used for Chaoborus was adapted from the one used for fish in reference [12], with the active respiration (when fish swim) and gamete loss (only a small fraction of gametes results in viable fish) removed. Below only the equations that were changed are presented.

Instars 1–2

$$\frac{dInsec{t}_{1–2}}{dt}=Cons-Def-Resp-Exc-Mort-Pred-Promo+Recruit$$

(1)

Instars 3–4

$$\frac{dInsec{t}_{3–4}}{dt}=Cons-Def-Resp-Exc-Mort-Pred+Promo-Emerg$$

(2)

where Insect_1–2 (g/m³) is the biomass of Instars 1–2, Insect_3–4 (g/m³) is the biomass of Instars 3–4, Cons (g/m³.d) the consumption of particulate organic matter (POM), phyto- and zooplankton, Def (g/m³.d) the defecation of unassimilated food, Resp (g/m³.d) the respiratory loss, Exc (g/m³.d) the excretion of dissolved organic matter, Mort (g/m³.d) the non-predatory mortality, Pred (g/m³.d) the consumption of Chaoborus by fish, Promo (g/m³.d) the promotion from Instars 1–2 to Instars 3–4, Recruit (g/m³.d) the recruitment from viable eggs to Instars 1–2 and Emerg (g/m³.d) the emergence from Instars 3–4 to adult insects. Initial values for these differential equations were obtained by calibration (see below).

Emergence of aquatic insects represents a loss for the system. The model assumes that emergence is determined by the net rate of growth, considered as the sum of consumption and the loss terms other than mortality.

$$Emerg\left(t\right)=KEm\times \left(Cons\left(t\right)-Def\left(t\right)-Resp\left(t\right)-Exc\left(t\right)\right)$$

(3)

where Emerg (g/m³.d) is the emergence from Instars 3–4 to adult insects, KEm (unitless) the fraction of growth that goes to emergence, Cons (g/m³.d) the consumption of POM, phyto- and zooplankton, Def (g/m³.d) the defecation of unassimilated food, Resp (g/m³.d) the respiratory loss and Exc (g/m³.d) the excretion of dissolved organic matter.

Often, there is synchrony in insect emergence and in the model, this is assumed to be cued to temperature, with additional forcing as twice the emergence that would ordinarily be computed. During recruitment, egg biomass is lost from the adult insects and is transferred to the Instars 1–2 biomass, which is, in other words, the biomass gained from successful hatching.

If T (t) > (0.8 × TOpt) and T (t) < (TOpt − 1) then

$$Recruit\left(t\right)=PctEgg\times Emerg\left(t\right)\times 2$$

(4)

$$\mathbf{Else}:Recruit\left(t\right)=0$$

where Recruit (g/m³.d) is the recruitment from viable eggs to Instars 1–2, PctEgg (unitless) the fraction of adult biomass that is in the carried eggs and Emerg (g/m³.d) the emergence from Instars 3–4 to adult insects.

Finally, the model calibration was conducted following the same stepwise procedure as the one used in reference [12]. Chaoborus was added to the previously calibrated model and the AQUATOX default parameter values were selected as starting point for the model calibration [24]. The biomass concentrations measured in the experimental lake on 2 May 2000, were used as initial biomass concentrations and were not changed during calibration. Finally, a sensitivity analysis was performed, and the most influential parameters tuned.

2.3. Step 3: Simulating Multitrophic Interactions

The objective of this paper was to understand the direct and indirect effects of EE2 on a freshwater ecosystem by creating links between aquatic biodiversity and ecological functioning. In order to do that, the role of seasonal cycle alterations on the food web, resulting in an increase or decrease in number of a specific population, needs to be studied.

According to Duffy et al. [19], understanding how biodiversity affects functioning of ecosystems requires integrating diversity within trophic levels (horizontal diversity) and across trophic levels (vertical diversity). Distinguishing these two dimensions can help clarify how ecosystem functioning may be affected separately or simultaneously by consumptive interactions across trophic levels and competitive processes within levels. Besides, it is recommended to group species that perform similar roles in an ecosystem process into “functional groups” [24].

In the ecosystem model presented in this paper, species were divided based on their trophic status (e.g., their place in the food web as producers, decomposers and predators) and their seasonal cycles (growth bloom for plankton, spawning for fish or emergence for insects). Then, a global sensitivity analysis similar to reference [25] was used to better understand the multitrophic interactions occurring in the studied aquatic ecosystem.

The sensitivity of populations was evaluated by multiplying (by five) or dividing (by five) the initial biomass of a specific functional group. The consequences on the entire food web were then interpreted from the biomass dynamics of the other groups. For example, the initial biomass of rotifers was divided by five and then the biomass dynamics over time of the other zooplankton (horizontal diversity) and other trophic levels (vertical diversity) were analyzed. The functional groups selected for the model (Figure 1) were studied one by one.

To allow the comparison of all simulations, increasing and decreasing the initial biomass of every single functional group, the resulting simulated biomass were normalized according to the averaged biomass obtained over time. Finally, the above simulations were performed with and without Chaoborus in the food web. Indeed, since the ecosystem model presented in this paper was extended by adding Chaoborus to the food web, it is important to assess its role in the studied freshwater ecosystem.

2.4. Step 4: Predicting EE2 Indirect Effects

EE2 was added to the experimental lake at environmentally relevant concentrations; 4.8–6.1 ng/L in the epilimnion and around 2 ng/L in the hypolimnion (decrease due to sorption on sediments and bioaccumulation) [11]. In this whole-ecosystem experimental study, evidence of direct effects of EE2 on the abundance of fishes was found but little evidence could be found of direct effects of this synthetic estrogen on lower-trophic-level organisms [13]. The authors also observed the most notable negative indirect effects on the food web were declines in the biomass of the top predator lake trout due to a decline in its food supply (fathead minnow, pearl dace and juveniles of white sucker). Besides, while the increases in diverse invertebrate prey species were small, interpreted together they provide a strong indication of indirect effects on the lower food web due to observed declines in predation by small-bodied fishes.

Therefore, in this paper, the extended ecosystem model is used to confirm and take further the hypotheses on EE2 indirect effects on a freshwater food web [13]. First, the direct EE2 effects on fathead minnow that [12] simulated as a decrease in gamete production, increase in gamete mortality and/or increase in fish mortality, need to be confirmed as potential hypotheses with the new model including Chaoborus. Then, EE2 indirect effects on the other fishes (horizontal diversity) and the lower trophic levels (vertical diversity) will be interpreted with the objective to identify a top-down or bottom-up effect. Finally, indirect EE2 effects will be compared with and without Chaoborus in the food web.

3. Results

3.1. Extended Ecosystem Model

The addition of Chaoborus to the ecosystem model built and calibrated in [12] created new multitrophic interactions. Indeed, Chaoborus feed on POM, phytoplankton C (chlorophytes, dinoflagellates and Cyanophyta) and zooplankton (rotifers, Cladocera and copepods), while fish (fathead minnow and pearl dace) feed on Chaoborus (Figure 3). The most influential parameters for Chaoborus are presented in S1 (Supplementary Material: Table S1). Those values were obtained after the model was fit to the experimental data, by fine-tuning of the parameters calibrated in [12]. Then, the ecological experts validated the simulated biomass dynamics, also including the two lower data values at the end of September (Figure 5).

It can be noted on the graph that, in the spring, only instars 3–4 are present in the experimental lake. Starting mid-June, the beginning of insect emergence, they are leaving the aquatic environment for the terrestrial ecosystem, as can be seen with a drop of the instars 3–4 biomass. In July, the instars 1–2 biomass starts increasing because the adult insects are releasing eggs into the lake, which turn into larvae. Mid-July, instars 1–2 biomass decreases while instars 3–4 biomass increases because of promotion. Finally, instars 3–4 biomass keep increasing due to larvae growth, up to the end of September when temperatures and food availability start to fall and consequently, biological metabolisms slow down.

The most influential parameters presented in S1 (Table S1) for Chaoborus are mainly related to food consumption, with some parameters for reproduction, growth, emergence, excretion, respiration and mortality. The addition of Chaoborus also had an impact on calibrated values of some parameters of the functional groups connected to Chaoborus through multitrophic interactions (for more details, see Tables S2–S4 in S1).

After fine tuning of the calibrated parameters, the simulation graphs obtained for all trophic levels were compared with the ones in Clouzot and Vanrolleghem [12]. No significant changes can be noticed after Chaoborus was added to the model. In the example of zooplankton given in Figure 6, it clearly appears that the trend is similar to the one obtained with the previous ecosystem model. When comparing the Mean Square Residual Error (MSRE) values between the two ecosystem models, the simulation results were actually closer to the experimental data with the extended model (MSRE values were lower).

In conclusion, the ecosystem model developed in Clouzot and Vanrolleghem [12] adapted well to the addition of Chaoborus, with parameter values mostly changed based on the new trophic interactions. In other words, the extended ecosystem model can be used with confidence to further study the role of seasonal cycle alterations on the food web, as an aspect of the relationship between biodiversity and ecosystem functions.

3.2. Seasonal Cycles, Multitrophic Interactions and Biodiversity

For each trophic level, the simulated time-series results of the calibrated model are presented. These results correspond to the best fit of the calibrated model with the observations presented in [12]. Consultation with the scientists who gathered the experimental data confirmed that the calibrated model succeeded in simulating relatively well the main ecosystem trends observed during the open water season.

In order to assess the role of seasonal cycle alterations on the food web, the individual seasonal cycles of each modeled functional group are added to the simulation graphs and the impact of multitrophic interactions discussed. The results from the global sensitivity analysis are used to help better understand how ecosystem functioning may be affected separately or simultaneously by competitive processes within trophic levels (horizontal diversity) and consumptive interactions across trophic levels (vertical diversity).

3.2.1. Primary Producers—Phytoplankton

Seasonal Cycles and Multitrophic Interactions

For the three modeled groups of phytoplankton, algae blooms anticipated by experts are indicated on the simulation graph (Figure 7). It can be noted that actual algae blooms did not appear on the simulation graph (no significant increase of biomass), which is consistent with the experimental data.

Phytoplankton are at the bottom of the food web and are eaten by numerous predators (i.e., zooplankton, Chaoborus, fathead minnow and pearl dace). Therefore, a working theory can be that the increase of phytoplankton biomass during a bloom leads to more food for the predators, increasing their consumption rate. Due to those multitrophic interactions happening in freshwater ecosystems, the effects of the phytoplankton seasonal cycles (algae blooms) are not observed in the experimental data and also are not present in the model simulation. The physical mechanisms that would lead to blooms in the model would be a change in temperature, light and/or nutrient availability, which is overcome here by predatory pressure. Those results confirm that both individual seasonal cycles and multitrophic interactions need to be integrated into the conversation about biodiversity and ecosystem functioning (more details in Supplementary Material S2).

Horizontal Diversity

The global sensitivity analysis of phytoplankton shows that each phytoplankton group biomass had an indirect impact upon the other phytoplankton groups (i.e., diatoms in Figure 8). When their initial biomass concentration was increased, the biomass concentrations of the other phytoplankton increased as well. In that scenario, the group for which the biomass was increased became the main food source for their predators and the other phytoplankton were eaten in less quantity (which explains their increase in biomass). On the contrary, when the phytoplankton group for which the initial biomass was decreased became a limited food source, the predators eat the other phytoplankton (which explained their decrease in biomass).

Vertical Diversity

During the global sensitivity analysis of Chaoborus, the only phytoplankton impacted were the ones consumed by the insects (group C), with a significant decrease after the initial biomass concentration of Chaoborus was increased. This result highlights consumptive interactions across trophic levels and predatory pressure.

Regarding the results obtained for zooplankton, the effects on phytoplankton were more complex, probably because different zooplankton species interact differently with phytoplankton. When the initial biomass concentration of any zooplankton group was increased (see phyto A++, diatoms++ and phyto C++ in Figure 9), all phytoplankton increased in biomass. Rotifers were chosen as an example because they were not a food source for zooplankton in the studied system and yet they were impacted in the same way by a change in zooplankton biomass concentrations. When looking closer to the other results of the global sensitivity analysis, it can be noted that phosphorus (P) concentrations increased as well (from 0.001–0.004 to 0.002–0.007 g/m³), which could explain the increase of phytoplankton (for nitrogen, see Clouzot and Vanrolleghem [12]).

Another interesting result is that when the initial biomass concentration of any zooplankton group was decreased, differences appeared between diatoms and phytoplankton A/C. Indeed, diatoms biomass concentrations decreased (see diatoms—in Figure 9) while no significant changes were observed for group A/C (not presented in Figure 9). This result could be connected to the notion of functional groups explained previously (Section 2.3.), which suggest that species performing similar roles in an ecosystem process are grouped together. When looking at the calibrated parameters used in the model for phytoplankton (see [12]), it seems that both group A and C assimilate better P than nitrogen (N; in terms of nutrient affinity), while diatoms assimilate P and N the same. This difference in assimilation could be the reason why the phytoplankton groups were not impacted in the same way after zooplankton biomass was decreased.

Finally, no significant effects were observed on phytoplankton during the global sensitivity analysis of the modelled fish species (fathead minnow, pearl dace, white sucker and lake trout). The possible hypotheses are that if there is a top-down effect within the studied freshwater ecosystem, it does not go down to the bottom of the foodweb, as compensatory mechanisms are responsible for keeping overall phytoplankton biomass the same.

3.2.2. Primary Consumers—Zooplankton

Seasonal Cycles and Multitrophic Interactions

For the three modeled groups of zooplankton, blooms anticipated by experts are indicated on the simulation graph (Figure 10). Similar to phytoplankton, the rotifers bloom did not appear on the simulation graph, which was consistent with the experimental data. Rotifers were only eaten by Chaoborus, which increased in biomass in August (Figure 5), during the rotifers bloom. Therefore, the food consumption of Chaoborus increased, which can explain why the peak for rotifers bloom did not appear on the simulation graph.

With regards to copepods and Cladocera, the blooms anticipated by experts were also not present in the simulation results. Instead, biomass concentrations decreased steadily. Copepods and Cladocera are food sources for Chaoborus, fathead minnow, pearl dace and white sucker. Besides, they feed on phytoplankton, which have low biomass concentrations. Consequently, their individual seasonal cycles were hidden by the multitrophic interactions occurring in the freshwater ecosystem.

Horizontal Diversity

The global sensitivity analysis of the zooplankton simulation results shows that the three groups had a similar impact on the other groups of the same trophic level (i.e., Cladocera in Figure 11). For example, in the case of Cladocera, with their initial biomass concentration as the one in the experimental lake, the copepods biomass concentrations were first higher and then lower (Figure 11, see copepods). With less Cladocera, the copepods become the main food source for predators feeding on zooplankton, but they also get more food because of reduced competition with Cladocera. Until mid-July, copepods are in their bloom phase, so the increased predatory pressure did not have an impact yet, the food availability having a stronger effect. However, when their anticipated bloom is over, the predatory pressure is stronger, and their biomass concentrations decrease. Similar effects can be noted for rotifers but at a lower intensity (Figure 11, see rotifers).

In the scenario of increasing the initial biomass concentration of Cladocera, the copepods concentrations increase until mid-July and then decrease (Figure 11, see copepods ++). Cladocera became the main food source for predators and thus, copepods could grow until their bloom was done (mid-July). Once again, the effects were lower for rotifers, but they lasted longer (Figure 11, see rotifers ++).

Vertical Diversity

During the global sensitivity analysis of Chaoborus, significant effects on zooplankton were only observed when the initial biomass concentration was increased (see rotifers++, Cladocera++ and copepods++ in Figure 12). In that case, the predatory pressure of Chaoborus was higher, which explains the decrease in concentration of the three zooplankton groups.

With regards to phytoplankton, since they are a food source for zooplankton, changes in any phytoplankton biomass concentrations directly impact zooplankton biomass dynamics. Higher concentrations of zooplankton were observed when any phytoplankton initial biomass concentration was increased. On the opposite, less abundant phytoplankton resulted in lower concentrations of zooplankton. Finally, changes in any fish initial biomass concentration did not result in significant changes in zooplankton biomass dynamics, which could be explained by compensatory mechanisms or migration of zooplankton to avoid potential increased predatory pressure.

3.2.3. Secondary and Tertiary Consumers—Fish

Seasonal Cycles and Multitrophic Interactions

For the four modeled fish (fathead minnow, pearl dace, white sucker and lake trout), anticipated spawning periods are indicated on the simulation graph (Figure 13). The same biomass concentration pattern can be noted, with an increase of adult biomass before spawning, followed by a decrease during spawning, along with an increase of juvenile biomass. An increase of biomass also appeared in the summer, due to high temperatures and food availability, and then a decrease in the fall, due to a decrease in temperature and food availability. Fish being at the top of the food web and feeding on lower trophic levels, individual seasonal cycles clearly appeared on the simulation graphs, along with effects of ecological interactions.

Those same graphs (Figure 13) can also be interpreted with multitrophic interactions in mind, adding a layer of complexity inherent to any freshwater ecosystem. First, the initial increase of fathead minnow biomass concentrations can be linked to high concentrations of its preys at the same time (Cladocera and copepods) and the second increase after spawning, to high concentrations of another prey, Chaoborus. The decrease at the end can then be explained by a decrease in biomass of its prey (Chaoborus) and an increase of its predator (lake trout). Similar observations can be made for pearl dace, with spawning happening earlier. Next, the increase in lake trout biomass concentrations after August can be connected to the increase of all its prey (fathead minnow, pearl dace and juvenile white sucker). Finally, white sucker did not seem to be influenced much by the rest of the ecosystem, with biomass concentrations being relatively constant through the whole open water season.

Horizontal Diversity

In the modeled ecosystem, the secondary consumers were fathead minnow, pearl dace and white sucker, while lake trout was a tertiary consumer. Therefore, horizontal diversity could only be studied at the level of secondary consumers. The global sensitivity analyses of fathead minnow and pearl dace both demonstrated similar effects on each other (i.e., fathead minnow in Figure 14). For example, a decrease of the fathead minnow initial biomass concentration resulted in lower concentrations of pearl dace (see pearl—in Figure 14), which can be explained by a change of predatory pressure. Since lake trout eats interchangeably fathead minnow and pearl dace, a reduced fathead minnow biomass leads lake trout to eat more pearl dace, resulting in a decrease of their biomass. In the case of an increase of fathead minnow concentrations, less food is available for pearl dace, which leads to a decrease of their concentrations (see pearl ++ in Figure 14). In conclusion, a combination of predatory pressure and food availability can explain the observed impact on horizontal diversity. Similar effects were observed for the juveniles of white sucker, which are another food source for lake trout. Conversely, because lake trout did not feed on the adults of white sucker, no significant changes were observed when changing the initial biomass of fathead minnow or pearl dace.

Vertical Diversity

During the global sensitivity analyses of phytoplankton and zooplankton, no significant effects were observed on the secondary consumers (fathead minnow, pearl dace and white sucker), except for juvenile white sucker, with lower concentrations when one of its preys (copepods) decrease. Both fathead minnow and pearl dace had a diversified diet, which gave them a better resilience when changes happened in the ecosystem. On the other hand, significant effects were observed for fathead minnow, pearl dace and juvenile white sucker when changes were applied to their only predator, lake trout. When the lake trout initial biomass concentration was increased, its preys’ concentrations decreased, and vice-versa.

Other significant changes happened when the initial biomass concentration of Chaoborus was increased (Figure 15), as a result of multitrophic cascade effects. First, lower concentrations of fathead minnow were observed (same effect with pearl dace; see fathead [A]/[J]++ in Figure 15). Due to higher concentrations of Chaoborus, the zooplankton biomass was lower due to increased predatory pressure (Figure 11) and they became a limited food resource for the small fishes (fathead minnow and pearl dace). At the same time, a strong decrease in white sucker biomass can be noticed because of lower food resources and higher predatory pressure from lake trout on the juveniles (see white [A]/[J]++ in Figure 15). Finally, the lake trout biomass concentrations also decreased due to its preys becoming limited (see trout [A]/[J]++ in Figure 15). These results demonstrate a strong connection between the secondary consumers but also across the different trophic levels. Besides, Chaoborus appeared as a key species in the studied system, where changes in its biomass concentrations impacted the whole ecosystem.

3.2.4. Omnivorous, Mid-Trophic Level Group—Chaoborus

Seasonal Cycles and Multitrophic Interactions

The simulation graphs for Chaoborus (Figure 5), with the different phases highlighted, were previously discussed in Section 3.1, along with the multitrophic interactions involving Chaoborus. Noteworthy, when comparing the simulation graphs of Chaoborus and fathead minnow, similar trends were observed (Figure 5 and Figure 13a). As with fish, individual seasonal cycles of Chaoborus appeared on the simulation graphs.

Vertical Diversity

All zooplankton and phytoplankton C (chlorophytes, dinoflagellates and cyanophytes) are a direct food source for Chaoborus and thus, followed the same trend. For diatoms and phytoplankton A (chrysophytes and cryptophytes), the same direct effect on Chaoborus was observed because they had a direct effect on phytoplankton C (Figure 8). A reverse effect was observed with fathead minnow and pearl dace (fathead minnow in Figure 16b) because, in this interaction, Chaoborus was the food source. Again, multitrophic interactions were directly connected to the ecosystem dynamics, and thus its functioning.

When comparing the results of the global sensitivity analysis of the previous model [12], which did not consider Chaoborus, and the extended ecosystem model presented in this paper, it seemed that the inclusion of insects created higher resilience for the freshwater ecosystem. For example, the sensitivity analysis of Cladocera shows the effects on white sucker and lake trout were lower in intensity when Chaoborus was added to the modeled ecosystem (Figure 17). Similar results were obtained for the different functional groups selected for the model. Therefore, while adding an omnivorous, mid-trophic level group such as Chaoborus, which is connected to many other trophic levels, added complexity to freshwater ecosystems, it could also lead to an enhanced resilience of the ecosystems when changes were applied.

3.3. Indirect Effects of EE2

The experimental results obtained from Lake 260 showed that the strongest direct effect of EE2 was on fathead minnow, with a collapse of the fish species in the second year of adding EE2 to the lake [11]. Discussions with experts in ecotoxicology and endocrine disruption helped identify the parameters that should be modified in the ecosystem model. It was previously decided that a mix of reduced gamete production, increased gamete mortality and an increase of both adults and juvenile mortality was a potential hypothesis for explaining the collapse of both adult and juvenile fathead minnow due to EE2 addition [12]. Therefore, the same combination of parameters was applied to the extended ecosystem model.

In comparison with reference [12], the addition of Chaoborus to the freshwater ecosystem did not impact the EE2 indirect effect observed on the other fish species, i.e., a decrease in biomass of pearl dace and white sucker, accompanied with a change in lake trout biomass. Since lake trout could no longer feed on fathead minnow, they turned to their other food sources.

Regarding Chaoborus itself, both the simulation results and the experimental data for instars 3 and 4 after adding EE2 to the lake show higher Chaoborus concentrations than the results without the hormone (Figure 18a). The increase in Chaoborus biomass can be explained by the decrease in fish predatory pressure due to the collapse of the fathead minnow population. Similar effects were observed with Cladocera and copepods (Figure 18b), but at a lower degree due to the additional predatory pressure from the Chaoborus population.

For rotifers, no significant changes in biomass concentrations were observed in the ecosystem model simulations after adding EE2 (results not shown). However, it was found in the experimental lake that the average rotifer biomass did increase, even if they do not constitute a major part of the diet of the fish species impacted by EE2 [13]. The authors raised the hypothesis of complex multitrophic interactions that could explain this result, such as changes in the diet of zooplankton predators. Intraguild competition and predation among omnivorous zooplankton predators including minnows, Chaoborus and cyclopoid copepods may lead to complex interactions. For example, the effect of decreased fish predation may have been partially offset by increases in Chaoborus predation on other invertebrates. Declines in Tropocyclops, which feed on rotifers, or changes in the behavior of other invertebrate predators may also have affected rotifers and other small taxa. Therefore, one of the reasons for the rotifers increase not being captured by the simulation results may be that the model cannot simulate any changes of diet or habitat due to changes in the ecosystem, such as resource availability or predatory pressure. Besides, a simplification of the model was to consider particulate detritus as the main diet of rotifers, which is relatively close to reality but not entirely accurate.

Finally, the simulation results before and after adding EE2 confirm there are no indirect effects of the synthetic hormone on phytoplankton. This result is consistent with the observation made during the global sensitivity analysis of fish and the vertical impact on phytoplankton (Section 3.2.1). The top-down effect confirmed for EE2 on the freshwater ecosystem thus does not go down to the bottom of the freshwater food web, showing the resilience of this trophic level.

4. Discussion

4.1. Chaoborus Contribution to Ecosystem Functioning

A dependence diagram of the extended ecosystem model was created to better illustrate the Chaoborus contribution to the freshwater ecosystem functioning (Figure 19). Each trophic level is identified with a specific color, such as blue for fish. Biomass transfer through food consumption is represented with arrows and percentages indicate the fraction of the overall food consumption for each functional group. For instance for Chaoborus, its diet was composed of 10% rotifers, 40% Cladocera and copepods, 40% phyto C (chlorophytes, dinoflagellates, cyanophytes) and 10% particulate detritus (PD).

Energy and materials move between the biotic (i.e., species) and the abiotic (i.e., inorganic and organic nutrients) compartments, as well as into and out of the system. Ecosystem processes are quantified by measuring rates of these movements (e.g., plant production), while ecosystem functioning is quantified by measuring the magnitudes and dynamics of ecosystem processes. It has been shown that most ecosystem processes are driven by the combined biological activities of many species [24]. Additionally, because species can vary dramatically in their contributions to ecosystem functioning, the specific composition or identity of species in a community is important.

One way to describe the average number of times energy is transferred as it moves from basal resources to top predators is the food chain length (FCL), directly related to ecosystem functioning [19]. Since horizontal and vertical diversity interact, adding the omnivorous (i.e., feeding from more than one trophic level) Chaoborus to the food chain can qualitatively change diversity effects at adjacent levels (Figure 19). For instance, Hanazato [26] showed that in ponds containing the zooplankton predator Chaoborus, the zooplankton community structure became dominated by rotifers. In ponds without Chaoborus, zooplankton became dominated by cladocerans (a taxon that is competitively superior to rotifers).

Another interesting result is presented in Hanazato [27]: novel aspects of Chaoborus in driving ecosystem functioning in lakes and reservoirs. Chaoborus larvae seem to account for a significant fraction of both the hypolimnetic and sediment oxygen demands, and effectively trap nutrients in water and sediment, and in this way they enhance internal nutrient loading.

Finally, the simulation results obtained with the extended ecosystem model demonstrated that adding an omnivorous, mid-trophic level group such as Chaoborus, enhanced resilience (Figure 17). Conversely, missing such a species means that the actual resilience of an ecosystem and its functioning cannot be properly simulated. The stability of an ecosystem corresponds to its ability to maintain a comparable functioning in the presence of perturbations that drive it away from its original state [28]. Resilience, in turn, is the ability or the time taken by a system to recover from a change due to perturbations and is linked to stability.

4.2. EE2 Effects on Ecosystem Functioning

The EE2 indirect effects on the studied freshwater food web, summarized in a dependence diagram (Figure 20), were connected to multitrophic interactions, such as food resources availability, predatory pressure and compensatory mechanisms. It clearly appeared that the removal of fathead minnow from the food web changed the dynamics. For example, before EE2 addition, 80% of the lake trout’s diet consisted of fathead minnow and pearl dace. After EE2 was added, 80% of the lake trout diet was now only obtained from pearl dace, which added predatory pressure on this species, resulting in a decrease of the biomass. Following a similar logic, Chaoborus benefits from the removal of fathead minnow. While without EE2, Chaoborus represented 40% of the diet of both fathead minnow and pearl dace, it was now only consumed by pearl dace (40% of its diet), which explains the increase in Chaoborus biomass.

Those results aligned with the well understood drivers of the biodiversity and ecosystem functioning (BEF) link: selection effects, complementarity in resource use and species interactions [21]. Resources may be directed towards pathways that enable persistence in the new environment at the expense of other organismal processes, like growth, reproduction and foraging behavior. Changes in organismal processes consequentially impact upon the ecosystem processes they mediate, like primary production.

In the case of the synthetic hormone EE2, the main direct observed effect was the collapse of fathead minnow [13]. Such species at higher trophic levels nearly always reveal effects that span through the food web. However, the magnitude and direction of these effects are highly variable and are difficult to predict since these species exhibit many complex, indirect, non-additive and behavioral interactions [20]. The simulation results obtained from the extended ecosystem model confirmed a top-down cascading effect of EE2 on the freshwater food web, as highlighted in red in Figure 20. In their review, Duffy et al. [19] highlighted growing evidence that cascading impacts of predators on primary producers often occur through trait-mediated indirect effects, specifically by modifying behavior rather than via changes in herbivore density. Besides, since large predators are naturally low in species diversity, a few extinctions may result in loss of the entire top predator trophic level, with disproportionately large effects on ecosystem properties and processes.

Trophic cascades refer to the phenomenon that top predators suppress the biomass of intermediate consumers, which in turn releases consumers at the next lower trophic level from predation pressure, with this alternate suppression-and-release effect propagating down to the producer level. In the case of EE2, the main direct effect was observed on fathead minnow, one of the species at higher trophic levels, with cascaded effects to the top predator (i.e., lake trout) and to lower trophic levels (down to zooplankton). Jabiol et al. [29] suggested that multitrophic interactions complicate assessments of the functional consequences of biodiversity change and that even relatively modest effects of change within single trophic levels can ramify and amplify through ecological networks to influence ecosystem functioning.

Finally, increasing evidence is found that the key means by which species influence ecosystem functions (EFs) is through their functional traits (phenotypic attributes that direct niche exploitation) [30]. It has been hypothesized that beyond the lower end of a species richness gradient, the main driver of EFs is the community functional structure. Indeed, Chair et al. [24] highlighted the ability of competing species to replace or compensate for one another and thus minimize, at higher diversity, the ups and downs in functioning. Such observations are confirmed by the simulation results obtained with the extended model presented here, where a clearly higher ecosystem resilience was reached when including Chaoborus, i.e., increasing this ecosystem’s biodiversity.

4.3. In the Context of Ecosystem-Based Management (EBM)

Unlike the biodiversity to ecosystem functions (BEF) relationship, Teixeira et al. [15] shows there is less established biodiversity to ecosystem services (BES) research. Whether or not biodiversity benefits from the protection of ES, and vice versa, the authors consider BES as valuable scientific knowledge to turn the concept of ecosystem services into a practical conservation tool in the formulation of day-to-day policies at national or regional scales. Indeed, the simulation results of the extended ecosystem model obtained after adding EE2 to the studied lake suggested potential species movement and behavior changes in response to variations of food resource and predatory pressure, which can, in turn, influence ecosystem functioning and the provision of ecosystem services [31].

As Queirós et al. [21] pointed out, ecosystem models are synthetic mathematical descriptions of ecosystem processes joined together, guided by a mechanistic understanding of their regulating environmental drivers and biota, which can be used to project changes in the overall ecosystem properties. In this way, ecosystem models provide a platform where empirical findings can be used to investigate large-scale questions and to offer a holistic view of ecosystems where the impacts of conservation, management and global scenarios can be assessed. The simulation results obtained with the extended ecosystem model demonstrate that adding an omnivorous, mid-trophic level group such as Chaoborus, enhances resilience. Conversely, missing such species means that the actual resilience of an ecosystem and its functioning cannot be properly simulated.

In conclusion, the extended ecosystem model is a reliable and robust tool that can be used to simulate what is happening in freshwater food webs after perturbations, such as the addition of the synthetic hormone EE2. This can help better understand the link between biodiversity and ecosystem functioning in the context of EBM.