Impact of Nutrients, Temperatures, and a Heat Wave on Zooplankton Community Structure: An Experimental Approach

Shallow lakes are globally the most numerous water bodies and are sensitive to external perturbations, including eutrophication and climate change, which threaten their functioning. Extreme events, such as heat waves (HWs), are expected to become more frequent with global warming. To elucidate the effects of nutrients, warming, and HWs on zooplankton community structure, we conducted an experiment in 24 flow-through mesocosms (1.9 m in diameter, 1.0 m deep) imitating shallow lakes. The mesocosms have two nutrient levels (high (HN) and low (LN)) crossed with three temperature scenarios based on the Intergovernmental Panel on Climate Change (IPCC) projections of likely warming scenarios (unheated, A2, and A2 + 50%). The mesocosms had been running continuously with these treatments for 11 years prior to the HW simulation, which consisted of an additional 5 ◦C increase in temperature applied from 1 July to 1 August 2014. The results showed that nutrient effects on the zooplankton community composition and abundance were greater than temperature effects for the period before, during, and after the HW. Before the HW, taxon richness was higher, and functional group diversity and evenness were lower in HN than in LN. We also found a lower biomass of large Cladocera and a lower zooplankton: phytoplankton ratio, indicating higher fish predation in HN than in LN. Concerning the temperature treatment, we found some indication of higher fish predation with warming in LN, but no clear effects in HN. There was a positive nutrient and warming interaction for the biomass of total zooplankton, large and small Copepoda, and the zooplankton: phytoplankton ratio during the HW, which was attributed to recorded HW-induced fish kill. The pattern after the HW largely followed the HW response. Our results suggest a strong nutrient effect on zooplankton, while the effect of temperature treatment and the 5 ◦C HW was comparatively modest, and the changes likely largely reflected changes in predation. Water 2020, 12, 3416; doi:10.3390/w12123416 www.mdpi.com/journal/water Water 2020, 12, 3416 2 of 19


Introduction
Lakes, the overwhelming majority of which are shallow, cover only 3% of the world's surface [1,2] but are biodiversity hotspots that provide significant ecosystem goods and services [3,4]. Today, shallow lakes are severely threatened by eutrophication, climate change, and their combined effects [5,6]. Global warming not only increases the average temperature but also enhances the frequency, intensity, and duration of extreme events such as heat waves [7][8][9][10]. Although no universal definition of "heat wave" (HW) exists [8], periods of abnormally hot weather affect the structure and function of freshwater systems. The frequency of HW is increasing [10], and HWs may affect both the physico-chemical conditions [11] and biota in lakes, as well as their interactions [12][13][14].
Zooplankton play an important role in the food web of shallow lakes, being placed at the intermediate trophic level linking fish and phytoplankton, and they contribute to maintaining clear water conditions and respond to eutrophication [15][16][17]. They are also sensitive to increasing temperatures, directly affecting their ontogeny and physiology and indirectly altering their sources of food [18][19][20]. Furthermore, several comparative studies have reported that increasing temperatures shape zooplankton communities through changes in fish communities, affecting the predation pressure on zooplankton (e.g., [21][22][23]). In heated mesocosms without fish, Strecker et al. [24] found a decrease in zooplankton biomass, especially that of large-sized cladocerans, and they suggested it to reflect a taxonomic shift in the composition of phytoplankton towards larger sizes, including toxic and less edible cyanobacteria, as seen in several studies [25]. A pan-European mesocosm experiment also showed a decrease in large-sized zooplankton along an increasing temperature gradient [26]. Temperature itself may also affect body size. The metabolic theory of ecology (MTE) suggests that body size decreases with increasing temperature [27], and declining body size is proposed to be a universal response to warming [28][29][30][31]. A study covering lakes at a latitude from 6 • to 74 • recorded a reduction in Cladocera body size with increasing temperatures [32], but in this and other studies [5,23] it was attributed to a parallel-occurring increase in fish predation.
Whether nutrient or temperature is of key importance in shaping the zooplankton community is debatable, but a recent cross-system analysis from mesocosms and lakes found an overall stronger effect of nutrients than of temperature on various biota in lakes and streams [33].
Several previous studies have revealed significant adverse effects of warming on freshwater ecosystems (e.g., [34][35][36][37]), but there are only a few studies with a focus on the impact of HW on lakes, particularly on the zooplankton. In this study, we explored the response of zooplankton community composition, biomass, diversity, and size structure to differences in nutrient concentration, temperature, and to a one-month artificial HW in mesocosms run at different temperatures according to the Intergovernmental Panel on Climate Change (IPCC) warming scenarios (A2 and A2 + 50%) crossed with low and high nutrient loading. The mesocosms used in our study are part of the longest running (since 2003) warming lake experimental facility in the world [38,39]. We hypothesized that: (1) eutrophication would have an overall stronger effect on the zooplankton community structure and functioning than relevant increases in temperature/HW due to global warming, (2) eutrophication would potentially increase the zooplankton biomass due to higher food availability unless planktivorous fish predation in the absence of piscivorous fish would prevent it. (3) Enhanced temperature and the HW would result in reduced zooplankton biomass, not least the large-bodied specimens, due to increased fish predation at the higher temperatures, unless the HW led to fish kill or fish spawning were affected negatively by higher temperature, as seen in some studies [34], and (4) the size diversity of the zooplankton community would be significantly reduced by an increase in nutrients and temperature/HW at both low and high nutrient concentrations due to selective predation on large-bodied specimens.

Experimental Setup
A mesocosm experimental facility in Lemming (Central Jutland, Denmark; 56 • 14 N, 9 • 31 E), which has been running since 2003, was used for our study. The setup contains 24 outdoor cylindrical stainless-steel tanks with a diameter of 1.9 m and a depth of 1.5 m in which the water is fully mixed by paddles. The sediment consists of a mixed 0.1 m layer of washed sand and a 0.1 m layer of nutrient-rich sediment (mainly coming from a nearby natural freshwater pond) that has been passed through a net (mesh size: 1-2 cm) to avoid uncontrolled vegetation and vertebrates. It is a flow-through system receiving tap water 6 times per day and has a water residence time of approximately 2.5 months.
The experimental design included a combination of three temperature scenarios and two different nutrient levels (four replicates): unheated (ambient, hereinafter A0) and the IPCC scenarios A2 (+3 • C) and A2 + 50% (+4.5 • C; hereinafter A2+) [40]. In the high-nutrient mesocosms (hereinafter HN), nutrients were added to half of the mesocosms on a weekly basis as Na 2 HPO 4 and Ca(NO 3 ) 2 solutions containing 2.7 mg phosphorus m −2 day −1 and 27.1 mg nitrogen m −2 day −1 . The low-nutrient mesocosms (hereinafter LN) received unenriched tap water (total phosphorus (TP) concentration of 0.003 mg L −1 and total nitrogen (TN) concentration of 2.6 mg L −1 ). With the two different loadings we aimed at creating shallow lakes in the two typical states-a clearwater, macrophyte-dominated state and a turbid phytoplankton or filamentous algae-dominated state [16]. The temperature scenarios are within the range of the IPCC scenarios [40].
At initiation of the experiment in 2003, plankton, benthic animals, and sediments (including resting stages of biota) were collected from five nutrient-poor to nutrient-rich lakes, mixed and added in equal amounts to all mesocosms. Macrophytes (mainly Elodea canadensis Michx and Potamogeton crispus L.) appeared naturally in all the LN, while the HN generally had no or only sparse vegetation. Male Gasterosteus aculeatus (three-spined sticklebacks) captured from natural populations were stocked in each mesocosm in near to natural densities according to nutrient concentrations: 1 individual in the LN and 12 in the HN. Later, female three-spined sticklebacks were added to the HN to allow breeding. Three-spined stickleback is a planktivorous fish that is often used in such kinds of experiments (e.g., [34,41] as it is robust, and individuals stay small. Due to the size of the mesocosms it is not possible to add also piscivorous fish, as the predator-prey interactions would not be realistic. When observed, dead fish were replaced to maintain a constant fish density. A more detailed description of the experimental setup can be found in Liboriussen et al. [38].

Heat Wave Treatment
The current experiment was undertaken from 1 June to end of August 2014. From 1 July to 1 August, a HW was created with an extra +5 • C increment of temperature in the heated A2 and A2+ scenarios (ca. +8 • C and +9.5 • C above ambient temperature, respectively), followed by a return to the normal heating. No heating was applied in the ambient temperature mesocosms. In the remaining text, the HW simulation experiment is referred to as "Before the HW", "During the HW", and "After the HW" (Figure 1).

Sampling and Laboratory Analysis
Temperature was measured at 30 min intervals and alkalinity was determined by end-point titration with 0.1 N HCl according to Gran [42] and measured once every 2 or 3 days. Sampling for determination of nutrients (TP and TN) and chlorophyll a (Chl-a) was done twice a week during the HW and otherwise once a week. Standard methods were used for chemical variables [43]. Chl-a was estimated using ethanol extraction [44].
Zooplankton samples were taken on 26 and 30 June representing the before period of the HW, 2, 7, 14, 21, and 28 July (the HW period) and 4 and 28 August (the after period of the HW). Sampling was conducted with a 1 m long tube sampler (covering the whole water column) at 6 randomly selected sites and pooled. Eight liters of the composite sample were filtered through a 50 µ m mesh and preserved in 4% Lugol's solution for further analysis. Identification of all zooplankton taxa was carried out to genus or taxon level, whenever possible. Copepods were separated into cyclopoids and calanoids (including adults and copepodites) and nauplii (without further determination). Counting and identification of the samples were done using a stereomicroscope (Leica M125) for cladocerans and copepods and an inverted microscope (Leica DMI 4000 B) for rotifers and nauplii. We measured the body size of about 25 individuals of each taxon, when possible, and calculated the body weight from length-weight allometric relationships [45][46][47][48]. For rotifers, we used the geometric formulae suggested by Ruttner-Kolisko [49]. We calculated biovolumes and converted them to dry weight according to Dumont et al. [45], Ruttner-Kolisko [49] and Malley et al. [50]. Phytoplankton samples were taken with the same frequency as the zooplankton samples, identification was undertaken to at least genus level, followed by counting and conversion to biomass as described in Filiz et al. [51]. Macrophyte abundance was quantified as percent volume of the water column inhabited by plants (plant volume inhabited (PVI, %)). Percentage cover and height of the submerged plants were assessed, allowing estimation of the proportion of the water column occupied by submerged plants (for further details see Davidson et al. [52]).

Stability Analyses
Following the procedure in Hillebrand et al. [53], we measured resistance in zooplankton total biomass as functional stability during the disturbance. Considering the fact that our HW lasted a month (extended pulse disturbance), we included data from all five samples taken during the HW for the resistance analysis. Resistance (as the differences between disturbed and undisturbed communities) was calculated by log response ratio (LRR). LRR was applied to zooplankton total biomass between undisturbed (just before the HW) and disturbed (during the HW) mesocosms within the treatments. Slope of regression over time was calculated in R with the 'lm' function in the "stats" package [54] to see the change in resistance during the time that HW was applied. One sample t-test was used to check if there was a significant difference between the benchmarks and the functional stability (resistance) with the 't.test' function in R [53,54]. In order to gain benchmark

Sampling and Laboratory Analysis
Temperature was measured at 30 min intervals and alkalinity was determined by end-point titration with 0.1 N HCl according to Gran [42] and measured once every 2 or 3 days. Sampling for determination of nutrients (TP and TN) and chlorophyll a (Chl-a) was done twice a week during the HW and otherwise once a week. Standard methods were used for chemical variables [43]. Chl-a was estimated using ethanol extraction [44].
Zooplankton samples were taken on 26 and 30 June representing the before period of the HW, 2, 7, 14, 21, and 28 July (the HW period) and 4 and 28 August (the after period of the HW). Sampling was conducted with a 1 m long tube sampler (covering the whole water column) at 6 randomly selected sites and pooled. Eight liters of the composite sample were filtered through a 50 µm mesh and preserved in 4% Lugol's solution for further analysis. Identification of all zooplankton taxa was carried out to genus or taxon level, whenever possible. Copepods were separated into cyclopoids and calanoids (including adults and copepodites) and nauplii (without further determination). Counting and identification of the samples were done using a stereomicroscope (Leica M125) for cladocerans and copepods and an inverted microscope (Leica DMI 4000 B) for rotifers and nauplii. We measured the body size of about 25 individuals of each taxon, when possible, and calculated the body weight from length-weight allometric relationships [45][46][47][48]. For rotifers, we used the geometric formulae suggested by Ruttner-Kolisko [49]. We calculated biovolumes and converted them to dry weight according to Dumont et al. [45], Ruttner-Kolisko [49] and Malley et al. [50]. Phytoplankton samples were taken with the same frequency as the zooplankton samples, identification was undertaken to at least genus level, followed by counting and conversion to biomass as described in Filiz et al. [51]. Macrophyte abundance was quantified as percent volume of the water column inhabited by plants (plant volume inhabited (PVI, %)). Percentage cover and height of the submerged plants were assessed, allowing estimation of the proportion of the water column occupied by submerged plants (for further details see Davidson et al. [52]).

Stability Analyses
Following the procedure in Hillebrand et al. [53], we measured resistance in zooplankton total biomass as functional stability during the disturbance. Considering the fact that our HW lasted a month (extended pulse disturbance), we included data from all five samples taken during the HW for the resistance analysis. Resistance (as the differences between disturbed and undisturbed communities) was calculated by log response ratio (LRR). LRR was applied to zooplankton total biomass between undisturbed (just before the HW) and disturbed (during the HW) mesocosms within the treatments. Slope of regression over time was calculated in R with the 'lm' function in the "stats" package [54] to see the change in resistance during the time that HW was applied. One sample t-test was used to check if there was a significant difference between the benchmarks and the functional stability (resistance) with the 't.test' function in R [53,54]. In order to gain benchmark values for resistance, instead of 0 (referring to no change between the mesocosms), we used data from just before the HW. Since the treatments (A2 and A2+) had already been held at higher than ambient temperatures for 11 years prior to the HW, the initial conditions between the treatments and their control were not the same. We did not calculate recovery or resilience since the treatments were resistant to the disturbance.
Linear mixed models (LMMs) were used to determine the treatment effects on zooplankton biomass, taxonomic composition, community parameters and size structure using the 'lmer' function in the "lme4" package in R [63]. Sampling date was used as a random effect, and the fixed effects were different temperatures (A2, A2+) and nutrient treatments (LN, HN). Due to the limited number of samples before and after the HW, we ran separate analyses for the three periods (before, during, and after the HW). After running the models, residuals were checked by Shapiro-Wilk's test (p > 0.05) for variance normality and homogeneity. If model residuals did not meet the normality assumption, box-cox transformation was applied to the dataset [64]. Interaction effects of different temperature and nutrient levels were selected according to the Akaike Information Criterion (AIC). If the models with interaction did not have lower AIC than those of the models without interaction, we did not take into consideration the interaction models in the analyses. To ascertain patterns between selected explanatory variables (temperature, PVI, Chl-a, alkalinity) and response variables (taxonomic grouping of zooplankton, taxon richness and evenness, and size diversity), co-inertia analysis (COI) was performed with the R package "ade4" [65] after testing for collinearity (8 being the accepted threshold, [66]) with the function 'vifstep' in the "usdm" package [67]. COI multivariate analysis can couple the two sets of data (here explanatory and response), following a Principle Component Analysis (PCA). Significance of the test acquired via free permutation test RV (number of random matching = 1000) [68,69].

Results
Daily average water temperatures for the three temperature scenarios from 1 June to 1 September are shown in Figure 1. Zooplankton sampling was undertaken from 26 June to 28 August. Before the HW in this time window, the average temperature was 17.3 ± 0.4 • C in A0, 19.8 ± 0.4 • C in A2, and 20.9 ± 0.4 • C in A2+. During the HW, the average temperature was 20.1 ± 2.2 • C in A0, 27.4 ± 2.2 • C in A2, and 28.7 ± 2.2 • C in A2+, and after the HW it was 18.2 ± 2.5 • C in A0, 21.9 ± 2.5 • C in A2, and 23.7 ± 2.5 • C in the A2+. During the zooplankton sampling period, average TP and TN concentrations were 0.404 ± 0.02 mg/L and 2.77 ± 0.13 mg/L, respectively, in HN, and 0.012 ± 0.001 mg/L and 0.284 ± 0.02 mg/L, respectively, in LN (see details in Table 1). Major changes were observed in the PVI in some of the treatments during the experimental period ( Figure 2). Linear regression revealed a significant decline with time in LN-A2 (p < 0.007, R 2 = 0.25) and HN-A0 (p < 0.02, R 2 = 0.19) and an increase in HN-A2 (p < 0.05, R 2 = 0.27), the latter due to an increase in filamentous algae (data not shown). During the HW 9-14 July, fish kills were observed in HN-A2 and HN-A2+ between 9-14 July, averaging 4.5 and 12 indiv.), respectively in the two treatments, while no fish kills were observed in HN-A0 and all the LN mesocosms.
Water 2020, 12, x FOR PEER REVIEW 6 of 20 Major changes were observed in the PVI in some of the treatments during the experimental period ( Figure 2). Linear regression revealed a significant decline with time in LN-A2 (p < 0.007, R 2 = 0.25) and HN-A0 (p < 0.02, R 2 = 0.19) and an increase in HN-A2 (p < 0.05, R 2 = 0.27), the latter due to an increase in filamentous algae (data not shown). During the HW 9-14 July, fish kills were observed in HN-A2 and HN-A2+ between 9-14 July, averaging 4.5 and 12 indiv.), respectively in the two treatments, while no fish kills were observed in HN-A0 and all the LN mesocosms.  Rotifera biomass and zooplankton taxon richness were higher in HN and the biomass of large Cladocera, and also the grazing potential of zooplankton, expressed as the zooplankton: phytoplankton (Zoop: Phyto) biomass ratio, was lower in HN than in LN ( Table 2, Figures 3 and 4). HN was dominated by small-bodied zooplankton taxa ( Table 2, Figures 3 and 4), and the phytoplankton community largely consisted of Microcystis sp. [51]. HN also had lower FGD and FGE (Table 2, Figure 3).
Concerning the temperature effects, numbers of small Cladocera and small Copepoda were

Community Structure before the Heat Wave
The initial zooplankton communities showed differences in structure between LN and HN, while there were only subtle differences among the temperature treatments within LN  Rotifera biomass and zooplankton taxon richness were higher in HN and the biomass of large Cladocera, and also the grazing potential of zooplankton, expressed as the zooplankton: phytoplankton (Zoop: Phyto) biomass ratio, was lower in HN than in LN ( Table 2, Figures 3 and 4). HN was dominated by small-bodied zooplankton taxa ( Table 2, Figures 3 and 4), and the phytoplankton community largely consisted of Microcystis sp. [51]. HN also had lower FGD and FGE (Table 2, Figure 3).
Concerning the temperature effects, numbers of small Cladocera and small Copepoda were significantly higher in LN-A2 and taxon evenness higher in LN-A2+ than in LN-A0 (Table 2, Figures 3  and 4). Temperature effects in HN were only seen for taxon richness in A2 (Table 2). Moreover, we found a significant negative interaction between nutrients and A2 warming for small Copepoda, and a significant positive interaction with A2+ warming for total zooplankton biomass ( Table 2).

Community Structure during and after the Heat Wave
During and after the HW, the nutrient effect largely followed the pattern observed before the HW. However, during the HW, we also found a significantly lower size diversity as well as taxon evenness and taxon diversity in HN than in LN, and a higher biomass of small cladocerans ( Table 2, Figures 3 and 4). Considering the temperature effects, there was a significantly higher Rotifera biomass in LN-A2+ than in LN-A0, accompanied by a significantly lower size diversity, Zoop: Phyto biomass ratio, taxon evenness and diversity as well as lower FGD and FGE. In HN-A2+, the biomasses of total zooplankton and small Copepoda were significantly higher than in HN-A0 (Table 2). In addition, there was a positive nutrient and warming interaction in HN-A2+ for the biomass of total zooplankton, large Copepoda, small Copepoda, and the Zoop: Phyto ratio during the HW ( Table 2).
The pattern for the nutrient effect after the HW largely followed the HW response, except that the nutrient effect on small cladocerans was no longer significant. Concerning temperature, the biomass of large Cladocera, size diversity, and taxon evenness were significantly lower in LN-A2+ than in LN-A0 (Table 2, Figures 3 and 4). Taxon richness was significantly lower in both HN-A2 and HN-A2+ than in HN-A0. Moreover, the biomasses of total zooplankton and large Copepoda and the Zoop: Phyto biomass ratio were significantly higher in HN-A2+ than in HN-A0 (Table 2, Figures 3 and 4). Interaction effects of nutrient and A2+ warming were significant for most variables, and there was a positive interaction effect for total zooplankton, large Cladocera, large Copepoda, and small Copepoda biomasses as well as for size diversity, the Zoop: Phyto biomass ratio, taxon evenness and taxon diversity. Taxon richness, however, was, negatively affected by the interaction between nutrient and A2+ warming ( Table 2).

Co-Inertia Analysis
The relationships between the selected explanatory variables (temperature, Chl-a, alkalinity, PVI) and the zooplankton community response variables (large and small Cladocera and Copepoda, Rotifera, taxon richness and evenness of zooplankton, and size diversity) were further evaluated in the COI analyses, which were run separately for HN and LN. Dissolved oxygen was excluded in this analysis as it generally was high (super-saturation except for a short term period during a cooling event during the middle of the heat wave, see [70]). Significant relationships (p < 0.001) were found, with RV-coefficients of 0.40 and with the variance explained by the first two axes of 94.0% and 5.1% for LN; for HN, RV-coefficients were much lower, 0.14, as was the variance explained by the first two axes, namely 67.4% and 26.2% ( Figure 5). In Figure 5, panels (c) and (f) show the concordance between the explanatory and response variables and the distributions before, during, and after the HW are displayed by color gradients. The relative location of samples within the two datasets (explanatory and response) are presented by arrows in Figure 5, panels (c) and (f). Even though the grouping of the treatments that was exposed to HW (A2 and A2+ during the HW) was not clearly distinguishable from the others for either LN or HN, these treatments were grouped with A2 and A2+ from after the HW in LN (Figure 5c).

Stability Analyses
The stability analysis revealed no differences compared to the benchmark during the HW in LN-A2 (thus showing resistance to the HW disturbance), whereas LN-A2+ deviated, showed less resistance, as from the second week of the HW and significantly so in the last week of sampling (Figure 6a,b). HN-A2 similarly deviated from the benchmark after the third week, with a significant change occurring in the fourth and fifth weeks of the HW, thus became less resistant. HN-A2+ showed no significant change from the benchmark but a similar trend as in HN-A2 (Figure 6c,d), thus remaining resistant. The slope of change tended to be positive for all treatments during the HW, but only significantly so for HN-A2 (Figure 6e). (c,f) represent the concordance between panels a and b (color codes: cyan gradient-before period (A0: light cyan, A2: cyan, A2+: dark cyan), yellow gradient-during period (A0: light yellow, A2: yellow, A2+: dark yellow), pink gradient-after period (A0: light pink, A2: pink, A2+: dark pink)). Arrows in (f,c) display the relative position of samples within the two datasets (explanatory and response).
In the LN ordination plots, PVI was an important explanatory variable, not least of all for large and small Copepoda, which were strongly associated with PVI. Large Cladocera and taxon richness, small Cladocera and size diversity were also associated with PVI and/or low Chl-a and with temperature. Taxon evenness was negatively associated with temperature, while Rotifera were positively associated with temperature and Chl-a (Figure 5a,b). The pattern in HN followed that in LN, although the relationship was weaker.

Stability Analyses
The stability analysis revealed no differences compared to the benchmark during the HW in LN-A2 (thus showing resistance to the HW disturbance), whereas LN-A2+ deviated, showed less resistance, as from the second week of the HW and significantly so in the last week of sampling (Figure 6a,b). HN-A2 similarly deviated from the benchmark after the third week, with a significant change occurring in the fourth and fifth weeks of the HW, thus became less resistant. HN-A2+ showed no significant change from the benchmark but a similar trend as in HN-A2 (Figure 6c,d), thus remaining resistant. The slope of change tended to be positive for all treatments during the HW, but only significantly so for HN-A2 (Figure 6e).

Discussion
Our results clearly showed that the effects of eutrophication on the zooplankton community composition and abundance were much larger than those observed for the temperature treatments in all three periods (before, during, and after the HW), supporting our first hypothesis. This concurs with a recent cross-system analysis from mesocosms and lakes, finding an overall stronger effects of nutrients than of temperature on various biota [33]. However, our results may also reflect the fact that the variation in nutrients far exceeded the variation in temperature, even though the temperature

Discussion
Our results clearly showed that the effects of eutrophication on the zooplankton community composition and abundance were much larger than those observed for the temperature treatments in all three periods (before, during, and after the HW), supporting our first hypothesis. This concurs with a recent cross-system analysis from mesocosms and lakes, finding an overall stronger effects of nutrients than of temperature on various biota [33]. However, our results may also reflect the fact that the variation in nutrients far exceeded the variation in temperature, even though the temperature increase was relatively large. The COI analysis, moreover, revealed that PVI was a key explanatory variable in LN, in particular for large and small Copepoda, which were strongly associated with PVI (Figures 2 and 5). Large Cladocera, taxon richness, small Cladocera, and size diversity were also associated with PVI and/or low Chl-a. It is well established that macrophytes act both as a habitat and a refuge against fish predation, and the biofilm on the plant surface provides food for some cladocerans and copepods (e.g., [71,72]). Although the effect of PVI in HN largely followed the pattern in LN, it was weaker, which is to be expected as the PVI was overall much lower in HN, except for HN-A0 ( Figure 2) and in eutrophic lakes in general [16].
We found no major difference in total zooplankton biomass among the treatments, it being higher only in the HN mesocosms heated to the A2+ warming scenario during and after the HW. We would have expected (Hypothesis 2) an increase in zooplankton biomass from LN to HN conditions, as typically seen in shallow Danish lakes [73]. In HN, fish have been allowed to breed freely since 2006, which likely has created a high, though variable, predation pressure on the zooplankton [20] and a shift to Rotifera dominance as observed in eutrophic shallow Danish lakes [73]. High fish predation is known to alter the zooplankton community in favor of rotifers [73,74] and this may explain the overall lower biomass and the lower Zoop: Phyto biomass ratio in HN compared with LN [35]. Our mesocosms were too small to have also piscivorous fish included, which may have contributed to the modest effect of nutrients on the zooplankton biomass (Hypothesis 2). However, the deteriorated food quality caused by cyanobacteria dominance in these mesocosms [51] may also have been a contributory factor [75][76][77].
Before the HW, taxon diversity, evenness, and size diversity were not affected by nutrients. Taxon richness, however, was higher, and functional group diversity and evenness were lower in HN. This pattern was also found both during and after the HW, but in addition taxon diversity, evenness, and size diversity were lower in HN than in LN, emphasizing the effect of nutrients on community structure and function especially exposed to the HW (Hypotheses 1, 3 and 4).
Prior to the HW, we found only few significant temperature treatment effects. Small Cladocera and small Copepoda were higher in LN-A2 and large Cladocera lower in LN-A2+ compared with the controls (Table 2), suggesting a higher predation pressure with warming reflects size-selective predation (Hypothesis 3). Such a temperature effect on predation is observed both in cross-system mesocosm studies and in field data [23,26]. We also found a significantly lower taxon richness in HN-A2 and an increase in taxon evenness in LN-A2+ compared with their respective controls, but otherwise there was no observable warming effect.
During the HW, major temperature treatment effects were found in LN-A2+ (higher rotifer biomass; lower taxon diversity, taxon evenness, size diversity, FGD, FGE, and Zoop: Phyto ratio), and these changes persisted for some of the variables after the HW ( Table 2). The higher contribution of rotifers and higher Zoop: Phyto ratio suggesting increasing predation [19] likely reflect the higher temperature, as we did not find any significant reduction in PVI in LN-A2+ (Figure 2), which otherwise could have led to higher predation risk [71].
Furthermore, HN-A2+ showed a different pattern during and after the HW from the before period, having higher total zooplankton biomass (during and after the HW) and higher biomass of large cladocerans and Zoop: Phyto ratio than HN-A0. Contrary to the findings in LN, the changes in HN-A2+ indicate less fish predation. Accordingly, we observed fish kill in HN-A2 and HN-A2+, 9-14 days into the HW, but being much larger in HN-A2+ (Hypothesis 3). The remaining fish population may also have suffered warming effects as Moran et al. [78] experimentally demonstrated a negative effect of warming on three-spined stickleback survival and reproduction, with a 60% reduction of stickleback biomass at a temperature increase of 4 • C above ambient summer temperatures. Moreover, in their experiment the combined effect of heating and enhanced nutrient loading led to a further stickleback population decline, attributed to indirect effects of warming such as low oxygen concentrations. Unfortunately, we have no fish data from this study (apart from data on fish kill), but an earlier study by Šorf et al. [79] in the same mesocosms revealed that the lowest adult stickleback catch per unit effort (CPUE) and the highest recruitment occurred in the unheated HN mesocosms, suggesting that stickleback recruitment declined at the higher temperatures. Concurring with these results, Mehlis and Bakker [80] found that while less sperm fertilized more eggs at higher temperature (25 • C), the overall percentage of motile sperm significantly decreased. A negative effect on fish due to HW might, therefore, partly explain why the zooplankton biomass increased and crustaceans became more abundant in HN-A2+ during and after the HW (Hypothesis 3), while no such effect is to be expected in LN due presence of only one male stickleback.
The overall modest effect of temperature (apart from the highest temperature, A2+) concurs with results from a series of mesocosm experiments conducted in Liverpool, UK (e.g., [34,81]) as well as an earlier study of zooplankton in our mesocosms [20]. From the Liverpool experiment, McKee et al. [81] concluded that the effect of increased temperature caused by global warming is not expected to significantly alter ecosystem functional structures and processes during the next century. However, microevolution may be a contributing factor, considering that the zooplankton communities of our mesocosms have been facing warming for 11 years prior to this experiment. This is supported by a study in our mesocosm systems by Van Doorslaer et al. [82], which showed rapid (within a year) microevolution of S. vetulus. They found that clones exposed for 1 year to the A2+ treatment had about a five times higher survival rate at 26 • C than that of clones exposed to the ambient regime. Moreover, age at first reproduction and offspring numbers were affected, and the authors concluded that the S. vetulus populations would probably not be able to maintain themselves under global warming. Microevolution may thus help buffer changes in community structure under global warming. Moreover, although some zooplankton taxa are considered to be sensitive to temperature changes [83], cosmopolitan taxa such as C. sphaericus, A. rectangula, B. calcyflorus, and K. quadrata dominated the zooplankton community in the high-nutrient mesocosms [84][85][86][87][88], which may explain the resistance to increased temperature. Additionally, C. sphaericus was the most abundant Cladocera taxon in the HN-A2+ mesocosms and has an advantage by having a short egg development time, allowing fast development of populations under favorable environmental conditions [89]. Moreover, several previous experimental studies have shown tolerance of selected cladocerans to temperatures up to 30 • C [90,91]. An investigation by Donze [92] demonstrated an effect of heat on the survival of cladocerans at 30 • C with 0% mortality, while at 34 • C mortality was around 75%. The temperature variation in our mesocosms did not exceed these temperatures during the experimental period (Figure 1), and cladocerans might, therefore, have shown resistance to the heating scenarios (A2 and A2+ according to IPCC), though heat-shock proteins in zooplankton may somehow prevent cellular damage and constitute a molecular defense mechanism against all kinds of stressors [93,94].
The resistance analysis and the change in the resistance during HW revealed that all the treatments tended to have over-performance (i.e., higher biomass than in the benchmark for a given treatment) by the end of the HW disturbance (Figure 6a,d), which may be explained by the increased resource availability [51]. However, our results mostly showed insignificant changes in the different treatments, indicating overall high resistance to the HW disturbance. In LN-A2+ and HN-A2, however, the resistance in the beginning of the HW did not last until the end of the HW period.
We acknowledge that our study has some limitations regarding elucidation of the HW effect, as the HW was run without true controls to account for seasonal effects. Differences observed among the three periods (before, during, and after the HW) should therefore be interpreted with caution. However, the rather modest effect of the HW, despite a 5 • C increase in temperature resulting in, respectively, an 8 and 9.5 • C increase above the temperature in the ambient mesocosm, suggests that our conclusion of high nutrient but lower temperature and HW effects is valid. However, fish kill in HN-A2 and not least in HN-A2+ enhanced the effect of the HW. Modest effects of the HW were also found for phytoplankton in the same experiment [51], whereas temperature changes had a substantial effect on the biomass of bacterioplankton, heterotrophic flagellates, ciliates [95], and macrophytes [13], oxygen metabolism (gross ecosystem production and ecosystem respiration), and bacterioplankton production [96], as well as on greenhouse gases dynamics [70]. These results suggest that microbial communities and ecosystem processes are more sensitive to short-term HWs than the larger-bodied and more slow-growing organisms, such as zooplankton, which are likely to exhibit a delayed or weak response, unless a major fish kill occurs.

Conclusions
Our results indicate that the nutrients' effect on zooplankton community composition was much higher than that of temperature for three periods before, during, and after the heat wave (HW). They further showed that the HW had only modest effects on the zooplankton community. Our results indicate that a rise in the number of HWs and global temperatures may potentially weaken the capacity of maintaining or establishing a clear water state by shifting the zooplankton community towards dominance of small-sized taxa with a lower grazing potential, at least in systems that lack piscivorous fish, provided that the zooplanktivorous fish present are not subjected to fish kill at the higher temperature of HWs as seen in the HN-high temperature mesocosms in our experiment. Reduction of the external nutrient loading can, therefore, somehow compensate for the effect of warming due to climate change and help prevent fish kills during heat waves.