Phytoplankton Community Response to Nutrients, Temperatures, and a Heat Wave in Shallow Lakes: An Experimental Approach

Phytoplankton usually responds directly and fast to environmental fluctuations, making them useful indicators of lake ecosystem changes caused by various stressors. Here, we examined the phytoplankton community composition before, during, and after a simulated 1-month heat wave in a mesocosm facility in Silkeborg, Denmark. The experiment was conducted over three contrasting temperature scenarios (ambient (A0), Intergovernmental Panel on Climate Change A2 scenario (circa +3 ◦C, A2) and A2+ %50 (circa +4.5 ◦C, A2+)) crossed with two nutrient levels (low (LN) and high (HN)) with four replicates. The facility includes 24 mesocosms mimicking shallow lakes, which at the time of our experiment had run without interruption for 11 years. The 1-month heat wave effect was simulated by increasing the temperature by 5 ◦C (1 July to 1 August) in A2 and A2+, while A0 was not additionally heated. Throughout the study, HN treatments were mostly dominated by Cyanobacteria, whereas LN treatments were richer in genera and mostly dominated by Chlorophyta. Linear mixed model analyses revealed that high nutrient conditions were the most important structuring factor, which, regardless of temperature treatments and heat waves, increased total phytoplankton, Chlorophyta, Bacillariophyta, and Cyanobacteria biomasses and decreased genus richness and the grazing pressure of zooplankton. The effect of temperature was, however, modest. The effect of warming on the phytoplankton community was not significant before the heat wave, yet during the heat wave it became significant, especially in LN-A2+, and negative interaction effects between nutrient and A2+ warming were recorded. These warming effects continued after the heat wave, as also evidenced by Co-inertia analyses. In contrast to the prevailing theory stating that more diverse ecosystems would be more stable, HN were less affected by the heat wave disturbance, most likely because the dominant phytoplankton group cyanobacteria is adapted to high nutrient conditions and also benefits from increased temperature. We did not find any significant change in phytoplankton size diversity, but size evenness decreased in HN as a result of an increase in the smallest and largest size classes simultaneously. We conclude that the phytoplankton community was most strongly affected by the nutrient level, but less sensitive to changes in both temperature treatments and the heat wave simulation in these systems, which have been adapted for a long time to Water 2020, 12, 3394; doi:10.3390/w12123394 www.mdpi.com/journal/water Water 2020, 12, 3394 2 of 21 different temperatures. Moreover, the temperature and heat wave effects were observed mostly in LN systems, indicating that the sensitivity of phytoplankton community structure to high temperatures is dependent on nutrient availability.

better conservation and restoration strategies [54]. As Tilman [55] stated, more diverse systems tend to have greater ecosystem stability and greater resistance to a disturbance, reflecting interspecific complementarity, higher resource use efficiency, fewer diseases, and increased nutrient stores and supply rates in the long term. Having simpler food webs and lower biodiversity, eutrophic systems can be expected to show a lower level of ecosystem stability in response to warming than non-eutrophic systems. Limited studies of warming effects on oligotrophic lakes can show opposite results. For instance, while Yvon-Durocher et al. [56] presented higher phytoplankton richness with warming, Verbeek et al. [57] found lower richness in experimental enclosures. The multifaceted concept of ecosystem stability includes several features (e.g., resistance, resilience, persistence, robustness) to better characterize stability [58,59]. Resistance, the amplitude of change in the biomass of the ecosystems in response to a disturbance [60], recovery, the ability to return to their pre-disturbance state, and resilience, the time needed to return to their pre-disturbance state [61], are important aspects of the ecosystem stability to understand the system response.
In this study, we investigated the response of phytoplankton biomass, community structure and size structure to different temperatures, before, during, and after a 1-month experimentally induced heat wave at two contrasting nutrient loadings. The mesocosms used are established at the longest running warming experimental facility in the world (running since 2003) [62,63]. We hypothesized that (1) nutrients play a key role in determining phytoplankton community structure and biomass (2) heat waves would change the community structure, phytoplankton richness, and evenness, and likely lead to an increase in cyanobacteria and total phytoplankton biomass, with a stronger effect at high rather than low temperatures and high nutrient concentrations, (3) the mesocosms with low nutrient concentrations would show higher stability to heat waves compared with the nutrient enriched mesocosms due to their higher phytoplankton diversity, and (4) cell size distribution of phytoplankton community would be significantly reduced in the long term heated treatments and also during the heat wave, regardless of the nutrient scenarios.

Experimental Setup
The study was conducted in a mesocosm experimental facility (Lemming, Central Jutland, Denmark; 56 • 14 N, 9 • 31 E), which has been running since 2003. It is the longest running heated freshwater mesocosm setup in the world and encompasses of 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 mesocosms are continuously supplied with tap water and have a water residence time of approximately 2.5 months [62]. The sediment of enclosures, which were added in 2003, consists of a mixed 0.1 m layer of washed sand and a 0.1 m layer of lake sediment. Lake sediment collected from both nutrient-rich, phytoplankton-dominated and nutrient-poor, macrophyte-dominated Danish lakes and sieved through a net (mesh size: 1-2 cm) to remove vegetation and macroinvertebrates before it was added to all mesocosms equally ( Figure 1).
The experiment design combined three temperature treatments with two nutrient levels (high nutrient, HN; and low nutrient, LN) in a factorial design (four replicates of each). Nutrients were added weekly to the HN as Na 2 HPO 4 and Ca(NO 3 ) 2 solutions (2.7 mg/m 2 day −1 phosphorus and 27.1 mg/m 2 day −1 nitrogen), as in the years prior to our HW experiment. The temperature treatments consisted of unheated ambient mesocosms (where water temperatures fluctuated with the air temperature, hereinafter referred to A0) and mesocosms heated according to IPCC's A2 and A2+50% (hereinafter A2+) scenarios [64]. A2 was 2-4 • C warmer and A2+ 4-6 • C warmer than the ambient tanks ( Figure 1). The seasonal variations between the temperature treatments are based on future climate predictions for Denmark reduced to monthly resolution using the reference period 1961-1990 [65]. near to natural densities relative to nutrient concentrations: 1 individual in LN and 12 in HN. Later, female three-spine stickleback were introduced to HN to allow breeding. Since planktivorous fish three-spined stickleback is robust and stay small, it is often used for the mesocosm experiments [66]. A more detailed description of the experiment setup can be found in Liboriussen et al. [62]. We applied a 1-month (1 July-1 August) heat wave (hereinafter HW) in this experiment to the long-term heated mesocosms, run at two contrasting nutrient levels. The already-heated mesocosms (A2 and A2+) were exposed to an additional increase of 5 °C above the ambient temperature (in A0), leading to a circa +8 °C (A2) and +9.5 °C (A2+), respectively, higher temperature than in the A0 temperatures during the HW (A0s were not exposed to HW at any time since the mesocosm establishment in 2003). When the HW experiment was stopped, the temperature differences in the heated mesocosms was lowered to pre-HW values. The current experiment was undertaken from 26 June to 28 August 2014, which included three periods: before, during, and after the HW (Figure 1).

Sampling and Analyses of Physical and Chemical Variables
For phytoplankton analyses, 10 L water samples were taken with a tube sampler (V = 3.6 L) from different randomly chosen parts of each mesocosms and pooled. A total of 50 ml was taken from the sampled water, fixed with 1 ml Lugol's solution, and stored. Samples taken on 26 and 30 June represented the before period of the HW, samples from 2, 7, 14, 21, and 28 July the HW period (during), and samples from 4 and 28 August the after HW period. The same person performed the counting and identification. Phytoplankton counts followed the protocol developed in the WISER project (EU FP7, No.226273). For each phytoplankton taxon, at least 10 individuals were measured, except some rare species, using the Leica image analysis program to calculate phytoplankton size and biovolume. Zooplankton samples were taken with the same frequency as the phytoplankton samples. The zooplankton data used are from Işkın et al. [67].
Sampling for chlorophyll a (Chl-a) and nutrient concentrations were conducted twice a week during the 1-month HW period and once a week during the remaining experimental periods. Macrophyte abundance was measured as plant volume inhabited (PVI, %), as described in Davidson et al. [68]. Chl-a was determined by filtering water on Whatman GF⁄C filters and undertaking Macrophytes (mainly Elodea canadensis Michx and Potamogeton crispus L.) appeared naturally in all LN mesocosms, while in the HN mesocosms, submerged plants were generally sparse or absent. The distribution and coverage of macrophytes were manipulated before initiating the experiment by removing floating-leaved macrophytes and duckweed to prevent their dominance. Male Gasterosteus aculeatus (three-spine sticklebacks) captured from natural populations were stocked in each tank in near to natural densities relative to nutrient concentrations: 1 individual in LN and 12 in HN. Later, female three-spine stickleback were introduced to HN to allow breeding. Since planktivorous fish three-spined stickleback is robust and stay small, it is often used for the mesocosm experiments [66]. A more detailed description of the experiment setup can be found in Liboriussen et al. [62].
We applied a 1-month (1 July-1 August) heat wave (hereinafter HW) in this experiment to the long-term heated mesocosms, run at two contrasting nutrient levels. The already-heated mesocosms (A2 and A2+) were exposed to an additional increase of 5 • C above the ambient temperature (in A0), leading to a circa +8 • C (A2) and +9.5 • C (A2+), respectively, higher temperature than in the A0 temperatures during the HW (A0s were not exposed to HW at any time since the mesocosm establishment in 2003). When the HW experiment was stopped, the temperature differences in the heated mesocosms was lowered to pre-HW values. The current experiment was undertaken from 26 June to 28 August 2014, which included three periods: before, during, and after the HW (Figure 1).

Sampling and Analyses of Physical and Chemical Variables
For phytoplankton analyses, 10 L water samples were taken with a tube sampler (V = 3.6 L) from different randomly chosen parts of each mesocosms and pooled. A total of 50 ml was taken from the sampled water, fixed with 1 ml Lugol's solution, and stored. Samples taken on 26 and 30 June represented the before period of the HW, samples from 2, 7, 14, 21, and 28 July the HW period (during), and samples from 4 and 28 August the after HW period. The same person performed the counting and identification. Phytoplankton counts followed the protocol developed in the WISER project (EU FP7, No.226273). For each phytoplankton taxon, at least 10 individuals were measured, except some rare species, using the Leica image analysis program to calculate phytoplankton size and biovolume. Zooplankton samples were taken with the same frequency as the phytoplankton samples. The zooplankton data used are from Işkın et al. [67].
Sampling for chlorophyll a (Chl-a) and nutrient concentrations were conducted twice a week during the 1-month HW period and once a week during the remaining experimental periods. Macrophyte abundance was measured as plant volume inhabited (PVI, %), as described in Davidson et al. [68]. Chl-a was determined by filtering water on Whatman GF⁄C filters and undertaking spectrophotometric analysis after ethanol extraction [69]. Standard methods were used for chemical variables [70]. Measurements of temperatures were conducted at 30 min intervals using probes and averaged to 24 h means (see [62]).

Statistical Analyses
Genus evenness was calculated according to the Pielou's evenness index (J) from the abundance data using the "vegan" package [71]. Richness was an expression of the total number of all identified genera. Resource use efficiency (RUE) was calculated as ln[Chl-a:TP] (units were taken as µg Chl-a:mgTP) as an indicator of changes in phytoplankton biomass in response to the HW at low and high nutrient concentrations. Although both Chl-a and phytoplankton biomass can be used to calculate RUE, we choose to use Chl-a as it is an estimate independent of the microscopic analysis [72]. The difference between initial nutrient conditions of low and high nutrient mesocosms (for both TP and total nitrogen (TN)) were tested by analyses of variances (ANOVA) in R with the "aov" function in "stats" package [73].

Effects of Nutrients, Warming and Heat Wave
Linear mixed models (LMMs) were used to determine treatment effects on the biomasses of total phytoplankton and the phytoplankton groups, genus richness and evenness, RUE, grazing pressure, and size structure using the "lmer" function in the "lme4" package [74] in R. Models were run separately for each period (before, during, and after the HW), as we had no true control treatment. Mesocosm replicates and time were used as random effects in all the models, while the fixed effects were temperatures (A0, A2 and A2+) and nutrient concentrations (LN, HN). Models were run to test nutrient effects (LN-A0 vs HN-A0), temperature effects (A0 vs A2 and A2+ within LN and HN), and interaction effects (nutrient vs. A2 and A2+). Variance normality and homogeneity of model residuals were checked applying Shapiro-Wilk's test (p > 0.05). If model residuals did not meet the assumption of normality, datasets were scaled [75] following a box-cox transformation by using "car" package in R (family = bcPower, [76]). Models with an interaction effect between temperature and nutrient levels were selected according to the Akaike Information Criterion (AIC). If the AICs of models with interaction were higher than those of the models without interaction, they were not taken into consideration. Since Euglenophyta and Xanthophyta groups had many zero observations, we did not include them in the tests, to avoid potential violation of the analysis.
Co-inertia Analysis (COI) was used to reveal the patterns between explanatory (TP, TN, temperature, biomasses of Cladocera, Copepoda and Rotifera, and PVI) and phytoplankton response variables (biomasses of total phytoplankton, Chlorophyta, Cyanobacteria, and Bacillariophyta, RUE, genus richness, and evenness), separately for LN and HN. After testing for collinearity (accepted threshold was 8, [77]) with the function "vifstep" in the "usdm" package [78], we performed co-inertia (COI) analyses using the explanatory and response dataset simultaneously with the function "dudi.pca" in R "ade4" package [79]. The optimizing criterion in the COI analysis is the resulting sample scores (explanatory scores and response scores), which are the most covariant [80]. Significance of the COI analyses was tested with a free permutation test RV (number of random matching = 1000, significance threshold p = 0.05) [81]. Biomass of Rotifera, Cladocera and Copepoda were used individually instead of total zooplankton biomass since the feeding efficiency; thus, the grazing pressure differ for each group [82].

Stability Analyses
Resistance, resilience, and recovery were calculated for phytoplankton total biomass following Hillebrand et al. [59], who defined the resistance as "the ability to withstand the perturbation". The difference between perturbed and the control community was calculated using the log response ratio (LRR, resistance = ln(disturbed treatment/control treatment)) of the phytoplankton total biomass. Since our disturbance lasted a month (extended pulse disturbance), we used not only data from the first sampling event but from all five during the HW to calculate resistance. Additionally, we calculated the change in resistance during the disturbance over time as regression slope (resistance = intercept + slope × time) of phytoplankton biomass LRR in R programming with the "lm" function in "stats" package [73], which was also used for resilience calculations after HW stopped. Recovery was calculated, applying the same calculations of resistance. Resilience and recovery were calculated only for nonresistance treatments over time (LN-A2 and LN-A2+) for the last sampling date (28.08) to assess if there was a recovery [59].
In our study, the benchmark for resistance and recovery was not equal to 0 (implying no change between the mesocosms) because our experiment had been running for 11 years prior to the HW experiment, and treatments were already formed by adding nutrients and increasing temperatures for 11 years. Therefore, the benchmarks were calculated from the sampling date just before the HW (30 June, by LRR of disturbed (A2 and A2+) to control community (A0)). One sample t-test was used to check [59] if there was a significant difference between the benchmarks and the stability indicators with the "t.test" function in the "stats" package [73].

Phytoplankton Size Structure
To reveal the changes in phytoplankton size structure, two parameters of size distribution (size diversity and size evenness) were computed (according to [83][84][85]). The size diversity index was calculated based on the Shannon-Wiener diversity index [86] adapted to a continuous variable using individual size measurements (as length). Since the method uses a continuous probability density function for the probability estimation, negative values may occur, and they indicate extremely low size diversity. Since size diversity expresses both variability and regularity of the size distribution, size evenness (Je) was also computed to distinguish the regularity of the size distribution. It was calculated by dividing the size diversity exponential by its possible maximum for a given size range and ranged between 0 and 1. Both size metrics were based on individual abundance (not on biomass), which give a unique value per size distribution among treatments. For further investigation, which was performed only in HN since size evenness was significantly changed only in these treatments, phytoplankton abundances were grouped into four size classes (hereinafter SC) according to individual size measurements (SC1: <5 µm, SC2: 5-20 µm, SC3: 20-100 µm, and SC4: >100 µm) and their changes were tested by LMMs to determine treatment effects in different periods as described above.

Results
The temperatures of three periods and the average TP and TN concentrations for all treatments during the experimental period (26 June to 28 August) are given in Table 1.

Initial Conditions
Since the mesocosm facility has been running for 11 years prior to the HW experiment, there could have been differences between the nutrient conditions among the temperature treatments within LN and HN (Table 1). ANOVA analyses, however, showed no significant differences in TN and TP among the three temperature treatments neither in LN nor in HN (p > 0.05). Desired temperature differences were ensured throughout the study regardless of the nutrient levels. A2 treatments were 2-4 • C and A2+ treatments were 4-6 • C higher than A0s and during the HW A2 and A2+, treatments were heated by an additional 5 • C (Table 1, Figure 1), while no heating was done in A0. Due to a natural HW in July, the unheated treatment (A0) was, however, higher during the period of the HW than in the periods before and after the HW (Table 1, Figure 1).

Phytoplankton Community Composition Throughout the Study
In LN, Chlorophyta was the dominant phytoplankton group in terms of biomass throughout the study period for all temperature treatments. Their contribution ranged from 35% to 91% except for LN-A2 after the HW when Cryptophyta dominated (54% contribution). Moreover, in LN-A2 during the HW, Chlorophyta and Cryptophyta co-dominated (36% and 33%, respectively) the phytoplankton community. A more diverse community was found in LN (all treatments) where Chlorophyta, Cyanobacteria, Bacillariophyta and Cryptophyta contributed importantly to biomass in all temperature treatments during all three periods, except LN-A2+ after the HW, which was dominated by Chlorophyta ( Figure 2). In LN-A0 Pediastrum boryanum and Chlorella sp., in LN-A2 Cryptomonas sp. and Plagioselmis lacustris, and in LN-A2+, the green algae species Botryococcus braunii were the dominant species. In HN, Cyanobacteria were the dominant group throughout the study in all temperature treatments, with biomass contributions ranging from 35% to 75%, except for HN-A0, where Chlorophyta was the dominant group and constituted 48% and 49% of the total biomass before and during the HW, respectively ( Figure 2). Besides Cyanobacteria, there were large contributions of Chlorophyta and Bacillariophyta in all the HN treatments. Micractinium sp. and Microcystis sp. were the two dominant taxa in A0, while Microcystis sp. dominated A2 and A2+ in HN.
Water 2020, 12, x FOR PEER REVIEW 7 of 21 A2+, treatments were heated by an additional 5 °C (Table 1, Figure 1), while no heating was done in A0. Due to a natural HW in July, the unheated treatment (A0) was, however, higher during the period of the HW than in the periods before and after the HW (Table 1, Figure 1).

Phytoplankton Community Composition Throughout the Study
In LN, Chlorophyta was the dominant phytoplankton group in terms of biomass throughout the study period for all temperature treatments. Their contribution ranged from 35% to 91% except for LN-A2 after the HW when Cryptophyta dominated (54% contribution). Moreover, in LN-A2 during the HW, Chlorophyta and Cryptophyta co-dominated (36% and 33%, respectively) the phytoplankton community. A more diverse community was found in LN (all treatments) where Chlorophyta, Cyanobacteria, Bacillariophyta and Cryptophyta contributed importantly to biomass in all temperature treatments during all three periods, except LN-A2+ after the HW, which was dominated by Chlorophyta ( Figure 2). In LN-A0 Pediastrum boryanum and Chlorella sp., in LN-A2 Cryptomonas sp. and Plagioselmis lacustris, and in LN-A2+, the green algae species Botryococcus braunii were the dominant species. In HN, Cyanobacteria were the dominant group throughout the study in all temperature treatments, with biomass contributions ranging from 35% to 75%, except for HN-A0, where Chlorophyta was the dominant group and constituted 48% and 49% of the total biomass before and during the HW, respectively (

Nutrient Effects
The LMM analyses revealed that high-nutrient conditions, regardless of temperature treatments and periods, increased total phytoplankton, Chlorophyta, Bacillariophyta, and Cyanobacteria biomasses but decreased the grazing pressure and genus richness compared to low nutrient conditions ( Table 2, Figures 3 and 4). Nutrient enrichment did not change RUE and Chrysophyta biomass before the HW, whereas it led to a significant increase in these variables during and after the HW and increased Cryptophyta biomass during the HW. In contrast, effects on Dinophyta and genus evenness were not significant in any of the periods.

Nutrient Effects
The LMM analyses revealed that high-nutrient conditions, regardless of temperature treatments and periods, increased total phytoplankton, Chlorophyta, Bacillariophyta, and Cyanobacteria biomasses Water 2020, 12, 3394 8 of 21 but decreased the grazing pressure and genus richness compared to low nutrient conditions ( Table 2, Figures 3 and 4). Nutrient enrichment did not change RUE and Chrysophyta biomass before the HW, whereas it led to a significant increase in these variables during and after the HW and increased Cryptophyta biomass during the HW. In contrast, effects on Dinophyta and genus evenness were not significant in any of the periods.
the HW (Table 2), excluding the Cryptophyta biomass, which was higher in HN-A2+ than HN-A0 and the interaction effect between nutrient and A2+ warming increased its biomass (Figure 4). With the HW effect, A2 and A2+ warming started to show significant changes in LN (Table 2). During the HW in LN-A2, genus richness was significantly lower, Chrysophyta biomass was significantly higher, and cyanobacteria tended to be higher than in LN-A0. During the HW in LN-A2+, there was a significant higher biomass of total phytoplankton, Chlorophyta, Chrysophyta, and higher RUE, while the grazing pressure was significantly lower than LN-A0. In addition, Dinophyta and cyanobateria tended to be higher during the HW in LN-A2+ than LN-A0 (Table 2, Figures 3 and 4). During the HW, no significant differences were observed for HN-A2 and HN-A2+ compared to HN-A0, except a tendency of raise in cyanobacteria in HN-A2. The interaction effect between nutrient and A2+ warming was positive for grazing effect and negative for Chlorophyta biomass and RUE during the HW.   Table 2. p values of LMMs (bold values indicate p < 0.05, heat wave period (during) is highlighted, (*) refers to interaction models and NA refers to higher AIC of models with interactions. Note that the 'Nutrient' rows include results from comparison of HN-A0 with LN-A0 in all three periods not exposed to HW, (+) and (−) were used for higher and lower effects than in A0, respectively).

Crypto-Phyta
Chryso-Phyta  RUE were lower in both HN-A2 and HN-A2+ than HN-A0, Dinophyta tended to be lower in HN-A2, and the grazing effect was significantly higher in HN-A2+ than HN-A0 after the HW. The interaction effects of nutrient and A2 warming was negative for the biomass of Cryptophyta and Dinophyta and for RUE, while it was positive for the Chrysophyta biomass. The interaction between nutrient and A2+ warming was negative for Total phytoplankton, Chlorophyta, Bacillariophyta biomasses, and RUE and it was positive for the grazing pressure after the HW (  Figures 3 and 4). Ordination plots of COI analyses are given in Figure 4. Panels (a) and (d) show the plots for explanatory variables and panels (b) and (e) the plots for response variables in LN and HN, respectively. COI analyses found significant relationships (p < 0.001) between response and explanatory variables for both LN and HN. RV coefficients were higher in LN (0.44) than in HN (0.33). The variance explained by the first two axes were also higher in LN (91.4% and 6.0% for factor 1 and 2, respectively) than in HN (78.0% and 16.4% for factor 1 and 2, respectively). Panels (c) and (f) display the concordance between the explanatory and response variable panels including the datasets 3.4. Temperature Treatments: Before, during, and after the Heat Wave A2 and A2+ warming in both LN and HN did not differ from the control treatments (A0) before the HW (Table 2), excluding the Cryptophyta biomass, which was higher in HN-A2+ than HN-A0 and the interaction effect between nutrient and A2+ warming increased its biomass (Figure 4). With the HW effect, A2 and A2+ warming started to show significant changes in LN (Table 2). During the HW in LN-A2, genus richness was significantly lower, Chrysophyta biomass was significantly higher, and cyanobacteria tended to be higher than in LN-A0. During the HW in LN-A2+, there was a significant higher biomass of total phytoplankton, Chlorophyta, Chrysophyta, and higher RUE, while the grazing pressure was significantly lower than LN-A0. In addition, Dinophyta and cyanobateria tended to be higher during the HW in LN-A2+ than LN-A0 (Table 2, Figures 3 and 4). During the HW, no significant differences were observed for HN-A2 and HN-A2+ compared to HN-A0, except a tendency of raise in cyanobacteria in HN-A2. The interaction effect between nutrient and A2+ warming was positive for grazing effect and negative for Chlorophyta biomass and RUE during the HW.
After the HW, we continued to observe most changes in LN (Table 2). Genus richness was significantly lower and Cryptophyta biomass was higher in LN-A2 than LN-A0, while Cyanobacteria biomass tended to be higher. In LN-A2+, total phytoplankton, Chlorophyta biomass and RUE were significantly higher and grazing effect was lower than LN-A0 after the HW (Figures 3 and 4). In this period, A2 and A2+ warming caused significant changes in HN. While Bacillariophyta biomass and RUE were lower in both HN-A2 and HN-A2+ than HN-A0, Dinophyta tended to be lower in HN-A2, and the grazing effect was significantly higher in HN-A2+ than HN-A0 after the HW. The interaction effects of nutrient and A2 warming was negative for the biomass of Cryptophyta and Dinophyta and for RUE, while it was positive for the Chrysophyta biomass. The interaction between nutrient and A2+ warming was negative for Total phytoplankton, Chlorophyta, Bacillariophyta biomasses, and RUE and it was positive for the grazing pressure after the HW ( Table 2, Figures 3 and 4).
Ordination plots of COI analyses are given in Figure 4. Panels (a) and (d) show the plots for explanatory variables and panels (b) and (e) the plots for response variables in LN and HN, respectively. COI analyses found significant relationships (p < 0.001) between response and explanatory variables for both LN and HN. RV coefficients were higher in LN (0.44) than in HN (0.33). The variance explained by the first two axes were also higher in LN (91.4% and 6.0% for factor 1 and 2, respectively) than in HN (78.0% and 16.4% for factor 1 and 2, respectively). Panels (c) and (f) display the concordance between the explanatory and response variable panels including the datasets representing before, during, and after the HW. The grouping of the data points (locations) during the HW was not distinguishable from the other periods for either LN or HN. However, the locations of A0, A2, and A2+ in the before period together with A0 in the during and after periods of the HW (treatments not exposed to HW) differ from the locations of A2 and A2+ of during and after periods for both LN and HN (Figure 5c,f).
Water 2020, 12, x FOR PEER REVIEW  12 of 21 representing before, during, and after the HW. The grouping of the data points (locations) during the HW was not distinguishable from the other periods for either LN or HN. However, the locations of A0, A2, and A2+ in the before period together with A0 in the during and after periods of the HW (treatments not exposed to HW) differ from the locations of A2 and A2+ of during and after periods for both LN and HN (Figure 5c,f). According to the (a) and (b) plots in Figure 5, Cyanobacteria biomass was strongly related to temperature in LN. Total phytoplankton and Chlorophyta biomass and RUE were related positively to TP and TN, and negatively to Cladocera. In LN, genus evenness was correlated positively with Copepoda and PVI, and negatively with genus richness. In HN, Cyanobacteria were strongly correlated with TP. Total phytoplankton biomass was strongly positively correlated with TN and negatively with PVI. Copepoda was strongly negatively related to Bacillariophyta and RUE, while Cladocera affected Chlorophyta negatively, but for both, the arrows are very short. PVI affected the genus evenness positively and the genus richness was negatively related to the temperature, TP, and Rotifera biomass in HN (Figure 5d,e).

Ecosystem Stability
LN-A2 showed higher values than the pre-HW benchmark starting from the second week (third sampling) of the HW. The values became significant at the end of the HW (Figure 6a), indicating poor resistance. Similarly, LN-A2+ had higher values than pre-HW benchmark, which became significant after the second week (Figure 6b). HN-A2 showed the lowest resistance among all the treatments immediately after the start of the HW, but approached the benchmark value after the third week (forth sampling, Figure 6c). HN-A2+ was most resistant and stayed close to its benchmark during the According to the (a) and (b) plots in Figure 5, Cyanobacteria biomass was strongly related to temperature in LN. Total phytoplankton and Chlorophyta biomass and RUE were related positively to TP and TN, and negatively to Cladocera. In LN, genus evenness was correlated positively with Copepoda and PVI, and negatively with genus richness. In HN, Cyanobacteria were strongly correlated with TP. Total phytoplankton biomass was strongly positively correlated with TN and negatively with PVI. Copepoda was strongly negatively related to Bacillariophyta and RUE, while Cladocera affected Chlorophyta negatively, but for both, the arrows are very short. PVI affected the genus evenness positively and the genus richness was negatively related to the temperature, TP, and Rotifera biomass in HN (Figure 5d,e).

Ecosystem Stability
LN-A2 showed higher values than the pre-HW benchmark starting from the second week (third sampling) of the HW. The values became significant at the end of the HW (Figure 6a), indicating poor resistance. Similarly, LN-A2+ had higher values than pre-HW benchmark, which became significant after the second week (Figure 6b). HN-A2 showed the lowest resistance among all the treatments immediately after the start of the HW, but approached the benchmark value after the third week (forth sampling, Figure 6c). HN-A2+ was most resistant and stayed close to its benchmark during the entire disturbance (Figure 6d). However, its performance was higher than the pre-HW benchmark at the beginning of the HW, but not significantly so. The slope during the HW indicating the change in the resistance was positive and significant for LN, whereas it tended to be negative for HN, but insignificant (Figure 6e). Since LN-A2 and LN-A2+ showed significantly low-resistance, recovery and resilience were calculated for these two treatments for the last sampling (28.08, Figure 6f). While LN-A2 did not return to its pre-HW conditions, LN-A2+ showed recovery and reached a value close to its benchmark. Resilience after the disturbance was negative but insignificant for both treatments (Figure 6g).

Phytoplankton Size Structure
Phytoplankton size diversity did not show significant responses to the treatments throughout the study period, and negative values were rare (only 3 of the 216 samples analyzed). Size diversity revealed overall high values (the different sizes classes contributed relatively equally to the size distribution, [87]). However, size evenness decreased in both HN-A2 and HN-A2+ during and after the HW and showed interaction effects between nutrient and both A2 and A2+ warming ( Table 2). As we found significant results for size evenness but not for size diversity in HN, we investigated the size classes (SCs) in HN further to identify within-class changes (Figure 7, Table 3). During the HW, all the SCs increased significantly in HN-A2, while SC1 and SC4 increased in HN-A2+. After the HW, SC2 significantly increased in HN-A2, and SC4 increased in HN-A2+ (Figure 7, Table 3). Table 3. p values of size class (SC) LMMs within HN (bold values indicate p < 0.05, heat wave period (during) is highlighted). insignificant (Figure 6e). Since LN-A2 and LN-A2+ showed significantly low-resistance, recovery and resilience were calculated for these two treatments for the last sampling (28.08, Figure 6f). While LN-A2 did not return to its pre-HW conditions, LN-A2+ showed recovery and reached a value close to its benchmark. Resilience after the disturbance was negative but insignificant for both treatments (Figure 6g).

Phytoplankton Size Structure
Phytoplankton size diversity did not show significant responses to the treatments throughout the study period, and negative values were rare (only 3 of the 216 samples analyzed). Size diversity revealed overall high values (the different sizes classes contributed relatively equally to the size distribution, [87]). However, size evenness decreased in both HN-A2 and HN-A2+ during and after the HW and showed interaction effects between nutrient and both A2 and A2+ warming ( Table 2). As we found significant results for size evenness but not for size diversity in HN, we investigated the  Table 3). During the HW, all the SCs increased significantly in HN-A2, while SC1 and SC4 increased in HN-A2+. After the HW, SC2 significantly increased in HN-A2, and SC4 increased in HN-A2+ ( Figure 7, Table 3).

Nutrient Effect
We found that the nutrient treatment overall had a stronger effect on the phytoplankton biomass and community than both the warming treatments and the HW simulation (

Nutrient Effect
We found that the nutrient treatment overall had a stronger effect on the phytoplankton biomass and community than both the warming treatments and the HW simulation (Table 2). This greater effect of nutrients concurs with recent cross-system analyses of lakes at different scales in Europe [88], but they may in our case also be attributed to a much larger gradient in nutrients than temperature (Table 1). High nutrient availability led to a significantly higher total phytoplankton biomass, irrespective of temperature and HWs. We found a lower biomass ratio of zooplankton to phytoplankton ( Figure 3) and a lower size of grazers [67] in HN than in LN, indicating reduced grazer control of phytoplankton [89], which likely reflects a higher fish predation. Dominance of less edible cyanobacteria in HN (Figure 4) may have contributed to lower grazing, as has been seen in other studies [90][91][92].
Phytoplankton genus richness was significantly lower in HN than in LN, regardless of temperature treatments and periods. Several studies (e.g., [21,22,93]) have suggested the existence of a unimodal relationship between species richness and nutrients, which is low under oligotrophic and eutrophic conditions. In our experiment, cyanobacteria group was represented by much fewer species than the Chlorophyta group, and considering the fact that HN contained a large biomass of cyanobacteria, the observed decrease in genus richness was to be expected. A number of studies have shown a positive relationship between phytoplankton species richness and RUE [1,94,95], which has been attributed to a more efficient use of resources in more species rich communities [96]. However, despite the lower richness, we found higher RUE in HN than in LN. Besides the lower grazing pressure in HN, another reason for this might be that the dominant Microcystis spp. are able to use TP more efficiently than other phytoplankton taxa, which concurs with studies finding strong effects of cyanobacteria on RUE in eutrophic lakes [91,97,98].

Effects of Warming, Nutrient and Temperature Interactions, and Heat Waves
Since the mesocosm facility had already been running for 11 years at three contrasting temperatures before our experiment started, it is difficult to fully disentangle the actual HW effects of our treatment from those of the long-term temperature treatments or seasonal effects. Before the HW, the phytoplankton communities in the warming scenarios (A2 and A2+) did not differ markedly from the community in their controls (A0) in both LN and HN. During the HW, however, temperature effects were observed (in LN more than in HN), potentially indicating an HW effect ( Table 2). The changes that we observed during the HW continued after the HW in both LN and HN, which indicate a lasting effect of the HW. This is supported by the COI analyses (Figure 5c,f) showing a clear separation of the treatments not exposed to HW (all the treatments before the HW and A0 treatments during and after the HW) from the treatments exposed to HW (A2 and A2+ treatments during and after the HW).
It is well-established that high temperatures favor cyanobacteria; however, in our experiment, cyanobacteria biomass did not significantly increase during the HW in either LN or HN. Yet, in LN-A2 and LN-A2+, there was a tendency to an increase in cyanobacteria contribution. However, the cyanobacteria biomass in these two heated treatments were not significantly different from the biomass in LN-A0, which reflected an increase in the cyanobacteria biomass also during the HW (Figure 4), coinciding with the natural HW in July of this year (Figure 1) [67]. COI plots (Figure 5a,b) showed the strong relationship between temperature and Cyanobacteria biomass in LN, implying an HW effect on cyanobacteria at low nutrient concentrations. Although Cyanobacteria biomasses in HN ( Figure 4) were higher in HN-A2 and HN-A2+ than in HN-A0, we did not find significant effects of either warming or HW in HN. The higher biomasses of cyanobacteria in HN-A2 and HN-A2+ were observed not only during the HW but for the whole study period, indicating the key role of nutrients. In the COI plots, Cyanobacteria biomass was strongly linked to TP in HN (Figure 5d,e) and to temperature in LN, indicating that the sensitivity of cyanobacteria to high temperatures is dependent on nutrient availability.
Few other mesocosm studies have treated the effects of warming on phytoplankton at low nutrient concentrations. Verbeek et al. [57] observed an increase of both cyanobacteria and Chlorophyta with rising temperatures under low nutrient conditions, which can be expected, since they have similar optimum growth rates (25-35 • C, [99]). Hence, other environmental factors (e.g., nutrients, grazing pressure, mixing frequency, adaptation features to stratified conditions) that determine whether cyanobacteria or Chlorophyta become dominant in hot summers [100]. In our experiment, the biomass of Chlorophyta increased significantly in LN-A2+ during the HW, and it was dominated by Botryococcus braunii, which has been shown to have the highest growth at 30 • C [101]. By forming big colonies, this species might be strongly resistant to grazing, which could explain the decline in the zooplankton:phytoplankton biomass ratio in LN-A2+ (Table 2, Figure 3). The zooplankton:phytoplankton biomass ratio was higher in LN than in HN, especially in LN-A0 and LN-A2 (Figure 3). Using the same experiment, Işkın et al. [67] found higher Cladocera and Copepoda biomasses in these treatments. As the COI plots revealed, PVI affected Copepoda and Cladocera biomass positively, coinciding with a decrease in phytoplankton genus richness and an increase in evenness (Figure 5b). We found a significant lower genus richness in LN-A2 during and after the HW, which potentially could be attributed to a higher grazing effect in LN-A2, where zooplankton biomass was probably higher due to a higher plant abundance [67].
In general, the HW only caused minor changes in HN. Even though we did not observe any significant change in HN-A2+, there were negative interaction effects between nutrient and A2+ warming on Chlorophyta biomass, RUE, and the grazing pressure during the HW. Together with the HW, nutrient and A2+ warming reversed the individual effects of nutrient and A2+ warming in HN-A2+. The robust change after the HW was a decrease in the Bacillariophyta biomass, which seemed linked to the enhanced grazing pressure as evidenced both from the COI plots (Figure 5d,e) and zooplankton data [67]. Işkın et al. [67] showed an increase of large bodied Copepoda and Cladocera in HN-A2+ after the HW following fish kills, and the COI plots revealed opposite relationships between Cladocera and Chlorophyta and between Copepoda and Bacillariophyta in HN.

Ecosystem Stability
Few studies exist on the effects of HW as a disturbance factor for phytoplankton community structure. We observed higher performance (than pre-disturbance conditions) in all treatments of total phytoplankton biomass during the HW (Figure 6a-d), which may be an indicator of increased biomass at high temperatures. However, the pronounced difference in the response of the LN compared with HN ecosystems was the timing. Thus, HN ecosystems responded quickly to the HW, from the very beginning of the disturbance, and recovered before the end of the HW, whereas LNs showed high resistance at the beginning of the disturbance, which did not last until the end. Overall, the LN systems were more sensitive to the HW in terms of total phytoplankton biomass, as shown in the resilience boxplots (Figure 6e), which supports our previous analyses ( Table 2). In LN-A2, recovery was not observed, which was to be expected based on the theory that more diverse ecosystems are more stable [55], genus richness being significantly lower in LN-A2 than in LN-A2+ in our experiment. However, this theory does not concur with the results from the HN mesocosms that also had low genus richness but exhibited higher resistance. As Reynolds [102] stated, at long-term disturbance, the survival chance of sensitive groups is much lower than that of adaptive groups. The HN treatments were already dominated by advantageous cyanobacteria, which is often a dominant group under high nutrient and high temperature conditions.

Phytoplankton Size Structure
Faster metabolic rates are expected at a higher temperature, which tends to favor organisms with smaller body sizes [13]. We expected a shift to smaller size in the LN heated mesocosms as a consequence of both the nutrient depletion and higher temperatures, but a significant decrease in size diversity or evenness did not occur during the HW or in LN in the A2+ treatments throughout the study. Size diversity did not change in HN either. Size evenness, however, decreased in the heated HN mesocosms (Table 2). Both the smallest (SC1) and the largest (SC4) groups were lower during the HW in both HN-A2 and HN-A2+ compared with HN-A0. While the higher SC1 values before and during the HW can be explained by Microcystis spp., diatoms with larger body sizes seemed to cause the increase in SC4, possibly as a consequence of the grazing pressure on Bacillariophyta (Figures 5 and 7).

Conclusions
We examined the response of phytoplankton communities to a 1-month heat wave nested within a long-term experimental study run at three temperature crossed with two nutrient levels. High nutrient treatments were mostly dominated by cyanobacteria, while low nutrient treatments were richer in genera and dominated by Chlorophyta throughout the study period. While nutrients had a strong impact both before, during, and after the heat wave, the temperature and heat wave effects were modest. Temperature and heat wave effects were observed mostly in LN systems, indicating that the sensitivity of the phytoplankton community structure to high temperatures depends on nutrient availability. The study revealed the increasing vulnerability of freshwater ecosystems phytoplankton community to warming as well as HWs with declining nutrient concentrations. Ecosystem stability in the high-and low-nutrient mesocosms was affected differently by the heat wave disturbance. High-nutrient tanks showing greater stability, most likely, they were already dominated by cyanobacteria. Size diversity did not shift from large to small with warming, as otherwise expected. Our results indicate the need for stronger nutrient reductions for water management in a warmer climate in order to maintain or obtain the requested status of minimum "good ecological state" according to the European Water Framework Directive.