# Microcosm experiment combined with process-based modeling reveals differential response and adaptation of aquatic primary producers to warming and agricultural run-off

Gregorio López Moreira Mazacotte, Bastian Polst, Elisabeth Gross, Mechthild

Schmitt-Jansen, Franz Hölker, Sabine Hilt

► To cite this version:

Gregorio López Moreira Mazacotte, Bastian Polst, Elisabeth Gross, Mechthild Schmitt-Jansen, Franz Hölker, et al.. Microcosm experiment combined with process-based modeling reveals differential re-sponse and adaptation of aquatic primary producers to warming and agricultural run-off. Frontiers in Plant Science, 2023, 14 (1120441), 10.3389/fpls.2023.1120441. hal-04133190

HAL Id: hal-04133190

https://hal.univ-lorraine.fr/hal-04133190

Submitted on 19 Jun 2023

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
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Jun Xu,
Institute of Hydrobiology (CAS), China

REVIEWED BY
Monica Pinardi,
Institute for Electromagnetic Sensing of the
Environment, Italy
Ingo Dreyer,
University of Talca, Chile

*CORRESPONDENCE
Gregorio A. López Moreira Mazacotte
✉ gregorio.lopezmoreira@igb-berlin.de
Bastian H. Polst
✉ bastian.polst@gmail.com

†These authors have contributed
equally to this work and share
first authorship

RECEIVED 10 December 2022

ACCEPTED 03 April 2023

PUBLISHED 19 June 2023

**CITATION**

López Moreira Mazacotte GA, Polst BH, Gross EM, Schmitt-Jansen M, Hölker F and Hilt S (2023) Microcosm experiment combined with process-based modeling reveals differential response and adaptation of aquatic primary producers to warming and agricultural run-off.

*Front. Plant Sci.* **14**:1120441.

doi: 10.3389/fpls.2023.1120441

**COPYRIGHT**

© 2023 López Moreira Mazacotte, Polst, Gross, Schmitt-Jansen, Hölker and Hilt. This is an open-access article distributed under the terms of the [Creative Commons Attribution License (CC BY)](https://creativecommons.org/licenses/by/4.0/). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Microcosm experiment combined with process-based modeling reveals differential response and adaptation of aquatic primary producers to warming and agricultural run-off

Gregorio A. López Moreira Mazacotte¹*†, Bastian H. Polst²*†,
Elisabeth M. Gross³,⁴, Mechthild Schmitt-Jansen²,
Franz Hölker¹ and Sabine Hilt¹

¹Department of Community and Ecosystem Ecology, Leibniz Institute of Freshwater Ecology and Inland Fisheries (IGB), Berlin, Germany, ²Department of Bioanalytical Ecotoxicology, Helmholtz Centre for Environmental Research-UFZ, Leipzig, Germany, ³Laboratoire Interdisciplinaire des Environnements Continentaux (LIEC) UMR 7360 CNRS, Université de Lorraine, Metz, France, ⁴LTSER Zone Atelier Bassin de la Moselle, Metz, France

Fertilizers, pesticides and global warming are threatening freshwater aquatic ecosystems. Most of these are shallow ponds or slow-flowing streams or ditches dominated by submerged macrophytes, periphyton or phytoplankton. Regime shifts between the dominance of these primary producers can occur along a gradient of nutrient loading, possibly triggered by specific disturbances influencing their competitive interactions. However, phytoplankton dominance is less desirable due to lower biodiversity and poorer ecosystem function and services. In this study, we combined a microcosm experiment with a process-based model to test three hypotheses: 1) agricultural run-off (ARO), consisting of nitrate and a mixture of organic pesticides and copper, differentially affects primary producers and enhances the risk of regime shifts, 2) warming increases the risk of an ARO-induced regime shift to phytoplankton dominance and 3) custom-tailored process-based models support mechanistic understanding of experimental results through scenario comparison. Experimentally exposing primary producers to a gradient of nitrate and pesticides at 22°C and 26°C supported the first two hypotheses. ARO had direct negative effects on macrophytes, while phytoplankton gained from warming and indirect effects of ARO like a reduction in the competitive pressure exerted by other groups. We used the process-based model to test eight different scenarios. The best qualitative fit between modeled and observed responses was reached only when taking community adaptation and organism acclimation into account. Our results highlight the importance of considering such processes when attempting to predict the effects of multiple stressors on natural ecosystems.

**KEYWORDS**

climate change, regime shifts, alternative stable states, stress-induced tolerance, mathematical modeling, multiple stressors, shallow lakes, phytoplankton succession
# 1 Introduction

The overall share of land used for crops and pastures is increasing worldwide (Winkler et al., 2021). However, global cropland per capita is steadily decreasing as the world's population continues to grow (FAO, 2021), leading to a rapid increase in the global use of pesticides and fertilizers (Sharma et al., 2019). In addition, climate change is increasing the average temperature of most inland water bodies (O'Reilly et al., 2015). Together, these multiple stressors may severely affect aquatic ecosystems in agricultural areas. Many of these agroecosystems are shallow ponds or slow-flowing streams and are characterized by the alternative dominance of submerged macrophytes, periphyton or phytoplankton, competing for light and nutrients as major resources (Vasconcelos et al., 2016). Shifts from submerged macrophyte dominance to the less desirable phytoplankton dominance have been reported with increased nutrient loading (Sayer et al., 2010). Warming is expected to lead to more severe phytoplankton dominance (modeling study by Mooij et al., 2007) and to increase periphyton growth (Mahdy et al., 2015; Kazanjian et al., 2018), resulting in an overall weakening of the resilience of macrophyte-dominated systems to nutrient loading impacts (Meerhoff et al., 2022). Additional stressors such as pesticides, however, may antagonistically interact with temperature and nutrient loading, modulating the competition for resources between aquatic primary producers (Polst et al., 2022). Understanding the overall response of these ecosystems to multiple stressors is therefore challenging (Jackson et al., 2016) but crucial to predict future changes in their functioning (Hilt et al., 2017) and services (Janssen et al., 2021).

A first experiment investigating responses of the different primary producers typical for shallow aquatic systems to multiple stressors revealed a direct positive effect of co-occurring nitrate and pesticides on phytoplankton, with no significant effect of warming (Allen et al., 2021). A subsequent study, involving also primary consumers, showed that pesticides and nitrate may interact synergistically to reduce macrophyte dominance, and that pesticides and warming may have especially strong but opposite effects on specific macrophyte species (Vijayaraj et al., 2022). Climate warming has even been shown to lower critical thresholds for shifts in dominance between aquatic primary producers (Polst et al., 2022). These differences may derive from the differential sensitivities to toxicants of different groups and species of primary producers (Fairchild et al., 1998; Giddings et al., 2013), some of which may be more tolerant or have the ability to acclimate or adapt to them at different temperature-dependent rates (Chalifour and Juneau, 2011; Larras et al., 2013). Additional trophic levels, the structure of the ecosystem and the trophic status of the water body may also play a role (Wendt-Rasch et al., 2004), further complicating the analysis of direct and indirect effects of multiple stressors that may lead to complex responses of shallow aquatic ecosystems.

Reductionist factorial laboratory approaches are useful to address related research questions. However, time and budget constraints, and the invasive nature of macrophyte sampling often limit both the number of testable stressor combinations and the frequency of sampling that would be needed to understand

community dynamics. Combining experiments with process-based models may be a helpful tool to reveal variable trajectories, to test for critical threshold values, and to disentangle indirect stressor effects. Such models, however, need to be developed based on the specific research questions and the experimental set-up, i.e., custom-tailored to avoid running into equifinality issues (López Moreira et al., 2021). The empirical results, in turn, allow for subsequent model refinement, calibration and validation (e.g., Kim et al., 2021), an iterative process that builds upon every new insight and gained expert knowledge (Jakeman et al., 2006).

In this study, we combined a microcosm experiment on multiple stressor effects in shallow aquatic ecosystems with simulations we ran with a custom-tailored process-based model. We developed this model to investigate the response of the different primary producer groups (submerged macrophytes, periphyton and phytoplankton) to the combined effects of nitrate and pesticides typically found in agricultural run-off (ARO) with or without effects of climate warming. We hypothesized: 1) that ARO differentially affects primary producer groups increasing the likelihood of phytoplankton dominance in shallow aquatic ecosystems; 2) that warming facilitates this process and 3) that custom-tailored process-based models can support mechanistic understanding of experimental results through scenario comparison. First, we developed a simple process-based mathematical model based on the known effects of herbicides, fertilizers and warming on phytoplankton, periphyton and macrophytes. To assess model predictions, we carried out a microcosm experiment testing the response of three submerged macrophyte species (*Myriophyllum spicatum*, *Potamogeton perfoliatus* and *Elodea nuttallii*) and mixtures of phytoplankton and periphyton typically occurring in European freshwaters. We exposed these systems to a gradient of an experimental ARO cocktail, an artificial mix of organic pesticides (an herbicide, an insecticide and a fungicide), copper-(II) sulfate ($CuSO_4$) and potassium nitrate ($KNO_3$) for 19 days. The ARO mix was applied in a two-factor factorial design (dose-response set-up) at two different temperatures (ambient: 22°C, warming: 26°C). Experimental results led us to refine the process-based model and test eight scenarios of differential herbicide sensitivities among primary producer groups, temperature dependence of the response, and development of tolerance to the herbicide for three cases (sets) of phytoplankton community composition. Comparing *in silico* simulations with experimental results allowed us to select the best fit and improved mechanistic understanding of differential sensitivities of phytoplankton, periphyton and submerged macrophytes to combined ARO and warming.

# 2 Materials and methods

## 2.1 Process-based model to simulate the microcosms

### 2.1.1 Governing and supplementary equations

To simulate the combined effects of ARO and warming on the different groups of aquatic primary producers, we developed a process-
based mathematical model, which we implemented in MATLAB (R2020a). The model comprises a series of ordinary differential equations (ODEs) similar to those proposed in previous works (Jäger et al., 2010; Vasconcelos et al., 2016). These equations describe the following state variables: the live carbon (C) biomass of phytoplankton (as volume concentration, eq. 1, specific to each phytoplankton group), periphyton, hereafter understood as the growth on the vertical inner surface of the microcosm vases (as areal density, eq. 2), epiphyton, i.e., periphyton growing on macrophyte surfaces (as areal density, eq. 3) and macrophytes (eq. 4). Because exudates and lysates resulting from growth and senescence, respectively, are known to support microbial communities (Kieft et al., 2021), the model also includes governing equations for the C content of the cellular exudates of all primary producer groups (as volume concentration, eqs. 5-8) and for their dead C biomass (consistent units, eqs. 9-12). These processes were included because the presence of pesticides in ARO and an increased temperature may have important direct effects on decomposers (e.g., heterotrophic bacteria and fungi) and, consequently, on the rates of nutrient recycling within the microcosms (Nielsen, 2006). Additional governing equations describe the areal density of dead phytoplankton cells accumulating in the sediments (eq. 13, specific to each phytoplankton group), and the volume concentration of dissolved inorganic phosphorus (eq. 14). The resulting system of ODEs is:

$$ \frac{dC_{\text{live,phyto}}}{dt} = (p_{\text{phyto}} - l_{\text{phyto}} - \frac{w_{\text{live,phyto}}}{H_V})C_{\text{live,phyto}} \quad (1) $$

$$ \frac{dC_{\text{live,peri}}}{dt} = (p_{\text{peri}} - l_{\text{peri}})C_{\text{live,peri}} \quad (2) $$

$$ \frac{dC_{\text{live,epi}}}{dt} = (p_{\text{epi}} - l_{\text{epi}})C_{\text{live,epi}} \quad (3) $$

$$ \frac{dC_{\text{live,macro}}}{dt} = (p_{\text{macro}} - l_{\text{macro}})C_{\text{live,macro}} \quad (4) $$

$$ \frac{dC_{\text{exu,phyto}}}{dt} = l_{\text{exu,phyto}} C_{\text{live,phyto}} - b_{\text{exu,phyto}} C_{\text{exu,phyto}} \quad (5) $$

$$ \frac{dC_{\text{exu,peri}}}{dt} = l_{\text{exu,peri}} C_{\text{live,peri}} \frac{A_{\text{peri}}}{V} - b_{\text{exu,peri}} C_{\text{exu,peri}} \quad (6) $$

$$ \frac{dC_{\text{exu,epi}}}{dt} = l_{\text{exu,epi}} C_{\text{live,epi}} \frac{A_{\text{epi}}}{V} - b_{\text{exu,epi}} C_{\text{exu,epi}} \quad (7) $$

$$ \frac{dC_{\text{exu,macro}}}{dt} = l_{\text{exu,macro}} C_{\text{live,macro}} \frac{1}{V} - b_{\text{exu,macro}} C_{\text{exu,macro}} \quad (8) $$

$$ \frac{dC_{\text{dead,phyto}}}{dt} = l_{\text{d,phyto}} C_{\text{live,phyto}} - \left( \frac{w_{\text{dead,phyto}}}{H_V} + b_{\text{dead,phyto}} \right) C_{\text{dead,phyto}} \quad (9) $$

$$ \frac{dC_{\text{dead,peri}}}{dt} = l_{\text{d,peri}} C_{\text{live,peri}} - b_{\text{dead,peri}} C_{\text{dead,peri}} \quad (10) $$

$$ \frac{dC_{\text{dead,epi}}}{dt} = l_{\text{d,epi}} C_{\text{live,epi}} - b_{\text{dead,epi}} C_{\text{dead,epi}} \quad (11) $$

$$ \frac{dC_{\text{dead,macro}}}{dt} = l_{\text{d,macro}} C_{\text{live,macro}} - b_{\text{dead,macro}} C_{\text{dead,macro}} \quad (12) $$

$$ \frac{dC_{\text{dead,sed,phyto}}}{dt} = w_{\text{live,phyto}} C_{\text{live,phyto}} + w_{\text{dead,phyto}} C_{\text{dead,phyto}} - b_{\text{dead,phyto}} C_{\text{dead,sed,phyto}} \quad (13) $$

$$ \frac{dP_{\text{d}}}{dt} = S_{\text{phyto}} + S_{\text{peri}} + S_{\text{epi}} + S_{\text{macro}} \quad (14) $$

where the parameters $p_i=\text{phyto, peri,epi,macro}$ and $l_i=\text{phyto,peri,epi,macro}$ (eqs. 1-4) are variable total unit gain and unit loss rates of C biomass of each group; $w_{\text{live,phyto}}$ (eqs. 1,13) and $w_{\text{dead,phyto}}$ (eqs. 9,13) are constant sinking velocities of live and dead phytoplankton cells, respectively; $H_V$ is the constant height of water in the microcosm of constant water volume V; the $l_{\text{exu,i=phyto,peri,epi,macro}}$ (eqs. 5-8) are the variable unit background exudation rates of each group; the $b_{\text{exu,i=phyto,peri,epi,macro}}$ (eqs. 5-8) are the variable unit biodegradation rates of exudates; the $b_{\text{dead,i=phyto,peri,epi,macro}}$ (eqs. 9-13) are the variable unit biodegradation

rates of dead biomass; $A_{\text{peri}}$ is the constant area of the microcosm surface where wall periphyton can grow; and $A_{\text{epi}}$ the variable area of macrophyte surfaces where epiphyton can grow, formulated as a function of macrophyte biomass. The time coordinate and all state variables (governing eqs. 1-14) are listed and described in Table 1. Primary auxiliary variables, i.e., those appearing in the governing equations of the model, are listed and described in Table 2. The formulations of all secondary auxiliary variables appearing in the equations of primary auxiliary variables are presented as part of the supplementary information (SI), as well as the chosen values of all model parameters (Tables S1-S6) and the initial values of all state variables (Table S7).

Note that, in the model, periphyton and epiphyton are treated as separate primary producer groups. While periphyton grows on the vertical inner glass surface of the microcosm, which is of constant area, epiphyton grows on macrophyte surfaces that change over time as macrophytes develop. We implemented a space limitation factor in the auxiliary equation for the growth rates of periphyton and epiphyton (eqs. S19 to S22). Additionally, the growth rates of all primary producers are also dependent on nutrient and light availability, and are affected by temperature. On this note, in general, all modelled processes are accelerated under warmer conditions, as described in the SI. The effect of the herbicide was modeled based on log-logistic dose-response curves that are common in toxicokinetic studies (e.g., Copin and Chèvre, 2015). These curves were made specific to each group of primary producers, as described in the SI (Figure S1).

## 2.2 Microcosm experiment

### 2.2.1 Microcosm setup

To mimic shallow aquatic ecosystems, we used microcosms based on the OECD Guideline 239 *Water-Sediment Myriophyllum Spicatum Toxicity Test* (Allen et al., 2021; Polst et al., 2022; Vijayaraj et al., 2022). The microcosms consisted of glass vases (diameter: 25 cm, height: 40 cm, manufacturer: Sandra Rich, Germany) containing a glass bowl (diameter: 14 cm, height: 8 cm, manufacturer: Sandra Rich, Germany) filled with 20% Kaolin (Imerys, France), 5% peat (<1 mm), 1% nettle powder, 74% quartz sand (grain size fraction<0.2 mm, manufacturer: Schicker Mineral GmbH, Germany) and a 2 cm layer of quartz sand (<0.2 mm) on top to prevent resuspension of the sediment.

Three submerged macrophyte species that are common in temperate eutrophic freshwater ecosystems (Hilt et al., 2018) were collected from nearby eutrophic water bodies: *Myriophyllum spicatum* and *Elodea nuttallii* from Lake Müggelsee (Germany) and *Potamogeton perfoliatus* from River Spree (Germany). Two apical shoots of 8 cm (*M. spicatum*, *E. nuttallii*) or 10 cm (*P. perfoliatus*) were planted into the sediments of each microcosm.

All microcosms were filled with 8 L of Volvic® mineral water (Danone Waters Deutschland GmbH, Germany) to ensure homogeneity and low initial nutrient concentrations (NO₃⁻: 7.3 mg L⁻¹). After one week of initial adjustment of the system, the following species of photoautotrophic microorganisms were added, that are typical of
<table><caption>TABLE 1 Description, units of measurement of the time coordinate and state variables, and reference to the governing equations.</caption><thead><tr><th>Variable</th><th>Description</th><th>Units</th><th>Formulation</th></tr></thead><tbody><tr><td colspan="4">Time coordinate</td></tr><tr><td>t</td><td>Time coordinate, with origin at the start of the simulated period and positive direction forward</td><td>s</td><td>–</td></tr><tr><td colspan="4">State variables</td></tr><tr><td>C<sub>live,phyto</sub></td><td>Volume concentration of live carbon biomass of phytoplankton</td><td>mg C m<sup>-3</sup></td><td>Eq. 1</td></tr><tr><td>C<sub>live,peri</sub></td><td>Areal density of live carbon biomass of periphyton</td><td>mg C m<sup>-2</sup></td><td>Eq. 2</td></tr><tr><td>C<sub>live,epi</sub></td><td>Areal density of live carbon biomass of epiphyton</td><td>mg C m<sup>-2</sup></td><td>Eq. 3</td></tr><tr><td>C<sub>live,macro</sub></td><td>Live carbon biomass of macrophytes</td><td>mg C</td><td>Eq. 4</td></tr><tr><td>C<sub>exu,phyto</sub></td><td>Volume concentration of carbon in phytoplankton exudates</td><td>mg C m<sup>-3</sup></td><td>Eq. 5</td></tr><tr><td>C<sub>exu,peri</sub></td><td>Volume concentration of carbon in periphyton exudates</td><td>mg C m<sup>-3</sup></td><td>Eq. 6</td></tr><tr><td>C<sub>exu,epi</sub></td><td>Volume concentration of carbon in epiphyton exudates</td><td>mg C m<sup>-3</sup></td><td>Eq. 7</td></tr><tr><td>C<sub>exu,macro</sub></td><td>Volume concentration of carbon in macrophyte exudates</td><td>mg C m<sup>-3</sup></td><td>Eq. 8</td></tr><tr><td>C<sub>dead,phyto</sub></td><td>Volume concentration of dead carbon biomass of phytoplankton</td><td>mg C m<sup>-3</sup></td><td>Eq. 9</td></tr><tr><td>C<sub>dead,peri</sub></td><td>Areal density of dead carbon biomass of periphyton</td><td>mg C m<sup>-2</sup></td><td>Eq. 10</td></tr><tr><td>C<sub>dead,epi</sub></td><td>Areal density of dead carbon biomass of periphyton</td><td>mg C m<sup>-2</sup></td><td>Eq. 11</td></tr><tr><td>C<sub>dead,macro</sub></td><td>Dead carbon biomass of macrophytes</td><td>mg C</td><td>Eq. 12</td></tr><tr><td>C<sub>dead,sed,phyto</sub></td><td>Areal density of dead carbon biomass of phytoplankton that is stored in the sediments</td><td>mg C m<sup>-2</sup></td><td>Eq. 13</td></tr><tr><td>P<sub>d</sub></td><td>Volume concentration of dissolved inorganic phosphorus</td><td>mg P m<sup>-3</sup></td><td>Eq. 14</td></tr></tbody></table>

shallow freshwater ecosystems: four planktonic species – including two cyanobacteria (*Chroococcus minutus* and *Anabaena* sp. PCC7210) and two green algae (*Scenedesmus obliquus* and *Desmodesmus subspicatus*) – and four benthic species – including a diatom (*Nitzschia palea*), two filamentous green algae (*Uronema* sp. and *Oedogonium* sp.) and a cyanobacterium (*Komvophorum* sp.). The species were selected and

grown in Volvic® water prior to inoculation, and then combined for an inoculum with a total biovolume of 1.25×10⁹ µm³. Three plastic strips (30×2.5 cm) were fixed in the sediment bowl and attached to the inner glass vase, providing a surface for periphyton growth (150 cm² each). An aeration system was constructed to ensure mixing of the water within each microcosm using air pumps.

<table><caption>TABLE 2 Description, units of measurement of primary auxiliary variables (appearing in the governing equations), and reference to their formulation.</caption><thead><tr><th>Variable</th><th>Description</th><th>Units</th><th>Formulation</th></tr></thead><tbody><tr><td>p<sub>phyto</sub></td><td>Fractional growth rate of carbon biomass of phytoplankton groups</td><td>s<sup>-1</sup></td><td>Eq. S15</td></tr><tr><td>p<sub>peri</sub></td><td>Fractional growth rate of carbon biomass of periphyton</td><td>s<sup>-1</sup></td><td>Eq. S19</td></tr><tr><td>p<sub>epi</sub></td><td>Fractional growth rate of carbon biomass of epiphyton</td><td>s<sup>-1</sup></td><td>Eq. S20</td></tr><tr><td>A<sub>epi</sub></td><td>Area of macrophyte surfaces where epiphyton can grow</td><td>m<sup>2</sup></td><td>Eq. S24</td></tr><tr><td>p<sub>macro</sub></td><td>Fractional growth rate of carbon biomass of macrophytes</td><td>s<sup>-1</sup></td><td>Eq. S25</td></tr><tr><td>I<sub>i=phyto,peri,epi,macro</sub></td><td>Fractional loss rate of carbon biomass of primary producers</td><td>s<sup>-1</sup></td><td>Eq. S28</td></tr><tr><td>I<sub>d,i=phyto,peri,epi,macro</sub></td><td>Fractional death rate of carbon biomass of primary producers</td><td>s<sup>-1</sup></td><td>Eq. S29</td></tr><tr><td>I<sub>exu,i=phyto,peri,epi,macro</sub></td><td>Fractional exudation rate of carbon biomass of primary producers</td><td>s<sup>-1</sup></td><td>Eq. S30</td></tr><tr><td>b<sub>dead,i=phyto,peri,epi,macro</sub></td><td>Fractional biodegradation rate of dead carbon biomass of primary producers</td><td>s<sup>-1</sup></td><td>Eq. S34</td></tr><tr><td>b<sub>exu,i=phyto,peri,epi,macro</sub></td><td>Fractional biodegradation rate of cellular exudates of primary producers</td><td>s<sup>-1</sup></td><td>Eq. S35</td></tr><tr><td>S<sub>phyto</sub></td><td>Source/sink of phosphorus, phytoplankton</td><td>mg P m<sup>-3</sup></td><td>Eq. S36</td></tr><tr><td>S<sub>peri</sub></td><td>Source/sink of phosphorus, periphyton</td><td>mg P m<sup>-3</sup></td><td>Eq. S37</td></tr><tr><td>S<sub>epi</sub></td><td>Source/sink of phosphorus, epiphyton</td><td>mg P m<sup>-3</sup></td><td>Eq. S38</td></tr><tr><td>S<sub>macro</sub></td><td>Source/sink of phosphorus, macrophytes</td><td>mg P m<sup>-3</sup></td><td>Eq. S39</td></tr></tbody></table>
The microcosms were placed inside two climate chambers, each containing a reference temperature sensor in one of the microcosms. Heating and cooling of the microcosms happened *via* air temperature changes. A 16h:8h light:dark cycle was applied using luminescent light with a mean of 77.2 ± 9.9 µmol m⁻² s⁻¹ measured at the water surface.

## 2.2.2 ARO composition

The artificial mixture of pesticides representing a characteristic agricultural run-off (ARO) was selected as in our previous studies (e.g., Allen et al., 2021), including three organic pesticides – terbuthylazine (selective chloro-s-triazine herbicide, PSII inhibitor, CAS number 5915-41-3), pirimicarb (fast-acting selective carbamate insecticide, AChE inhibitor, CAS number 23103-98-2) and tebuconazole (triazole fungicide, demethylation inhibitor, CAS number 107534-96-3) – copper (inorganic pesticide, as CuSO₄) and nitrate (fertilizer, as KNO₃). The initial concentrations of pesticides were based on dose-response assays and subsequent EC20 calculations (as described in Allen et al., 2021). An initial concentration of 9 mg L⁻¹ NO₃-N was selected as representative for nitrogen concentrations in small lakes in agricultural catchments during spring and summer (James et al., 2005; Xu et al., 2014). For the application in the microcosms, the organic pesticides were diluted in dimethyl sulfoxide (DMSO, Sigma-Aldrich, 67-68-5, final concentration in the microcosm 0.1%). CuSO₄ and KNO₃ were diluted in ultrapure water.

## 2.2.3 Two-factor factorial dose-response design

The experiment was conducted in two climate chambers (22°C, 26°C) and in each, we followed a dose-response setup for ARO with one control and five concentration levels (ARO 1x, 2x, 4x, 8x, 16x; enrichment factor of 2). The ARO mixture and concentrations were chosen based on former work by Allen et al. (2021). The ARO concentration used in their study was set as intermediate treatment concentration in our experiment (ARO 4x) with two higher and two lower concentrations for the gradient. Each treatment had five to seven replicates.

## 2.2.4 Time schedule

The bowls filled with sediment were prepared and pre-wetted with Volvic® mineral water and stored in dark conditions at 22°C for two days before macrophytes were planted and microcosms were filled with Volvic® mineral water. Afterwards, macrophytes were given two weeks to adjust to the conditions. The inoculum of microscopic primary producers and plastic strips were added one week before the start of the experiment. At the start of the experiment, treatments were applied by adding the ARO mix to the microcosms and increasing the temperature in one of the climate chambers to 26°C. Hereafter, nutrients in the form of a KNO₃ and KH₂PO₄ mixture as in Allen et al. (2021) were added twice a week, simulating repeated nutrient loading to compensate for a fast nutrient uptake by wall periphyton, and thus ensure sufficient nutrient availability to sustain further primary producer growth. On days 16, 17 and 18, we sampled periphyton and phytoplankton, and took water samples for chemical analysis (day

18). On day 19, to end the experiment, the lights were turned off and sampling of the macrophytes took place until day 23 (Figure S2).

## 2.2.5 Sampling methods

### 2.2.5.1 Biomass sampling of primary producers

Periphyton was brushed off from the plastic strips using a toothbrush and suspended in Volvic® water. This periphyton-suspension was then filtered with pre-weighted glass fiber filters (0.7 µm), dried at 60°C and weighted to determine the periphyton dry weight. To determine phytoplankton biomass, water samples were filtered using pre-weighted glass fiber filters (0.7 µm), dried at 60°C and weighted. Macrophytes were removed from their microcosm for the final sampling and the aboveground part of each macrophyte specimen was separately packed in paper bags, dried at 60°C for two days and then weighted to determine the dry weight of each stem.

### 2.2.5.2 Water sampling for pesticide concentrations

Two hours and 18 days after application of the treatments, 4 mL-water subsamples of three microcosms per ARO treatment were taken and frozen at -20°C. Later, the samples were filtered (0.2 µm) and analyzed for their concentration of the three pesticides used in the ARO mixture. Measurement of the pesticides was conducted using a UltiMate3000 HPLC System combined with an LTQ-OrbiTrap XL (Thermo Scientific, USA). These samples were then analyzed with an UltiMate3000 HPLC System (column: Phenomenex, Art.-No. 00B-4462-Y0) and an attached LTQ Orbitrap XL (Thermo Scientific) operated in positive ionization mode.

## 2.2.6 Statistical analyses

Non-parametric statistical tests were used to test for differences between treatments related to the ARO application or the temperature. These analyses were conducted using the software R (R Core Team, 2020; v4.0.0). For each of the three primary producer groups differences in biomass were tested *via* the Kruskal-Wallis test followed by Dunn's *post hoc* test without further correction due to the low number of treatments.

Further, biomass of each primary producer group was extrapolated to the whole biomass per microcosm, including periphyton growing on the inner walls of the microcosm. Based on this total microcosm primary producer biomass, the effect sizes (as Hedges' g) were calculated using the *esc* package (Lüdecke, 2019; v0.5.1). Using the mean biomass of each primary producer group as well as the total biomass of all primary producer groups together, pie charts were created for simplified presentation of the share of each primary producer group within each treatment.

## 2.3 Model refinement and assessment against experimental results

The information we progressively gained over the course of the experiment and after experimental results became available allowed for further refinement of the model to account for some processes
like the decay of toxicant concentrations over time (revealed by water sampling results), or to better represent processes like nutrient uptake by wall periphyton, nutrient recycling due to biodegradation, and trait-based phytoplankton community adaptation by splitting the phytoplankton compartment into two groups. Group α was conceived as a fast-growing phytoplankton group that was highly sensitive to the herbicide (i.e., low EC50 value), whereas group β was made up of slow-growing phytoplankton that were resistant to the herbicide (i.e., a much higher EC50 value than group α).

Moreover, the log-logistic dose-response curves were made variable in time to account for the potential development of tolerance (adaptation) to the herbicide and/or acclimation of the organisms to better cope with the combined stressors. We assumed these adaptation and acclimation processes to be directly linked to the generation time of the organisms, i.e., faster for microscopic primary producers. This was achieved by implementing dynamic parameters for the dose-response curves, as described in the SI. The potential of the refined model to simulate community adaptation, tolerance development and organism acclimation led us to design and test the model under different scenarios of herbicide sensitivity, temperature dependence of primary producer response and adaptation and/or acclimation.

### 2.3.1 Scenarios of differential sensitivities, temperature dependence, adaptation and/or acclimation

We designed eight different scenarios to test our hypotheses and to mechanistically understand experimental results. These scenarios related to: a) whether microscopic primary producers (phytoplankton group α, periphyton and epiphyton) were equally, more or less sensitive to the herbicide than macrophytes; b) whether these sensitivities were influenced by temperature or

not; and c) whether sensitivities decreased over time or not as a result of the development of tolerance to the herbicide by adaptation or acclimation of microscopic primary producers to environmental conditions (Table 3).

In scenarios A1 and A2, all primary producer groups were equally sensitive to the herbicide except for phytoplankton group β (conceived as slow-growing and resistant to the herbicide). In scenarios B1 and B2, macrophytes were more sensitive to the herbicide than all other primary producer groups. In scenarios C1 and C2, microscopic primary producers were more sensitive to the herbicide than macrophytes, with the exception of phytoplankton group β, which, again, was more resistant than all other groups.

In all aforementioned scenarios (A1-2, B1-2 and C1-2), we assumed that no development of tolerance to the herbicide by acclimation or adaptation to environmental conditions (Lips et al., 2022) occurred over the course of the exposure, i.e., sensitivities do not decrease over time. In scenarios D1 and D2, however, we account for this by implementing decreasing herbicide sensitivities over the course of the exposure for all microscopic primary producers. In these two scenarios, phytoplankton group α, periphyton and epiphyton start off being more sensitive to the herbicide than macrophytes, but end up becoming less sensitive than macrophytes by the end of the simulated period as a result of tolerance development. Here too, phytoplankton group β starts off being less sensitive than all other groups, and becomes increasingly tolerant over time.

In all scenarios indexed with “1”, we assumed that temperature had no effect at all on herbicide sensitivities, whereas in all scenarios indexed with “2”, sensitivities to the herbicide were assumed to be higher at the lower temperature (22°C) than at the higher temperature (26°C), following empirical evidence by Tasmin et al., 2013. Moreover, this is consistent with the fact that metabolic processes are in general accelerated at higher

<table><caption>TABLE 3 Scenarios of differential sensitivities, temperature dependence, adaptation and/or acclimation.</caption><thead><tr><td>Scenario</td><td>Herbicide sensitivity by group</td><td>Effect of temperature on herbicide sensitivity</td></tr></thead><tbody><tr><td>A1</td><td>Fast-growing phytoplankton (group α), wall periphyton, epiphyton and macrophytes are equally sensitive to the herbicide.</td><td>NO</td></tr><tr><td>A2</td><td>Slow-growing phytoplankton (group β) are less sensitive to the herbicide than all other groups. Sensitivities do not decrease over time for any group (no adaptation or acclimation).</td><td>Sensitivities are higher at the lower temperature (22°C)</td></tr><tr><td>B1</td><td>Macrophytes are more sensitive to the herbicide than all microscopic primary producers.</td><td>NO</td></tr><tr><td>B2</td><td>Slow-growing phytoplankton (group β) are less sensitive to the herbicide than all other groups. Sensitivities do not decrease over time (no adaptation or acclimation).</td><td>Sensitivities are higher at the lower temperature (22°C)</td></tr><tr><td>C1</td><td>Fast-growing phytoplankton (group α), wall periphyton and epiphyton are more sensitive to the herbicide than macrophytes.</td><td>NO</td></tr><tr><td>C2</td><td>Slow-growing phytoplankton (group β) are less sensitive to the herbicide than all other groups. Sensitivities do not decrease over time (no adaptation or acclimation).</td><td>Sensitivities are higher at the lower temperature (22°C)</td></tr><tr><td>D1</td><td>Initially, fast-growing phytoplankton (group α), periphyton and epiphyton are more sensitive to the herbicide than macrophytes. These groups become less sensitive to the herbicide over time to become less sensitive to macrophytes by the end of the exposure.</td><td>NO</td></tr><tr><td>D2</td><td>From the beginning of the exposure, low-growing phytoplankton (group β) are less sensitive to the herbicide than all other groups and become increasingly tolerant to it over time.</td><td>Sensitivities are higher at the lower temperature (22°C)</td></tr></tbody></table>
temperatures (Clarke and Fraser, 2004; Pedrosa Gomes and Juneau, 2017).

To further account for possible community adaptation and succession of phytoplankton to the applied stressors related to an increasing dominance of the more tolerant species (Blanck, 2002; Tlili et al., 2016; Lips et al., 2022), we ran three sets of simulations under all eight scenarios. Set 1: a single phytoplankton group that is fast-growing and highly sensitive to the herbicide, hereafter referred to as group α; Set 2: a single phytoplankton group that is slow-growing but much less sensitive to the herbicide, hereafter referred to as group β; Set 3: a mixed community of both phytoplankton groups α and β, as previously described.

### 2.3.2 Assessment of model simulation results

We assessed the goodness of fit between model results and the observations from the microcosm experiment following a pattern-oriented model validation approach that took the ARO treatment as predictor of final biomass values for each of the following three primary producer groups: phytoplankton, periphyton (wall growth in the model, periphyton strip growth in the experiment) and macrophytes. The resulting patterns were compared qualitatively in terms of the

similarity between simulated and observed response curves and quantitatively on the basis of the correlation coefficients calculated for all scenarios between model results and empirical observations.

## 3 Results

### 3.1 Microcosm experiment

#### 3.1.1 Biomass of different primary producer groups

At 22°C, all ARO treatments showed significantly lower macrophyte biomass than the controls (p-value < 0.05), but did not significantly differ among the five ARO treatments. At 26°C, macrophyte biomass across all ARO concentrations was not significantly different in relation to the controls or the respective ARO concentrations at 22°C (Figures 1, 2 and Table 4). At 22°C, phytoplankton biomass showed a significant increase compared with the controls at the second tested ARO concentration level (ARO 2x) and higher (apart from ARO 8x). At 26°C, phytoplankton biomass was significantly higher at ARO 4x and ARO 16x. No significant differences

<figure><img src="image_3.png" /><figcaption>FIGURE 1
Effect sizes (Hedges' g) and its standard deviation for the response of different primary producer groups (blue = macrophytes, green = phytoplankton, orange = periphyton) and ARO treatment levels at 22°C (top row) and 26°C (bottom row).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2
Response of different primary producer groups to multiple stressors by increasing concentrations of pesticides and nutrients in agricultural runoff (ARO) and warming (+4°C) observed in experimental microcosms (left) and simulated (right) for a mixed phytoplankton community (set 3) under the most plausible scenario (D2) for warm (26°C, red) and ambient (22°C, blue) temperature treatments. Asterisks indicate significant differences with respect to the control.</figcaption></figure>

between the respective ARO concentrations at the two temperatures were detected (Figures 1, 2). Periphyton biomass was not significantly affected by ARO neither at 22°C nor at 26°C, likely due to the high variation among replicates in the controls (Figure 2). Effect sizes for periphyton showed negative effects at the two highest tested ARO levels for both temperatures (Figure 1).

The average accumulated microcosm biomass for all primary producers decreased in the ARO treatments compared with the control at 22°C. At 26°C, the accumulated microcosm biomass of the ARO treatments was higher than their control at 26°C and the respective ARO treatments at 22°C. (Figure 3). In the control treatment, macrophytes dominated with periphyton having the

<table><caption>TABLE 4 Nominal concentrations of all ARO constituents and measured concentrations 2 hours after its application (n=3, ARO 8: n=1,< d.l. = below detection limit). Concentrations in µg L⁻¹.</caption><thead><tr><th>Treatment</th><th>Terbuthylazine</th><th>Pirimicarb</th><th>Tebuconazole</th><th>Copper</th><th>NO<sub>3</sub>-N</th></tr></thead><tbody><tr><td>Control</td><td>0 (&lt; d.l.)</td><td>0 (&lt; d.l.)</td><td>0 (&lt; d.l.)</td><td>0 (&lt; d.l.)</td><td>0 (&lt; d.l.)</td></tr><tr><td>ARO 1x</td><td>0.75 | 0.64 ± 0.01</td><td>3.75 | 3.82 ± 0.22</td><td>22.5 | 21.11 ± 1.21</td><td>10.5</td><td>2250</td></tr><tr><td>ARO 2x</td><td>1.5 | 1.28 ± 0.06</td><td>7.5 | 7.24 ± 0.14</td><td>45 | 38.07 ± 4.00</td><td>21</td><td>4500</td></tr><tr><td>ARO 4x</td><td>3 | 2.74 ± 0.06</td><td>15 | 15.56 ± 0.62</td><td>90 | 83.46 ± 6.56</td><td>42</td><td>9000</td></tr><tr><td>ARO 8x</td><td>6 | 5.70 ± 0.00</td><td>30 | 31.07 ± 0.00</td><td>180 | 171.80 ± 0.00</td><td>84</td><td>18000</td></tr><tr><td>ARO 16x</td><td>12 | 11.6±+0.59</td><td>60 | 63.64 ± 3.02</td><td>360 | 323.03 ± 23.06</td><td>168</td><td>36000</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIGURE 3</figcaption></figure>

Biomass proportions of the three primary producer groups upscaled to whole microcosm biomass. Size of the pies indicate relative total biomass in comparison with the controls (ARO 0) at 22°C.

second highest share. The share of both macrophytes and periphyton decreased with ARO treatments due to increasing phytoplankton abundance, with the lowest share of macrophytes on the accumulated microcosm biomass at ARO 4x (Figure 3). The share of phytoplankton increased at 26°C.

## 3.1.2 Change in pesticide concentrations

After two hours of exposure, the pesticide concentrations were approximately 10% lower than the nominal concentration (Table 4). One microcosm showed higher concentrations than planned (and wrong ratios between the three pesticides) and was, thereafter, excluded from further analysis. At the end of the experiment, between 55-90% of the initial pesticide concentrations were still detected (Table S8). No temperature-related differences in the final pesticide concentrations were found.

## 3.2 Process-based microcosm model

Our custom-tailored process-based model showed that, in all simulation sets, periphyton and macrophyte carbon biomass values per unit volume (mg C m⁻³) decreased along the axis of increasing ARO concentrations at both 22°C and 26°C. Whether this decrease was more or less pronounced, and the pattern itself of the response ultimately depended on how well phytoplankton developed, with a stronger decrease in benthic primary producers when phytoplankton performed better. For simulation set 1, a community of only fast-growing phytoplankton that are highly sensitive to the herbicide (group α only, Figure 4), none of the eight scenarios resulted in phytoplankton response patterns similar to those observed empirically, i.e., from the microcosm experiment (Figures 1, 2). Empirical and simulated phytoplankton response patterns to ARO were more comparable for simulation sets 2 (group β only, Figure 5) and 3 (mixed community of

phytoplankton groups α and β, Figure 6), but only under warm conditions. Correlation coefficients calculated for the 3 sets also showed that the best fits across all primary producer groups were obtained for sets 2 and 3 when herbicide sensitivities were assumed to be temperature dependent and to decrease over time as a result of tolerance development, i.e., scenario D2 (Table 5).

# 4 Discussion

To better understand the complex nature of primary producer interactions and responses to multiple stressors, we combined indoor microcosm experiments with a custom-tailored process-based model. We confirmed hypotheses 1) that ARO differentially affects primary producer groups increasing the likelihood of phytoplankton dominance in shallow aquatic ecosystems and 3) that custom-tailored process-based models can support mechanistic understanding of experimental results through scenario comparison. Hypothesis 2, a facilitation of shifts to phytoplankton dominance by warming, was only partially confirmed. Here we discuss the respective implications of our findings.

## 4.1 Agricultural run-off differentially affects primary producers and their interactions

Our results demonstrate that ARO has a differential effect on the tested primary producer groups. While we did not observe a full dominance of either primary producer group at the end of the microcosm experiment, an increasing proportion of phytoplankton and decreasing proportion of macrophytes in the total microcosm biomass clearly indicates an increasing risk for the system to shift to
<figure><img src="image_6.png" /><figcaption>FIGURE 4
Response of primary producer groups under the eight scenarios (A1-D2) for the case of a single, fast-growing phytoplankton group that is highly sensitive to the herbicide (simulation set 1, group α only) under ambient (22°C, blue) and warm (26°C, red) temperature conditions.</figcaption></figure>

full phytoplankton dominance. The overall non-linear increase of phytoplankton biomass with increasing ARO concentrations was most likely caused by nitrate (Scheffer et al., 1993; Olsen et al., 2015). Similar microcosm studies testing the same ARO mixture also found an increase in phytoplankton biomass with accompanied decrease in macrophyte biomass when nitrate and pesticides co-occurred (Polst et al., 2022; Vijayaraj et al., 2022). The increase in phytoplankton biomass may have contributed to the observed decline in macrophyte biomass by shading. Regime shifts with macrophyte loss due to nutrient loading are usually caused by increased shading by periphyton and phytoplankton (Phillips et al., 2016). Yet, at 22°C, a decline in macrophyte biomass was already observed at the lowest ARO concentration, which did not increase phytoplankton biomass, indicating a higher sensitivity of the tested submerged macrophyte species to ARO compared with phytoplankton. This contrasts the findings of Giddings et al. (2013) who reported a higher sensitivity of phytoplankton species compared with macrophytes for 4 out of 5 photosystem-II-inhibiting herbicides in a meta-analysis based on standardized single-species tests. In longer lasting experiments allowing for multiple generation cycles, microalgae community composition can change *via* acclimation, adaptation and selection of more tolerant species leading to a pollution-induced community

tolerance (Blanck, 2002; Lorente et al., 2015; Lips et al., 2022). Such processes could explain the higher sensitivity and the lower tolerance of macrophytes to ARO in comparison to phytoplankton and periphyton in our experiment. Our results thus suggest that differential sensitivities of primary producers to pesticides can modulate their competition for light and facilitate the loss of macrophytes, eventually leading to regime shifts in shallow lentic aquatic ecosystems.

## 4.2 Warming affects the risk of an ARO-induced regime shift to phytoplankton dominance

In our experiment, ARO at both tested water temperatures clearly changed the contributions of different primary producers to the total biomass towards a higher proportion of phytoplankton indicating the risk of an ARO-induced shift to full phytoplankton dominance. However, absolute macrophyte biomass did not respond to ARO in the warming (+4°C) treatment, similar to findings by Allen et al. (2021). Possibly, higher temperatures led to an enhanced macrophyte growth partially, compensating negative pesticide effects on macrophytes. Higher growth rates of
<figure><img src="image_7.png" /><figcaption>FIGURE 5
Response of primary producer groups under the (A1- D2) for the case of a single, slow-growing phytoplankton group that is highly resistant to the herbicide (simulation set 2, group *β* only) under ambient (22°C, blue) and warm (26°C, red) temperature conditions.</figcaption></figure>

macrophytes at higher temperatures are known from other studies (e.g., Zhang et al., 2019; Zhang et al., 2022). Moreover, Vijayaraj et al. (2022) and Polst et al. (2022) described a stronger decrease in macrophyte biomass and thus lower threshold levels for ARO-induced shifts to phytoplankton dominance in +4°C treatments. These results may be caused by a lower sensitivity of phytoplankton towards herbicides at elevated temperatures (Chalifour and Juneau, 2011; Tasmin et al., 2013). Periphyton also seems less sensitive to herbicides at elevated temperatures (Larras et al., 2013), but we could not confirm this neither in our experiment nor in related experiments using the same ARO mixture (Allen et al., 2021; Polst et al., 2022). Another recent study also found no effect of combined warming, nutrients and pesticides on periphyton and suggested that treatment effects compensated each other (Zhang et al., 2022). We conclude that warming facilitates macrophyte and phytoplankton growth, which can modulate macrophyte response thresholds to ARO. However, the increasing proportion of phytoplankton in the primary producer community under combined warming and ARO exposure indicates an overall increased risk of macrophyte loss, also supported by findings in Vijayaraj et al. (2022) and Polst et al. (2022).

## 4.3 Process-based models support mechanistic understanding of experimental results

Custom-tailored models of appropriate complexity, i.e., including only the processes that are strictly necessary to reproduce empirical observations, are helpful to shed light on the critical mechanisms behind ecosystem-wide responses (Vasconcelos et al., 2016; López Moreira et al., 2021). This is particularly true when addressing the effects of multiple stressors, for which the knowledge base is currently just starting to develop. Moreover, the nearly infinite amount of combinations of stressor levels limits the potential contribution of purely laboratory-based studies, e.g., microcosm experiments, where decisions are made based on several factors including space, materials and budget available, viability and phenology of model organisms, time constraints and the human effort that is needed. Considering that readily available models would need substantial modifications, recalibration and revalidation before use and may still remain unnecessarily complex for exploratory research, we suggest to calibrate and validate newly developed, custom-made simple
<figure><img src="image_8.png" /><figcaption>FIGURE 6</figcaption></figure>

Response of primary producer groups under the eight scenarios eight scenarios (A1- D2) for the case of a mixed phytoplankton community (simulation set 3, groups α and β) under ambient (22°C, blue) and warm (26°C, red) temperature conditions.

models with experimental data and to test scenarios which allow
subsequent model refinements and mechanistic interpretation of
model results.

Our pattern-oriented modeling effort allowed for the
reconstruction of primary producer temporal dynamics at
different levels of approximation to the empirical observations. It
also made it possible to assess the trajectories of all state variables
(exemplified in Figures S3, S4) and the most likely limiting
conditions to primary producer growth and their change over
time (exemplified for epiphyton in Figure S5). Because the
empirical data was mostly limited to end-point values, however,
model calibration and scenario selection was limited to final
biomass results. Correlation coefficients allowed for the
quantitative identification and selection of the most plausible
scenario that best fitted empirical observations.

4.3.1 Phytoplankton biomass

Assuming a community of only fast-growing phytoplankton that
are highly sensitive to the herbicide (set 1: group α only), correlation
coefficients for phytoplankton biomass where mostly negative
(Table 5), indicating that the model was clearly not able to
reproduce the pattern observed in the experiment, where final
biomass values were generally higher for the highest ARO

concentrations (Figures 2, 4). This complex phytoplankton
response was better captured under the assumption of a single,
slow-growing but relatively resistant phytoplankton group (set 2:
group *β* only) (Figure 5) and under the assumption of a mixed
phytoplankton community (set 3: both groups) (Figure 6), for which
we obtained several positive albeit not statistically significant
correlation coefficients (Table 5). However, this was generally true
only under warm conditions (26°C). The exceptions to this were
scenarios C1 and D1, for which correlation coefficients for
phytoplankton were negative under both temperature treatments
(both sets, although marginally less negative for set 3, mixed
phytoplankton community).

4.3.2 Periphyton and macrophyte biomass

For periphyton, statistically significant positive correlations were
obtained for all three sets under all eight scenarios at the higher
temperature of 26°C (Table 4). For the cool temperature treatment (22°
C), these positive correlations were statistically significant only under
scenarios B2, C1-2 and D1-2. Macrophyte responses to ARO were also
generally positively correlated across all scenarios, sets and temperature
treatments, but statistically significant only under scenarios B1-2 and
D1-2 at 22°C, and under scenarios B2, C1 and D1-2 at 26°C (Table 4).
However, the only scenario where the simulated macrophyte response
<table><caption>TABLE 5 Correlation coefficients between simulated and observed patterns under eight scenarios and three cases of phytoplankton community composition.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="2">Phytoplankton</th><th colspan="2">Periphyton</th><th colspan="2">Macrophytes</th></tr><tr><th>Cool</th><th>Warm</th><th>Cool</th><th>Warm</th><th>Cool</th><th>Warm</th></tr></thead><tbody><tr><td colspan="7">Set 1: Phytoplankton group α only</td></tr><tr><td>A1</td><td>-0.514</td><td>-0.532</td><td>0.582</td><td>0.876**</td><td>0.472</td><td>0.777*</td></tr><tr><td>A2</td><td>-0.514</td><td>-0.524</td><td>0.582</td><td>0.786*</td><td>0.472</td><td>0.615</td></tr><tr><td>B1</td><td>-0.502</td><td>-0.552</td><td>0.584</td><td>0.877**</td><td><strong>0.969**</strong></td><td>0.579</td></tr><tr><td>B2</td><td>-0.170</td><td>-0.549</td><td>0.748*</td><td><strong>0.927**</strong></td><td><strong>0.969**</strong></td><td>0.753*</td></tr><tr><td>C1</td><td>-0.522</td><td>-0.381</td><td>0.862**</td><td><strong>0.951**</strong></td><td>0.453</td><td>0.757*</td></tr><tr><td>C2</td><td>-0.522</td><td>-0.562</td><td>0.862**</td><td><strong>0.926**</strong></td><td>0.453</td><td>0.598</td></tr><tr><td>D1</td><td>0.013</td><td>-0.416</td><td>0.843**</td><td><strong>0.970**</strong></td><td>0.836**</td><td>0.849**</td></tr><tr><td>D2</td><td>0.013</td><td>-0.547</td><td>0.843**</td><td><strong>0.952**</strong></td><td>0.836**</td><td>0.765*</td></tr><tr><td colspan="7">Set 2: Phytoplankton group β only</td></tr><tr><td>A1</td><td>0.471</td><td>0.537</td><td>0.654</td><td><strong>0.926**</strong></td><td>0.471</td><td>0.777*</td></tr><tr><td>A2</td><td>0.471</td><td>0.506</td><td>0.654</td><td>0.798*</td><td>0.471</td><td>0.614</td></tr><tr><td>B1</td><td>0.537</td><td>0.644</td><td>0.655</td><td><strong>0.926**</strong></td><td><strong>0.969**</strong></td><td>0.579</td></tr><tr><td>B2</td><td>-0.509</td><td>0.660</td><td>0.805*</td><td><strong>0.968**</strong></td><td><strong>0.969**</strong></td><td>0.753*</td></tr><tr><td>C1</td><td>-0.522</td><td>-0.523</td><td>0.883**</td><td>0.869**</td><td>0.465</td><td>0.770*</td></tr><tr><td>C2</td><td>-0.522</td><td>0.539</td><td>0.883**</td><td><strong>0.968**</strong></td><td>0.465</td><td>0.606</td></tr><tr><td>D1</td><td>-0.491</td><td>-0.522</td><td>0.870**</td><td>0.915**</td><td>0.837**</td><td>0.849**</td></tr><tr><td>D2</td><td>-0.491</td><td>0.586</td><td>0.870**</td><td><strong>0.979**</strong></td><td>0.837**</td><td>0.766*</td></tr><tr><td colspan="7">Set 3: Mixed community of phytoplankton groups α and β</td></tr><tr><td>A1</td><td>0.477</td><td>0.518</td><td>0.590</td><td>0.878**</td><td>0.472</td><td>0.777*</td></tr><tr><td>A2</td><td>0.477</td><td>0.544</td><td>0.590</td><td>0.779*</td><td>0.472</td><td>0.615</td></tr><tr><td>B1</td><td>0.507</td><td>0.480</td><td>0.592</td><td>0.880**</td><td><strong>0.969**</strong></td><td>0.579</td></tr><tr><td>B2</td><td>-0.582</td><td>0.439</td><td>0.778*</td><td><strong>0.951**</strong></td><td><strong>0.969**</strong></td><td>0.753*</td></tr><tr><td>C1</td><td>-0.365</td><td>-0.344</td><td>0.861**</td><td>0.890**</td><td>0.459</td><td>0.767*</td></tr><tr><td>C2</td><td>-0.365</td><td>0.432</td><td>0.861**</td><td><strong>0.950**</strong></td><td>0.459</td><td>0.604</td></tr><tr><td>D1</td><td>-0.414</td><td>-0.365</td><td>0.850**</td><td><strong>0.939**</strong></td><td>0.837**</td><td>0.849**</td></tr><tr><td>D2</td><td>-0.414</td><td>0.505</td><td>0.850**</td><td><strong>0.971**</strong></td><td>0.837**</td><td>0.766*</td></tr></tbody></table>

pattern qualitatively resembled the observed one was scenario D2 (Figure 6).

### 4.3.3 Most plausible scenario

To identify the most plausible scenario, we first discarded all scenarios where final macrophyte biomass patterns were not significantly correlated under both temperature treatments. Scenarios A1-2, B1, and C1-2 were thus removed from further consideration. Subsequently, scenario D1 was also discarded, as it failed to produce positive correlations for phytoplankton biomass. Finally, considering that, between the two remaining scenarios, i.e., scenarios B2 and D2: a) phytoplankton responses were both

qualitatively and quantitatively fundamentally the same under scenarios B2 and D2; b) positive correlations for periphyton were higher and more statistically significant under scenario D2; and that, c) qualitatively, macrophytes were best simulated under scenario D2, quantitatively supported by statistically significant and highly positive correlation coefficients; we identified and selected scenario D2 as the most plausible scenario.

### 4.3.4 Comparisons among the three simulation sets

Because the model failed to produce good fits for the final phytoplankton biomass response pattern for set 1 (group α only),

<footnote>* Significant at a 90% confidence level (p-value<0.1).</footnote>

<footnote>** Significant at a 95% confidence level (p-value<0.05).</footnote>

<footnote>** (in bold) Significant at a 99% confidence level (p-value<0.01).</footnote>
with all correlation coefficients being negative at 26°C and only marginally positive for scenarios D1-2 at 22°C (Table 5), this simulation set was removed from any further consideration.

Under scenario D2, selected as the most plausible one, for simulation sets 2 (group *β* only) and 3 (mixed community of groups *α* and *β*), simulated and observed responses were positively correlated for both periphyton (lowest correlation coefficients were 0.870 and 0.850 for sets 2 and 3 at 22°C, respectively) and macrophytes (lowest correlation coefficients were 0.766 for both sets 2 and 3 at 26°C). All these correlations were statistically significant at a confidence level of either 90%, 95% or even 99% (Table 5). Phytoplankton responses, however, were only weakly positively correlated in the warming treatment (0.586 and 0.505 for sets 2 and 3 at 26°C, respectively) and even negatively correlated under ambient temperature (22°C). However, none of these correlations for phytoplankton were statistically significant.

In scenario D2, sensitivities to the herbicide of all primary producer groups were dependent on temperature (organisms were less sensitive at the higher temperature). Additionally, herbicide sensitivities of all microscopic primary producers also decreased over time as they became more tolerant to the herbicide. This resulted in a simulated pattern under scenario D2 that was positively correlated with the observed phytoplankton response for the warm temperature treatment. Under scenario D1, where sensitivities were not affected by temperature, a positive correlation could not be achieved for any of the sets at neither of the two temperature treatments.

In our study, the most plausible scenario was the one involving stress-induced phytoplankton tolerance development, including temperature-dependent organism acclimation and adaptation to ARO (scenario D2) as well as stress-induced succession in communities (simulation set 3, mixed phytoplankton community). The relevance of these processes has been shown under the concept of pollution-induced community tolerance (Blanck, 2002; Schmitt-Jansen et al., 2016; Lips et al., 2022). Another recent study reported observed changes in phytoplankton communities of small lentic waters following exposure to agricultural run-off (Wijewardene et al., 2021) but did not test for changes in tolerance to the pesticides. While designed for a specific experimental set-up, our model can be used for future hypothesis testing in experimental studies with a similar complexity (focusing on primary producer interactions), but can also be extended to include higher trophic levels.

## 4.4 Final remarks

Our microcosm experiment revealed a differential response of aquatic primary producers to the combined effects of warming and ARO, potentially leading to phytoplankton dominance, a less desirable ecosystem state. At both tested temperatures, phytoplankton was favored and increased its biomass and proportion in the primary producer community with increasing ARO concentrations due to the higher availability of nitrate and ability to adapt to pesticides. Macrophytes became less sensitive to

ARO with warming, but may still be weakened due to the increasing proportion of phytoplankton eventually shading them out. The development of a simple, process-based model allowed understanding the role of temperature dependence, tolerance development and community adaptation to combined ARO and warming, because a scenario including these processes ultimately led to the best fit between predicted and observed responses. Our results highlight the importance of considering stress-induced tolerance development, adaptation and succession when predicting mid- and long-term effects of toxicants on complex primary producer communities that include multiple species. Trait diversity may compensate for the direct negative effects of pesticides on individual species to keep ecosystems productive at a macroscopic scale, especially at higher temperatures and under nutrient-enriched conditions, where ARO may unbalance the system through indirect effects on vulnerable communities like macrophytes.

## Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author contributions

The experimental part of the study was mostly done by BP. SH and other colleagues at the IGB Berlin offered help during the final sampling activities of the experiment. The process-based model was developed and implemented by GLMM. SH, FH, MS-J, and BP all provided guidance for the development of the model. GLMM and BP contributed equally to the final manuscript. EG helped setting up the experiment, provided guidance and contributed to the discussion and editing of the manuscript as Principal Investigator of the CLIMSHIFT project, which funded this work. All authors contributed to the article and approved the submitted version.

## Funding

This work is part of the CLIMSHIFT French-German project funded by the DFG (SCHM 2555/5-1 to MS-J and HI 1380/8-1 to SH) and ANR (ANR-17-CE32-0013 to EG).

## Acknowledgments

Special thanks to Valeria Lucero Reimann Salazar, Simon Asmus and other members of IGB who helped during the microcosm experiment, and to Herwig Stibor, Joséphine Leflaive, Joey Allen and other members of CLIMSHIFT for providing valuable comments as part of their role in the project and for their help in the preparation of the microcosm experiment and the development of the process-based model.
## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated

organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1120441/full#supplementary-material

## References

Allen, J., Gross, E. M., Courcoul, C., Bouletreau, S., Compin, A., Elger, A., et al. (2021). Disentangling the direct and indirect effects of agricultural runoff on freshwater ecosystems subject to global warming: a microcosm study. *Water Res.* 190, 116713. doi: 10.1016/j.watres.2020.116713

Blanck, H. (2002). A critical review of procedures and approaches used for assessing pollution-induced community tolerance (PICT) in biotic communities. *Hum. Ecol. Risk Assessment: Int. J.* 8 (5), 1003–1034. doi: 10.1080/1080-700291905792

Chalifour, A., and Juneau, P. (2011). Temperature-dependent sensitivity of growth and photosynthesis of scenedesmus obliquus, navicula pelliculosa and two strains of microcystis aeruginosa to the herbicide atrazine. *Aquat. Toxicol.* 103 (1-2), 9–17. doi: 10.1016/j.aquatox.2011.01.016

Clarke, A., and Fraser, K. (2004). Why does metabolism scale with temperature? *Funct. Ecol.* 18 (2), 243–251. doi: 10.1111/j.0269-8463.2004.00841.x

Copin, P.-J., and Chèvre, N. (2015). Modelling the effects of pulse exposure of several PSII inhibitors on two algae. *Chemosphere* 137, 70–77. doi: 10.1016/j.chemosphere.2015.05.035

Fairchild, J. F., Ruessler, D. S., and Carlson, R. A. (1998). Comparative sensitivity of five species of macrophytes and six species of algae to atrazine, metribuzin, alachlor, and metolachlor. *Environ. Toxicol. Chem.* 17 (9), 1830–1834. doi: 10.1002/etc.5620170924

FAO (2021). "Land cover statistics. global, regional and country trends 2000-2019," in *FAOSTAT analytical brief series no. 37* (Rome: FAO).

Giddings, J. M., Arts, G., and Hommen, U. (2013). The relative sensitivity of macrophyte and algal species to herbicides and fungicides: an analysis using species sensitivity distributions. *Integrated Environ. Assess. Manage.* 9 (2), 308–318. doi: 10.1002/ieam.1387

Hilt, S., Brothers, S., Jeppesen, E., Veraart, A. J., and Kosten, S. (2017). Translating regime shifts in shallow lakes into changes in ecosystem functions and services. *BioScience* 67 (10), 928–936. doi: 10.1093/biosci/bix106

Hilt, S., Alirangues Nuñez, M. M., Bakker, E. S., Blindow, I., Davidson, T. A., Gillefalk, M., et al. (2018). Response of submerged macrophyte communities to external and internal restoration measures in north temperate shallow lakes. *Front. Plant Sci.* 9. doi: 10.3389/fpls.2018.00194

Jackson, M. C., Loewen, C. J. G., Vinebrooke, R. D., and Chimimba, C. T. (2016). Net effects of multiple stressors in freshwater ecosystems: a meta-analysis. *Global Change Biol.* 22 (1), 180–189. doi: 10.1111/gcb.13028

Jäger, C., Diehl, S., and Emans, M. (2010). Physical determinants of phytoplankton production, algal stoichiometry and vertical nutrient fluxes. *Am. Nat.* 175 (4), E91–E104. doi: 10.1086/650728

Jakeman, A., Letcher, R., and Nortonac, J. (2006). Ten iterative steps in development and evaluation of environmental models. *Environ. Model. Software* 21 (5), 602–614. doi: 10.1016/j.envsoft.2006.01.004

James, C., Fisher, J., Russell, V., Collings, S., and Moss, B. (2005). Nitrate availability and hydrophyte species richness in shallow lakes. *Freshw. Biol.* 50 (6), 1049–1063. doi: 10.1111/j.1365-2427.2005.01375.x

Janssen, A. B. G., Hilt, S., Kosten, S., Klein, J. J.M., de Paerl, H. W., and Van de Waal, D. B. (2021). Shifting states, shifting services: linking regime shifts to changes in ecosystem services of shallow lakes. *Freshw. Biol.* 66 (1), 1–12. doi: 10.1111/fwb.13582

Kazanjian, G., Velthuis, M., Aben, R., Stephan, S., Peeters, E. T.H.M., and Frenken, T. (2018). Impacts of warming on top-down and bottom-up controls of periphyton production. *Sci. Rep.* 8, 9901. doi: 10.1038/s41598-018-26348-x

Kieft, B., Li, Z., Bryson, S., Hettich, R. L., Pan, C., Mayali, X., et al. (2021). Phytoplankton exudates and lysates support distinct microbial consortia with specialized metabolic and ecophysiological traits. *PNAS* 118 (41), e2101178118. doi: 10.1073/pnas.2101178118

Kim, Y., Lee, Y.-S., Wee, J., Hong, J., Lee, M., Kim, J. G., et al. (2021). Process-based modeling to assess the nutrient removal efficiency of two endangered hydrophytes:

linking nutrient-cycle with a multiple-quotas approach. *Sci. Total Environ.* 763, 144223. doi: 10.1016/j.scitotenv.2020.144223

Larras, F., Lambert, A.-S., Pesce, S., Rimet, F., Bouchez, A., Montuelle, B., et al. (2013). The effect of temperature and a herbicide mixture on freshwater periphytic algae. *Ecotoxicol. Environ. Saf.* 98, 162–170. doi: 10.1016/j.ecoenv.2013.09.007

Lips, S., Larras, F., and Schmitt-Jansen, M. (2022). Community metabolomics provides insights into mechanisms of pollution-induced community tolerance of periphyton. *Sci. Total Environ.* 824, 153777. doi: 10.1016/j.scitotenv.2022.153777

López Moreira, M., Gregorio, A., Toffolon, M., and Hölker, F. (2021). Hitting the sweet spot of complexity: reasons why the development of new custom-tailored models is still warranted and should be encouraged in aquatic sciences. *J. Limnol.* 80 (3), 2035. doi: 10.4081/jlimnol.2021.2035

Lorente, C., Causapé, J., Glud, R. N., Hancke, K., Merchán, D., Muñiz, S., et al. (2015). Impacts of agricultural irrigation on nearby freshwater ecosystems: the seasonal influence of triazine herbicides in benthic algal communities. *Sci. Total Environ.* 503-504, 151–158. doi: 10.1016/j.scitotenv.2014.06.108

Lüdecke, D. (2019). _esc: effect size computation for meta analysis (Version 0.5.1) (Digital software publication). doi: 10.5281/zenodo.1249218

Mahdy, A., Hilt, S., Filiz, N., Beklioğlu, M., Hejzlar, J., Özkundakci, D., et al. (2015). Effects of water temperature on summer periphyton biomass in shallow lakes: a pan-European mesocosm experiment. *Aquat. Sci.* 77, 499–510. doi: 10.1007/s00027-015-0394-7

Meerhoff, M., Audet, J., Davidson, T. A., De Meester, L., Hilt, S., Kosten, S., et al. (2022). Feedback between climate change and eutrophication: revisiting the allied attack concept and how to strike back. *Inland Waters* 12 (2), 187–204. doi: 10.1080/20442041.2022.2029317

Mooij, W. M., Janse, J. H., De Senerpont Domis, L. N., Hülsmann, S., and Ibelings, B. W. (2007). Predicting the effect of climate change on temperate shallow lakes with the ecosystem model PCLake. *Hydrobiologia* 584, 443–454. doi: 10.1007/s10750-007-0600-2

Nielsen, J. (2006). "Microbial process kinetics," in *Basic biotechnology*, 3. ed, vol. pp. Eds. C. Ratledge and B. Kristiansen (New York: Cambridge University Press), 155–180.

Olsen, S., Chan, F., Li, W., Zhao, S., Søndergaard, M., and Jeppesen, E. (2015). Strong impact of nitrogen loading on submerged macrophytes and algae: a long-term mesocosm experiment in a shallow Chinese lake. *Freshw. Biol.* 60 (8), 1525–1536. doi: 10.1111/fwb.12585

O'Reilly, C. M., Sharma, S., Gray, D. K., Hampton, S. E., Read, J. S., Rowley, R. J., et al. (2015). Rapid and highly variable warming of lake surface waters around the globe. *Geophys. Res. Lett.* 42 (24), 10773–10781. doi: 10.1002/2015GL066235

Pedrosa Gomes, M., and Juneau, P. (2017). Temperature and light modulation of herbicide toxicity on algal and cyanobacterial physiology. *Front. Environ. Sci.* 5 (50). doi: 10.3389/fenvs.2017.00050

Polst, B. H., Hilt, S., Stibor, H., Hölker, F., Allen, J., Vijayararaj, V., et al. (2022). Warming lowers critical thresholds for multiple stressor-induced shifts between aquatic primary producers. *Sci. Total Environ.* 838 (Part 4), 156511. doi: 10.1016/j.SCITOTENV.2022.156511

R Core Team (2020). "R: a language and environment for statistical computing," in *R foundation for statistical computing* (Vienna, Austria: R Foundation for Statistical Computing). Available at: https://www.R-project.org/.

Sayer, C. D., Davidson, T. A., and Jones, J. I. (2010). Seasonal dynamics of macrophytes and phytoplankton in shallow lakes: a eutrophication-driven pathway fromplants to plankton? *Freshw. Biol.* 55, 500–513. doi: 10.1111/j.1365-2427.2009.02365.X
Scheffer, M., Hosper, S. H.H., Meijer, M.-L. L., Moss, B., and Jeppesen, E. (1993). Alternative equilibria in shallow lakes. *Trends Ecol. Evol.* 8 (8), 275-279. doi: 10.1016/0169-5347(93)90254-M

Schmitt-Jansen, M., Bley, L.-M., Krumbiegel, M.-L., and Rotter, S. (2016). Induced community tolerance of periphyton towards combined salt and toxic stress. *Freshw. Biol.* 61 (12), 2152-2161. doi: 10.1111/fwb.12799

Sharma, A., Kumar, V., Shahzad, B., Tanveer, M., Sidhu, G. P.S., Handa, N., et al. (2019). Worldwide pesticide usage and its impacts on ecosystem. *SN Appl. Sci.* 1, 1446. doi: 10.1007/s42452-019-1485-1

Tasmin, R., Shimasaki, Y., Tsuyama, M., Qiu, X., Khalil, F., Okino, N., et al. (2013). Elevated water temperature reduces the acute toxicity of the widely used herbicide diuron to a green alga, pseudokirchneriella subcapitata. *Environ. Sci. Pollut. Res.* 21, 1064-1070. doi: 10.1007/s11356-013-1989-y

Tlili, A., Berard, A., Blanck, H., Bouchez, A., Cássio, F., Eriksson, K. M., et al. (2016). Pollution-induced community tolerance (PICT): towards an ecologically relevant risk assessment of chemicals in aquatic systems. *Freshw. Biol.* 61 (12), 2141-2151. doi: 10.1111/fwb.12558

Vasconcelos, F. R., Diehl, S., Rodriguez, P., Hedström, P., Karlsson, J., Byström, P., et al. (2016). Asymmetrical competition between aquatic primary producers in a warmer and browner world. *Ecology* 97 (10), 25802592. doi: 10.1002/ecy.1487

Vijayaraj, V., Laviale, M., Allen, J., Amoussou, N., Hilt, S., Hölker, F., et al. (2022). Multiple-stressor exposure of aquatic food webs: nitrate and warming modulate the effect of pesticides. *Water Res.* 216, 118325. doi: 10.1016/j.watres.2022.118325

Wendt-Rasch, L., Van Den Brink, P. J., Crum, S. J. H., and Woin, P. (2004). The effects of a pesticide mixture on aquatic ecosystems differing in trophic status: responses of the macrophyte myriophyllum spicatum and the periphytic algal community. *Ecotoxicol. Environ. Saf.* 57, 383-398. doi: 10.1016/j.ecoenv.2003.09.010

Wijewardene, L., Wu, N., Qu, Y., Guo, K., Messyasz, B., Lorenz, S., et al. (2021). Influences of pesticides, nutrients, and local environmental variables on phytoplankton communities in lentic small water bodies in a German lowland agricultural area. *Sci. Total Environ.* 780, 146481. doi: 10.1016/j.scitotenv.2021.146481

Winkler, K., Fuchs, R., Rounsevell, M., and Herold, M. (2021). Global land use changes are four times greater than previously estimated. *Nat. Commun.* 12, 2501. doi: 10.1038/s41467-021-22702-2

Xu, Z., Zhang, X., Xie, J., Yuan, G., Tang, X., Sun, X., et al. (2014). Total nitrogen concentrations in surface water of typical agro- and forest ecosystems in China 2004-2009. *PloS One* 9 (3), e92850. doi: 10.1371/journal.pone.0092850

Zhang, P., Grutters, B. M.C., van Leeuwen, C. H. A., Xu, J., Petruzzella, A., van den Berg, R. F., et al. (2019). Effects of rising temperature on the growth, stoichiometry, and palatability of aquatic plants. *Front. Plant Sci.* 9, 1947. doi: 10.3389/fpls.2018.01947

Zhang, P., Wang, T., Zhang, H., Wang, H., Hilt, S., Shi, P., et al. (2022). Heat waves rather than continuous warming exacerbate impacts of nutrient. *Environ. Int.* 168, 107478. doi: 10.1016/j.envint.2022.107478PD