The biotic and abiotic factors driving freshwater phytoplankton community composition across a continental scale

The biotic and abiotic factors driving freshwater phytoplankton

community composition across a continental scale

Jelmer B. Klaassen University of Amsterdam IBED-FAME

Supervision of Jolanda Verspagen DOI: 10.5281/zenodo.3904475

Abstract

Phytoplankton are a crucial group within freshwater communities, as they form the base of the food web and are important primary producers. This makes them vital to the ecosystem functioning in lakes. However, certain phytoplankton groups can have adverse effects, such as species who are toxic to other organisms or humans. On those grounds, it is of societal and scientific importance to understand the abiotic and biotic factors that drive the abundance of different phytoplankton groups. Nevertheless, recent research mostly focused on toxic or abundant species and genera in small geographical ranges, while this study focused on environmental factors that drive a wider range of phytoplankton groups across a continental scale to supplement existing knowledge. Using datasets from 1235 U.S. lakes, the biomass of the 7 most important taxonomic phytoplankton groups, the biomass of phytoplankton groups with 4 specific traits and the phytoplankton genus richness are statistically analyzed to examine which of 19 ecological relevant factors explain the abundance of different freshwater phytoplankton groups best. The results show that there are large differences in the predicting factors among phytoplankton groups. Turbidity, Secchi-depth and nutrients are generally significant drivers in the models that explain the biomass of total phytoplankton and all of the taxonomic groups, except for chrysophyte biomass. Chrysophyte biomass is largely driven by temperature and Secchi-depth, while temperature, together with chlorophyll a and zooplankton abundance, is also a significant driver of genus richness. Trait-based phytoplankton biomasses are mostly explained by nutrients, depth, lake area and CO2 concentration, although differences exist between traits. In conclusion, new drivers that affect certain freshwater phytoplankton groups were identified, but most of these results confirm the findings in previous research who analyzed a small geographical scale, suggesting that identified drivers on a small geographical scale can be extrapolated. These various results can help water management, conservational management and freshwater ecologists to predict the effects of environmental changes on freshwater phytoplankton communities across a large geographical scale.

Introduction

Lakes are an important source of recreation and drinking water, and also vital to fisheries and the economy. An abundant group of crucial primary producers in lakes is phytoplankton. These microscopic plants float in the water column and they form the base of the food web (Reynolds, 1984), which makes their presence crucial to zooplankton, fish(eries) and humans. However, toxic blooms of some of these algae can have adverse effects on food webs, water quality and human health (Reynolds, 1984; Paerl & Huisman, 2009). Thus, it is of societal, ecological and economical importance to understand the ecology of different phytoplankton groups. However, to date, most research has focused on which factors drive these harmful cyanobacterial blooms, while the drivers of abundances of other taxonomic or functional phytoplankton groups are typically overlooked. Few studies that do include other groups restrict their analysis to limiting geographical regions. This study aims to embody a more comprehensive approach, in order to identify the factors that drive the abundance of different freshwater phytoplankton groups across a broad geographical scale.

Factors that are believed to be the most important drivers of total phytoplankton biomass and community composition are temperature and nutrients (Boynton, Kemp, & Keefe, 1982; Jensen, Jeppesen, Olrik, & Kristensen, 1994), but this may differ among groups. Important taxa of freshwater phytoplankton in temperate lakes are diatoms, cyanobacteria, chlorophytes, dinoflagellates, euglenophytes, cryptophytes and chrysophytes (Watson, McCauley, & Downing, 1997). With a growing interest in using a trait-based approach to better understand the (future) composition of phytoplankton communities (Litchman & Klausmeier, 2008), this study will also analyze fifty of the most abundant genera of freshwater phytoplankton on four traits to assess the functional response to environmental factors. These traits that define phytoplankton contributions to ecosystem functioning are the capacities for nitrogen fixation, mixotrophy, silica-requirement, and the ability to form colonies (Vogt, St‐Gelais, Bogard, Beisner, & del Giorgio, 2017). First off, the different phytoplankton groups and the possible drivers of their abundance will be briefly discussed.

Taxonomic groups

The most abundant group of freshwater phytoplankton are the cyanobacteria (blue-green bacteria), whose fossils are dated up to 3.5 billion years ago, making them the oldest known oxygen-producing organisms on Earth (Whitton & Potts, 2007). Due to their long evolutionary lifespan, cyanobacteria have evolved diverse adaptations to survive climatic and human-induced changes (Huisman, Matthijs, & Visser, 2005). Eco-physiological adaptations include the ability to grow in warmer climates, gas vesicles to float up to the water surface, high affinity for phosphorus and the unique ability in some genera to fix nitrogen (Carey, Ibelings, Hoffmann, Hamilton, & Brookes, 2012). These adaptations combined might explain why nutrients (nitrogen and phosphorus), pH and temperature are often found to correlate with freshwater cyanobacterial biomass (Beaulieu, Pick, & Gregory-Eaves, 2013; Rigosi, Carey, Ibelings, & Brookes, 2014; Cremona, Tuvikene, Haberman, Nõges, & Nõges, 2018). With increasing temperatures and eutrophication, it is thus of scientific and societal importance to find out the relations between environmental factors and the abundance of cyanobacteria to help water management to prevent ecologically harmful cyanobacterial dominance, which is already expanding (Kosten et al., 2012).

Another important group of phytoplankton is the class Bacillariophyceae (diatoms), which is a widespread group of 12.000 described freshwater phytoplankton species that are abundant and globally important in the carbon and silica cycle (Mann, 1999; Smetacek, 1999; Guiry, 2012). Through photosynthesis, diatoms are responsible for roughly 40% of the CO2 uptake by the ocean (Nelson, Tréguer, Brzezinski, Leynaert, & Quéguiner. 1995). The distinguishing feature of all diatoms is a silica-containing cell wall, the frustule (Round, Crawford & Mann, 2007). This makes them dependent on silica concentrations in their environment (Schöllhorn & Granéli, 1996), as silica-depletion can lead to a chlorophyte and cyanobacteria dominated community (Schelske & Stoermer, 1971). However, most research is done on marine diatoms. It is thus of scientific importance to identify the factors that steer freshwater diatom abundance.

Green algae consist of approximately 8000 mostly freshwater phytoplankton species of the phylum Chlorophyta and the group of charophytes (Hoek, Mann, Jahns, & Jahns, 1995; Leliaert et al., 2012), united by the presence of the accessory pigment chlorophyll b. Green algae have been ecologically important for hundreds of millions

of years (Falkowski et al., 2004) and are globally distributed (Hoek et al., 1995). However, previous research mainly focused on a few (marine) species, but the effect of potential drivers on the whole group of freshwater green algae is poorly understood.

Like the cyanobacteria, bloom-forming species of dinoflagellates, euglenophytes and chrysophytes have been the main focus point of previous research, while a holistic overview for these taxonomic groups was generally overlooked. The phylum Dinoflagellata consists of about 2500 extant protists, of which approximately 350 are bloom-forming freshwater species (Sandgren, 1988; Guiry, 2012). The phylum Euglenophyta contains more than 800 mostly freshwater species of unicellular biflagellates with a substantial part of them being heterotrophic (Hoek et al., 1995). They can tolerate a wide pH range (2-8) and elevated metal concentrations (Olaveson & Nalewajko, 2000), which enables them to withstand a wide range of often polluted habitats, from acid and alkaline waters to aerobic and anoxic waters, where they can form colored blooms which can be medically beneficial to humans (Barsanti, Vismara, Passarelli, & Gualtieri, 2001; John, Whitton, & Brook, 2011). Golden algae, or chrysophytes, of the class Chrysophycea comprise 761 known species that are mostly found in freshwater habitats (Hibberd, 1976; Guiry & Guiry, 2020) and some chrysophytes species need silica to grow (Kristiansen, 1986). These flagellates can produce toxins that mass kill fish populations (Shilo, 1971). Thus, a better understanding of the drivers that shape their abundance can help water managers to prevent mass dying of fish and consequently creating a healthier ecosystem.

The group of cryptophyta (or cryptomonads) contain 100 described freshwater lake species that share the presence of complex and highly characteristic flagellates (Hoek et al., 1995; Guiry, 2012). They are a substantial part of the phytoplankton community and are particularly important and abundant in deep layer populations around the chemocline and deep chlorophyll maximum (DCM) (Gervais, 1998). As a whole, the factors that drive the abundance of the group of cryptophytes is poorly understood, with previous research focusing on few species (Knisely & Geller, 1986; Lewitus & Caron, 1990; Gasol, García-Cantizano, Massana, Guerrero, & Pedrós-Alió, 1993; Gervais, 1997; Gervais, 1998).

Trait-based groups and genus richness

Besides taxonomic groups, this research will focus on trait-based ecology. An important trait for phytoplankton is mixotrophy. Mixotrophic phytoplankton are able to take up organic matter through phagotrophy, besides their ability to photosynthesize (Arenovski, 1994). Mixotrophy is mostly found in the phytoflagellates and is an important part of a planktonic food web. It is also important to understand the factors behind the abundance of mixotrophic algae as species that form harmful algal blooms are often mixotrophic (Anderson, Glibert, & Burkholder 2002). Most of this research is done in boreal lakes, but the factors that drive the abundance of mixotrophic phytoplankton on a large geographical scale is poorly understood.

Diatoms and some chrysophytes need silica to grow (Hoek et al., 1995). Thus, a relation between silica concentration and the abundance of the trait-based group of Si-requiring phytoplankton is expected (Egge & Aksnes, 1992), while phosphorus limitation and zooplankton grazing lowers the dominance of diatoms at high silica concentrations (Castenholz, 1961; Feminella & Hawkins, 1995; Egge, 1998; Lange, Liess, Piggott, Townsend, & Matthaei, 2011). Nonetheless, there is little trait-based research on the factors that drive the whole freshwater group of Si-requiring phytoplankton.

Some cyanobacteria have the unique capacity to convert molecular nitrogen into inorganic nitrogenous compounds, using the enzyme nitrogenase to convert molecular nitrogen (Berman-Frank, Lundgren, & Falkowski, 2003). Therefore, nitrogen is shown to limit the growth of N-fixing phytoplankton (Vanderhoef, Huang, Musil, & Williams, 1974), but phosphorus stimulates the growth of diazotrophic algae (Vanderhoef et al., 1974; Moisander, Steppe, Hall, Kuparinen, & Paerl, 2003; Vogt et al., 2017). This leads to a problem for water managers who aim to control the effects of eutrophication by reducing the input of nitrogen, as N-fixing phytoplankton will increase nitrogen concentrations again under phosphorus-sufficient conditions (Schindler et al., 2008). It is therefore important to understand the factors that determine the abundance of N-fixing phytoplankton. Many phytoplankton taxa can form colonies in different shapes or forms (Reynolds, 1984). One of the advantages of this trait is that this is a defense mechanism against grazing (Reynolds, 1984), so large zooplankton densities

might favor the (relative) abundance of colony-forming species. It has been shown that nutrients positively affect the biomass of colony-forming phytoplankton groups, dependent on the taxa (Naselli-Flores & Barone, 2000). Colony-formation as a defense mechanism that is induced by zooplankton grazing is studied extensively (Van Donk, Ianora, & Vos, 2011), but overall there is little knowledge on colony-formation as a trait in freshwater community compositions and the other factors that drive the abundances of this functional group.

Lastly, the diversity of phytoplankton is crucial to the resource use efficiency and productivity of an aquatic ecosystem (Ptacnik et al., 2008). Loss of diversity could be detrimental to primary production of a lake, so understanding the factors which drive the loss of genus richness could help management. Temperature, chlorophyll and lake area were found to correlate to freshwater genus richness on a large geographical scale (Stomp, Huisman, Mittelbach, Litchman, & Klausmeier, 2011), but similar research is scarce.

This study will test the drivers of different phytoplankton groups by analyzing the relationships of 19 environmental and biological factors as possible predictors and the abundances, biomasses and genus richness of phytoplankton groups as response variables (R Core Team, 2013). Using Random Forest analysis and multi-model interference, the importance of individual predictors and their interactive effect are tested and combined in a predictive model, following the ‘cookbook’ by Feld, Segurado, & Gutiérrez-Cánovas (2016). This kind of research has not been previously done on such a large and broad scale. Previous research analyzed smaller datasets across smaller geographical ranges with less taxonomic or functional groups. This study combines knowledge from previous research to find a holistic answer to the question which factors drive the abundance of different taxonomic and functional groups of freshwater phytoplankton. The results of this broad study could add to the existing knowledge about freshwater phytoplankton habitats and could help water management, conservational management and phycology research as anthropogenic stressors are becoming a larger problem.

Methods

US-EPA Dataset

This research uses a large-scale dataset to answer the question which factors influence the biomass of different phytoplankton groups. In 2012, the US Environmental Protection Agency (EPA) conducted a survey at >1000 lakes to monitor the conditions of American lakes (U.S. Environmental Protection Agency, 2016). The data were made public and can be downloaded online ( https://www.epa.gov/national-aquatic-resource-surveys/data-national-aquatic-resource-surveys). The extensive methods of this survey can be found in USEPA (2017). In short, field teams sampled in randomly selected freshwater bodies which were greater than 1 hectare, were at least 1 meter deep, and had a minimum of 0.1 hectare open water. Water samples to analyze phytoplankton, chlorophyll a, nutrients (e.g. total nitrogen and total phosphorus) and other water chemistry variables were collected in the euphotic zone with an integrated water sampler, which is 2 m long and has an inside diameter of 3.2 centimeters. The survey used multi-probe sondes to measure temperature, oxygen and pH at depth intervals, while a Secchi-disk measured the transparency of a lake. Zooplankton was collected with two

Wisconsin nets, one fine mesh (50 µm) and one rough mesh (150 µm). This research used the zooplankton data of the fine mesh to get more inclusive data. After this collection in the field, phytoplankton and zooplankton taxa were classified by taxonomists microscopically and the biovolume (μm3_{/mL) of a taxon was calculated as} their abundance (cells/mL) times their biovolume of a cell (μm3_{). Lastly, water samples were analyzed in the lab} on different chemical variables and nutrient concentration.

The dataset contains values of many physical, chemical and biological parameters, including phytoplankton biomasses. 20 of those abiotic variables (Table 1) were selected for this analysis because they are indicative for important ecological processes, such as eutrophication (nutrient concentrations), climate change (temperature or pCO2) or biological relations (zooplankton abundance or chlorophyll a). These factors were analyzed on their effect on the phytoplankton biomass (volume of cells/mL) of 7 important taxonomic groups. Furthermore, the 50 freshwater phytoplankton genera with the most total biomass in this dataset were classified based on previous literature on 4 traits proposed by Vogt et al. (2017): N-fixation, mixotrophy, Si-requirement and colony formation (Supplementary Table 1). The share of the total biomass with that trait will be used as a response variable to analyze the effect of the predictors. Lastly, this research will also look at factors that drive

genus richness per site. Analysis of these data will be executed by performing a multiple stressor analysis in Rstudio (R Core Team, 2013), following the step-wise procedure of Feld et al. (2016).

Table 1. The mean, standard deviation, minimum and maximum values of 19 physical, chemical and biological characteristics of 1235 American lakes. All predictor variables are measured by the EPA survey, except for CO2. CO2 concentrations are calculated according to Verspagen et al. (2014), using temperature, pH, acid neutralizing capacity and conductivity.

Predictors Mean (± SD) Min Max N

Physical Depth (m) 3.16 (± 3.56) 0.80 24.07 1223

Temperature (⁰C) 21.93 (± 10.95) 6.81 35.50 1204

Secchi-depth (m) 2.13 (± 2.31) 0.02 28.00 1062

Turbidity (NTU) 9.45 (± 27.28) 0.01 447.15 1185

Latitude (⁰) - 26.07 48.98 1234

Lake Area (ha) 1110.84 (± 9611.65) 1.03 167489.61 1138

Chemical Total phosphorus (μg/L) 114.75 (± 269.63) 4 3,636 1225

Total nitrogen (mg/l) 1.12 (±2.49) 0.01 54.00 1227 Silica (mg/L) 10.88 (± 29.71) 0.02 935.00 1227 Calcium (mg/L) 28.33 (± 45.78) 0.12 594.90 1229 pH 7.79 (± 0.99) 3.25 11.87 1199 CO2 (ppm) 2,599.87 (± 5,562.06) 0.03 79,910.66 1165 Oxygen (mg/L) 6.78 (± 2.55) 0.25 31.80 1193

Dissolved Organic Carbon (mg/L) 8.29 (± 18.54) 0.23 515.81 1227 Conductivity (μS/cm) 688.36 (± 3,161.70) 2.82 64,810.00 1074 Acid neutralizing capacity (μeq/L) 2,145.05 (± 2,694.44) -3,361.40 29,112.00 1225

Biological Chl a (μg/L) 26.12 (± 54.06) 0.00 764.64 1225

Zooplankton biomass (µg dry weight/L) 280.54 (± 557.55) 0.01 6,205.29 1234 Zooplankton density

(number of organisms/L)

1,647.29 (± 15,711.55)

0.20 542,177.05 1233

Phytoplankton Total phytoplankton biovolume (µm3_{/mL) 1.7E+07 (± 4.2E+07) 14.989 5.8E+08 1235} Cyanobacteria biovolume (µm3_{/mL) 7.6E+06 (± 2.7E+07) 0.00 4.3E+08 1217} Chlorophyte biovolume (µm3_{/mL) 3.0E+06 (± 1.1E+07) 0.00 1.3E+08 1207} Diatom biovolume (µm3_{/mL) 2.5E+06 (± 1.1E+07) 0.00 2.0E+08 1164} Cryptophyte biovolume (µm3_{/mL) 1.1E+06 (± 5.3E+06) 0.00 1.3E+08 1127} Euglenophyte biovolume (µm3_{/mL) 1.5E+06 (± 1.4E+07) 0.00 4.0E+08 796} Dinoflagellate biovolume (µm3_{/mL) 3.6E+06 (± 2.1E+07) 0.00 4.4E+08 743} Chrysophyte biovolume (µm3_{/mL) 1.1E+06 (± 1.1E+07) 0.00 2.3E+08 538} Mixotrophic % of the biovolume (% µm3_{/mL) 40.48 (31.87) 0.00 100 1235} Si-requiring % of the biovolume (% µm3_{/mL) 15.34 (23.34) 0.00 99.39 1235} N-fixing % of the biovolume (% µm3_{/mL) 18.78 (27.38) 0.00 100 1235} Colony-forming % of the biovolume (% µm3_{/mL) 53.20 (30.59) 0.08 100 1235}

Genus richness 29.58 (11.84) 2 69 1235

Statistical analysis

First off, every variable is screened for outliers, using the function boxplot(). Values are classified as outliers when the observation passes more than 1.5 times the range between the 25th and 75th percentiles, called the interquartile range. Outliers that are physically impossible or are a result of poor data management are excluded of the raw data. Following outlier management, variables are log10_{-transformed if necessary to} approach normality. Shapiro-Wilk tests are performed on all variables to check statistically for normality. As all the variables have different scaling, all variables are standardized using z-transformation, which makes it possible to compare effect sizes between variables. Standardization followed a log10_{-transformation, if} relevant.

Furthermore this analysis assumes and requires that factors are independent, which means that they do not correlate with one another. This can be quantified with Variance Inflation Factors (VIFs), using the function vifstep() from the package usdm (Naimi, 2015). This function automatically excludes variables with a high VIF until all remaining variables show a VIF < 8, while it also accounts for non-linear relationships. Factors that are expected beforehand to explain the biomass of certain phytoplankton groups are not excluded, even if they have a VIF > 8. Correlations between the most influential factors are visualized with a correlation matrix, using the function chart.Correlation() from the package PerformanceAnalytics (Peterson et al., 2018).

As only the most influential factors and their interactions are of interest for the quantification of their effects, collinearity analysis and stressor hierarchy analysis (Random Forest) indicate which variables to include in the Generalized Linear Models. The package randomForestSRC (Ishwaran & Kogalur, 2016) can fit a large number of models on bootstrapped subsets of the data, after which factors can be ranked on their effect size per response variable. These rankings are visualized by creating bar graphs. Using the function rsfrc() of the package, the analysis is performed while accounting for non-linear relationships, complex interactions and missing values, which makes it particularly suitable for this research.

In order to find the model that predicts the most variance with the least amount of factors, multi-model inference with the R package MuMIn (Barton, 2009) selects the final model for testing multiple-stressor effects. All the variables and interactions that were hypothesized a priori or had a high variables importance (VIMP) resulting from the RF analysis are used as input for the function dredge(). The parameter Akaike’s Information Criterion weight (AICw) indicates the model with the highest explanatory power. If there is no extensive evidence for a single model after the first run, the predicting factor with the least explaining power is left out of the function dredge(). This repeats until any loss of variables results in a statistically less powerful model, a process also called backwards elimination. If an interaction is part of the model, the inclusion of the independent factors is obligatory, even if they are insignificant.

The most parsimonious model per response variable that is suggested after multi-model interference is analyzed on its explaining power with a Generalized Linear Model (GLM), using the function glm() in Rstudio. A statistically valid model requires residuals to display normality and homoscedasticity, which was checked visually.

Results

After exclusion of outliers, all explaining variables except O2, pH, temperature and Secchi were

log10-transformed. The biomass of taxonomic groups were also log10-transformed, while the share of biomass of the trait-based groups and genus richness did not need transformation (Supplementary Table 2). VIF analysis then indicated that pH correlated too strongly (VIF > 8) to other factors and must be excluded from further analysis (Supplementary Table 3). However, as pH was input for calculating pCO2 and was hypothesized to play a role in driving certain phytoplankton biomasses, pH was not left out in certain models (Table 2 and Supplementary Tables 4, 12, 24 and 26). The correlation matrix with the most important factors, according to RF analysis (Table 3), indicated significant relations between numerous factors and showed normal distributions of the variables (Figure 1).

Figure 1. Correlation matrix of the most important predictors. Left plots show the relationships with the regression line in red. The right charts report the linear correlation coefficients. The diagonal bar graphs indicate the distribution of the following variables: logNTL, logTurbidity, Secchi, logDepth, logpCO2, temperature and logArea. *** < 0.001, ** < 0.01, * < 0.05 and . < 0.1.

The 7 taxonomic groups that are used in this research made up 99.28% of the total mass (Figure 2). Cyanobacteria are the most abundant group with 42.78%, while chrysophytes only made up 2.75% of total phytoplankton abundance. The 50 genera with the most biovolume that are used in the trait-based analysis made up 93,34% of the total biovolume (Supplementary Table 28).

Figure 2. Average distribution of the biovolume (µm3_{/mL) of different taxonomic phytoplankton groups,} sampled from 1235 U.S. lakes.

After standardization of all the variables and responses, a RF analysis was conducted for every response variable. Total phytoplankton biomass, Chl a and latitude were not taken into account in all models except for genus richness, as they are expected to be either consequences of biomass or are not relevant for

phytoplankton biomass analysis. Results from these RF analysis showed the importance of factors and their interactions. As an example, Table 2 displays the outcome of the RF analysis on total phytoplankton biomass. Factors are ranked on their variable importance (VIMP), which illustrates the power of explaining the response variable. In the case of total phytoplankton biomass, logTurbidity, logNTL and Secchi are the most important variables, while temperature and logZooplanktonBiomass explain relatively little.

Table 2. Importance of the 17 included variables that explain total phytoplankton biomass, which RF analysis ranked on their VIMP. The percentage VIMP per variable of the total VIMP is also calculated.

Variable: VIMP VIMP (% of Total)

logTurbidity 0.1746 24.3617 logNTL 0.1399 19.52 Secchi 0.0771 10.7576 logPTL 0.0478 6.6695 logDOC 0.0328 4.5765 logConductivity 0.0320 4.4649 logDepth 0.0258 3.5998 logZooplankton density 0.0246 3.4324 pH 0.0230 3.2092 Oxygen 0.0230 3.2092 logANC 0.0216 3.0138 logCalcium 0.0195 2.7208 logArea 0.0176 2.4557 logSilica 0.0174 2.4278 logpCO2 0.0152 2.1208 Temperature 0.0127 1.7720 logZooplankton biomass 0.0121 1.6883

CYANOBACTERIA CHLOROPHYTES DIATOMS DINOFLAGELLATES

Another outcome of the RF analysis were partial dependence plots (Supplementary Figures 5 to 16). These graphs plot the independent relationship of every variable with the response. For instance, the partial dependence plot on the factors that correlate with total phytoplankton biomass (Figure 3) shows the strong positive correlation of logTurbidity and logNTL with the total biomass and the strong negative correlation between the Secchi-depth and total phytoplankton biomass. Often, a high VIMP corresponds to a strong relation with the response variable in the partial dependence plots. The full results of every RF analysis can be found in the Supplementary, while only the meaningful results will be highlighted here.

Figure 3. Partial dependence plots from RF analysis show the independent effect of every factor on total phytoplankton biomass. The Y-axis indicates the standardized values of the log10_{-transformed total} phytoplankton biomass and the X-axis displays the standardized values of every factor.

In 8 of the 13 RF analyses, logTurbidity was the most important factor (Table 3). Also nutrients (NTL & log-PTL) and Secchi-depth were shown to be prominent drivers of many phytoplankton groups, while oxygen, temperature, log-silica, log-lake area and log-ANC generally were generally low-ranking variables. Figure 4 shows the importance of groups of key ecological variables: light, nutrients, carbon chemistry, lake

morphology, zooplankton, temperature and silica. Light (logTurbidity and Secchi) and nutrients are shown to be the main drivers of phytoplankton biomasses, although there are differences between groups. Euglenophyte (VIMP % = 35.74) and total phytoplankton biomass (VIMP % = 35.12) are relatively better explained by light, while light (VIMP % = 5.39) and carbon chemistry (logpCO2, logANC and sometimes pH) (VIMP % = 4,91) appear to hardly drive the biomass of Si-requiring phytoplankton. The biomass of Si-requiring phytoplankton is

relatively largely explained by lake morphology (logArea and logDepth) (VIMP % = 24,25) and zooplankton (VIMP % = 12,4128). Temperature and silica generally explain little of the response, although temperature is relatively important in explaining genus richness (VIMP % = 13.78).

Table 3. RF analysis ranked factors on their importance (VIMP) in explaining the response variable. Full results can be found in the Supplementary. Chloro = chlorophyte biomass; Chryso = chrysophyte biomass; Colony = colony-forming share of biomass; Crypto = cryptophyte biomass; Cyano = cyanobacterial biomass; Diatom = diatom biomass; Dino = dinoflagellate biomass; Eugleno = euglenophyte biomass; Mixo = mixotrophic share of biomass; N-fix = nitrogen-fixing share of biomass; Si-req = silica-requiring share of biomass; Temp = temperature; Total = total phytoplankton biomass; logTurb = log10_-turbidity;

VIMP Rank

Total Cyano Diatom Dino Eugleno Chryso Chloro Crypto Mixo N-fix Si-req Colony Genus Richness 1 logTurb logTurb logTurb logTurb logTurb Secchi logNTL logTurb logTurb logPTL logDepth logNTL logTurb 2 logNTL logNTL logCalcium logNTL Secchi logTurb logTurb Secchi logNTL logTurb logDOC pH Temp 3 Secchi logPTL Secchi logPTL logPTL Temp logpCO2 logNTL logpCO2 logNTL logNTL logTurb logPTL

Figure 4. RF analysis returned the variable importance (VIMP) for all the factors per response variable. In order to get a comprehensive view of the key results, only 7 key variables are visually displayed here: light, nutrients, carbon chemistry, lake morphology, zooplankton, temperature and silica. The factor light is calculated by adding up the relative percentages VIMP (Variable VIMP/Total VIMP*100) of turbidity and Secchi, which are both parameters for light availability. This is also done for nutrients with NTL, PTL and calcium and for carbon chemistry with pCO2, ANC and in some models pH. The summed relative VIMP percentages of the biomass and density of zooplankton is displayed as the factor zooplankton, while temperature and silica are given as the relative VIMP of their respective factor. Conductivity, calcium, latitude and DOC are not represented here. Following the RF analysis, multi-model interference and backwards elimination chose the most parsimonious model for every response variable, which can be found in the Supplementary. Table 4 displays the model that significantly explains 47.09% of the total phytoplankton biomass, using logNTL, logDepth, logTurbidity, Secchi, logDOC and the interaction between logTurbidity and Secchi as significant predictors. The full results of every model are found in the Supplementary.

Table 4. Most parsimonious general linear model to explain total phytoplankton biomass (n = 802, LOGLikelihood = -853.658, AIC = 1723.3, AICw = 0.617, Radj = 0.4709, p < 2.2e-16).

Summaries of the other models (Table 5) show that logTurbidity and logNTL are important in predicting the biomass of most taxonomic and trait-based phytoplankton groups. Secchi is identified as a significant predictor in all the taxonomic groups, whereas trait-based biomass is explained more by logpCO2, logDepth and logArea. Temperature is the most important factor in predicting genus richness, while zooplankton and logChla are also significant predictors. Genus richness significantly negatively correlates to latitude (Figure 5), although latitude was not returned as a predictor in the model

0,00 10,00 20,00 30,00 40,00 50,00 60,00 70,00 80,00 90,00 100,00 V IM P ( % o f to tal V IM P)

Explaining power of ecological important variables

Light Nutrients Carbon Chemistry Lake Morphology Zooplankton Temperature Silica

Factor: Estimate Standard Error T value P

Intercept 0.11409 0.04613 2.473 0.0136 logNTL 0.44304 0.05080 8.722 < 2e-16 logDepth 0.07820 0.03474 2.251 0.0247 logTurbidity 0.30376 0.07297 4.163 3.48e-05 Secchi -0.16259 0.09381 -1.733 0.0834 logDOC -0.09813 0.04298 -2.283 0.0227 logTurbidity:Secchi 0.15224 0.06326 2.407 0.0163

Table 5. Predicting models for log10 _{and z-transformed biomass of different phytoplankton groups or the genus} richness per site. Predictors with a ^ are not significant. The p of all models is < 2.2e-16, except for the model that predicts chrysophytes (p = 2.235e-10). Predictors are ranked here on their p-value from low to high. Full versions of the model can be found in the Supplementary.

Response N AIC Radj Predictors

Total Phytoplankton 802 1723.3 0.4709 0.44*logNTL + 0.30*logTurbidity + 0.15*logTurbidity:Secchi - 0.10*logDOC + 0.08*logDepth – 0.16*Secchi^

Cyanobacteria 774 1728.3 0.4582 0.34*logNTL + 0.29*logDepth + 0.19*pH – 0.24*Secchi + 0.25*logTurbidity + 0.13* Temperature

Diatoms 659 1656 0.2474 0.28*logTurbidity – 0.22*Secchi + 0.13*logZooBiomass + 0.13*logCalcium + 0.12*logDepth Dinoflagellates 375 964 0.2003 –0.23*Secchi + 0.29*logNTL + 0.09*logCalcium^

Euglenophytes 408 956 0.3884 0.27*logTurbidity – 0.36*Secchi + 0.22*logNTL + 0.12*Temperature – 0.10*Oxygen Chrysophytes 268 751.7 0.1536 –0.18*Temperature – 0.23*Secchi + 0.31*logTurbidity + 0.15*logZooDensity Chlorophytes 746 1949.4 0.2062 0.24*logNTL – 0.18*Secchi + 0.11*Temperature + 0.11*logZooDensity Cryptophytes 678 1769 0.2128 0.41*logNTL – 0.35*Secchi – 0.10*Temperature – 0.08*Oxygen – 0.11*logDOC

Mixotroph 802 2095 0.1814 –0.23*logTurbidity + 0.18*logpCO2 – 0.17*logArea – 0.12*logANC

N-fixing 802 2141.1 0.2028 0.29*logDepth – 0.18*logpCO2 + 0.17*logTurbidity + 0.17*logNTL + 0.16*logPTL + 0.09*logSilica + 0.09*Temperature

Si-requiring 802 2120.5 0.1144 0.21*logArea – 0.27*logNTL – 0.11*Temperature + 0.09*logZooBiomass + 0.12*logPTL – 0.07*logpCO2

Colony-forming 802 2137.2 0.1464 –0.15*logpCO2 + 0.16*logDepth + 0.13*logArea + 0.18*logNTL + 0.16*logPTL + 0.08*logZooBiomass

Genus Richness 802 2097.8 0.2103 0.24*Temperature + 0.29*logChla + 0.19*logZooDensity – 0.25*logPTL – 0.22*logTurbidity – 0.15*logZooBiomass + 0.09*logpCO2

Figure 5. Latitude on the x-axis (⁰) significantly negatively correlates to genus richness on the y-axis (genus number per site) (Radj = 0.01338, F = 17.72, df = 1232, p-value = 2.75e-05).

Discussion

Total phytoplankton biomass

Temperature and nutrients were expected to be the most important drivers of total phytoplankton biomass (Boynton et al., 1982; Jensen et al., 1994; Vogt et al., 2017). Multi-model interference suggested a model with logNTL, LogTurbidity, logDepth and the interaction between logTurbidity and Secchi as positive predictors of the total phytoplankton biomass (Table 4), while LogDOC and Secchi correlated negatively. Contrary to our hypothesis, temperature was an insignificant predictor with a low VIMP (Table 2), although the positive effect of nutrients was proven (Table 4). Under the interaction with mixing, depth is previously shown to affect the biomass of phytoplankton negatively (Jäger, Diehl, & Schmidt, 2008), although this model indicated a positive correlation (p = 0.0247) (Table 4). The negative impact of logDOC on phytoplankton biomass is in line with previous research (Carpenter, Cole, Kitchell, & Pace, 1998), as high concentrations of DOC cause shading and thus limited light. This is confirmed by the significant positive relation between DOC and turbidity

(Supplementary Figure 4).

Turbidity was an important predictor in this model and in most other analyses (Table 5), which could be explained by the negative relation between turbidity and light availability (Diehl, Berger, Ptacnik, & Wild, 2002). As turbidity is a measure for light, it was expected that phytoplankton biomass would decrease with increasing turbidity (Cloern, 1987). Surprisingly however, the partial dependent plots (Figure 3 and Supplementary Figures 5 to 15) show that turbidity generally correlated strongly with the total biomass and the biomass of taxonomic and trait-based groups. LogTurbidity is linked positively to logNTL (Supplementary Figure 1), which might explain the relation between turbidity and total phytoplankton biomass. Moreover, algal blooms can be a cause of turbidity (Järvenpää & Lindström, 2004) and another study showed that algal blooms can occur in a shallow estuary with high turbidity, if the critical depth is deeper than the water depth (Fichez, Jickells, & Edmunds, 1992).

Taxonomic groups

In all the models that predicted biomasses of taxonomic groups (Table 5), Secchi-depth negatively influenced the biomasses. The Secchi-depth is a parameter for water clarity and thus the penetration of light through water (Lee, Shang, Du, & Wei, 2018). Like expected, the Secchi-depth was low under turbid concentrations (Supplementary Figure 2) and the partial dependent plots (Supplementary Figures 5 to 11) demonstrated that phytoplankton biomasses decline with higher Secchi depths. This is in line with the general perception that algal biomass negatively affects the Secchi depth, because of the hyperbolic relation between Chlorophyll a and Secchi-depth (Edmondson, 1972; Canfield & Hodgson, 1983). This is also observed in this dataset

(Supplementary Figure 3). Based on these reasonings, phytoplankton therefore affects the water clarity, measured by turbidity and Secchi, and not vice versa.

Cyanobacteria were the most abundant group in the 1235 sampled lakes with 42.78% of the biomass (Figure 2). It was expected that nutrients, temperature and pH would drive cyanobacterial abundance (Beaulieu et al., 2013; Cremona et al., 2018), although a high pH is a consequence of cyanobacterial blooms (Verspagen et al., 2014). The RF analysis results indeed indicated that nutrients (logNTL and logPTL) are important factors in explaining total cyanobacterial biomass, together with light (logTurbidity and Secchi) (Supplementary Table 4). LogNTL, Secchi and logTurbidity were also included in the most parsimonious model, together with the expected drivers temperature and pH (Supplementary Table 5). The positive effect of temperature on cyanobacteria confirms the general finding that harmful cyanobacterial blooms are occurring more frequently with warmer climates (Scheffer, Rinaldi, Gragnani, Mur, & van Nes, 1997; Paerl & Huisman, 2009). Likewise, increases in eutrophication and turbidity will contribute to increasing frequencies of cyanobacterial blooms (O’neil, Davis, Burford, & Gobler, 2012), which is in agreement with the positive relationship of logNTL with cyanobacterial biomass. The inclusion of LogDepth in the model might be due to the presence of gas vacuoles in most cyanobacteria (Carey et al., 2012). This enables them to dominate in deep and weakly mixed lakes (Huisman et al., 2004).

Nutrients, temperature, zooplankton and silica were expected to drive the biomass of diatoms (Castenholz, 1961; Egge & Aksnes, 1992; Feminella & Hawkins, 1995; Egge, 1998; Montagnes & Franklin, 2001; Lange et al., 2011). Our model suggested a significant positive relationship of diatom biomass by logZooBiomass, logCalcium and logDepth, while light also determined diatom biomass (Supplementary Table 7). Maximum depth (Zmax) was previously found to drive the abundance of Si-requiring species (Vogt et al., 2017), but contrary to our hypothesis zooplankton had a positive effect. The positive influence of calcium might be due to stochasticity or the correlation with other nutrients, as other nutrients did not make the model.

The biomass of dinoflagellates was expected to be determined by nutrients, temperature, light and calcium (Berman & Rodhe, 1971; Serruya & Pollingher, 1977; Sandgren, 1988; Chapman & Pfester, 1995; Sukenik et al., 2002). Dinoflagellates, who are also able to form blooms, were indeed predicted by light, nutrients and calcium (Supplementary Table 9). These positive relationships between logNTL or logCalcium and dinoflagellate biomass are thus in line with previous research proposing that the vertical migration of the widely researched Peridinium is limited by nutrients (Berman & Rodhe, 1971) and that the whole group of bloom-forming dinoflagellates tend to prefer water with a high calcium concentration (Sandgren, 1988). The hypothesis that temperature would predict dinoflagellate biomass was not confirmed.

Although little is known about the factors that influence euglenophyte biomass, it was expected that nutrients, temperature, light, O2 and CO2 would play a role (Buetow, 1962; Buetow & Schuit, 1968; Kitaya, Azuma, & Kiyota, 2005). The model indeed supported the hypothesis that a higher logNTL or temperature would result in more euglenophyte biomass (Supplementary Table 11), although there was no evidence for an optimum temperature (Kitaya et al., 2005). The biomass of euglenophytes was also largely influenced by light (Secchi and logTurbidity), supporting the suggestion that light intensity steers the growth rate of euglenophytes (Kitaya et al., 2005). The theory that phosporus limits the growth of euglenophytes was not supported by the model (Buetow & Schuit, 1968), as logNTL probably explains a large share of the response variance that logPTL would explain. Lastly, the negative effect of a surplus in oxygen was demonstrated, confirming the laboratory research by Kitaya et al. (2005) that there is an optimum O2 concentration.

As some chrysophytes need silica, it was expected to limit their biomass. Nutrients, temperature, pH and alkalinity were also predicted to play a role in explaining chrysophyte abundance (Roijackers & Kessels, 1986; Johansson & Granéli, 1999), although there is little knowledge on the factors that drive the biomass of chryophytes. Temperature was found to be the most significant driver of chrysophyte biomass (Supplementary Table 13), where higher temperatures limit the abundance of chrysophytes. This implies that further climate warming would threaten the abundance of chrysophytes, but eutrophication does not seem to affect chrysophyte biomass as there was no evidence to support the theory that silica and nutrients would limit the growth of chrysophytes (Johansson & Granéli, 1999). Multi-model interference did suggest that light and zooplankton densities are significant factors. Further research could explore and validate these findings. Beforehand, the influence of light on chlorophyte biomass was expected, together with nutrients, temperature and CO2. Secchi-depth indeed related negatively to chlorophyte biomass (Supplementary Table 15), which is discussed earlier. Besides Secchi-depth, the model suggested that nitrogen was an important predictor, which was in line with previous research who saw that increasing nitrogen concentrations relate to more biomass of two green algae species (Piorreck, Baasch, & Pohl, 1984). Temperature also played a role in green algae abundance, as increasing temperatures related to more biomass. This large-scale observation is similar to experimental studies who noticed a higher growth rate of six freshwater green algae and the marine Ulva sp. under warmer conditions (Duke, Litaker, & Ramus, 1989; Dauta, Devaux, Piquemal, & Boumnich, 1990; Converti, Casazza, Ortiz, Perego, & Del Borghi, 2009). The prediction that CO2 concentrations would drive chlorophyte biomass was not observed, but there was significant evidence for a positive relationship with zooplankton densities.

Secchi, logNTL and temperature were also important in explaining cryptophyte biomass (Supplementary Table 17). There was prior evidence for the limitation of cryptophytes species by light and nitrogen (Knisely & Geller, 1986; Lewitus & Caron, 1990; Gasol et al., 1993; Gervais, 1997; Gervais, 1998), but the negative correlation between temperature and cryptophytes biomass is unprecedented. This might implicate that climate warming has negative effects on cryptophytes abundance. In addition, both oxygen and DOC concentrations negatively

influenced cryptophyte biomass. As mentioned when discussing the model that explains total phytoplankton biomass, high levels of DOC cause shading. This results in less productivity and thus less biomass.

Trait-based groups

Anabaena, Aphanizomenon and Planktothrix were the three most abundant genera found in this dataset (Supplementary Table 28) and these genera contributed most to the abundance of colony-forming phytoplankton. Large colony-forming cyanobacteria like Anabaena, Planktothrix and Aphanizomenon are commonly dominant in freshwater lakes (Dokulil & Teubner, 2000) and this dominance will intensify with warmer climates (Kosten et al., 2012), thus a positive relationship between temperature and the biomass of colony-forming phytoplankton was expected. Also nutrients, temperature, pH and zooplankton grazing were predicted to influence colony-forming phytoplankton (Naselli-Flores & Barone, 2000; Wang, Tang, Wang, & Smith, 2010; Beaulieu et al., 2013). The model indeed showed that the biomass of colony-forming

phytoplankton was predicted by CO2 concentration, nutrients and zooplankton (Supplementary Table 25) and lake morphology (logDepth and logArea) was also included as a driver of colony-forming biomass. It is known that under nutrient-rich conditions, pCO2 limitation resulted in smaller colonies of the large colony-forming cyanobacteria Microcystis spp. (Ma et al., 2014). This is in contrast with our results, where low CO2

concentrations correlated to more colony-forming biomass. Moreover, there was no evidence that changes in temperature would impact the biomass of colony-forming phytoplankton.

The biomass of mixotrophic species was in line with our hypothesis and other research driven by light and pCO2 (Granéli, Carlsson, & Legrand, 1999; Vogt et al., 2017; Hansson, Grossart, del Giorgio, St‐Gelais, & Beisner, 2019), but also lake area and ANC played a significant role (Supplementary Table 19). This research did not find evidence for the hypothesis that temperature and nutrients would affect mixotrophic species (Jansson, Blomqvist, Jonsson, & Bergström, 1996; Granéli et al., 1999; Pålsson & Granéli, 2004; Wilken, Huisman, Naus‐Wiezer, & Van Donk, 2013), although LogNTL was an important factor according to RF analysis but it was not included in the model (Supplementary Table 18 & 19).

Nutrients and temperature were included in the model that explains the biomass of N-fixing phytoplankton species (Supplementary Table 21), which was in line with our hypothesis that nutrients and temperature would influence the biomass of heterocysts (Beaulieu et al., 2013). N-fixing cyanobacterial blooms increase with increasing temperatures and nutrients (Supplementary Table 21), which may lead to adverse ecological effects with climate change. Increasing temperatures and eutrophication already led to an expansion of harmful cyanobacterial blooms and a decrease of phytoplankton diversity in lakes and oceans and this trend will continue (Elliott, Jones, & Thackeray, 2006; Paerl & Huisman, 2009; Kosten et al., 2012). These harmful blooms of cyanobacteria are large, toxic surface-growths that cause oxygen depletion and alter food webs (Paerl & Huisman, 2009). For humans, they pose a major threat to water supplies, fishing, recreation and health (Codd, 1995; Falconer, 1999; Huisman et al., 2005; Paerl & Otten, 2013).

LogDepth, logpCO2, logTurbidity and logSilica also significantly drove heterocyst biomass (Supplementary Table 21). As cyanobacteria have the unique ability to fix nitrogen (Berman-Frank et al., 2003), this model is similar to the one that predicts cyanobacterial biomass (Supplementary Table 5). CO2 (negatively) and silica concentrations (positively) affected the biomass of heterocysts, but no explanation exists in previous literature for these relationships. An important remark is that this study looked at the biomass of species that could fix nitrogen and not at the rate of fixation. Heterocysts will only fix nitrogen under anoxic conditions as it is very resource-expensive (Fay, 1992), but the lakes generally had enough oxygen (Table 1) and there was no relationship between O2 and heterocyst biomass. Future research should therefore look at the factors that drive the biomass of N-fixing species under anoxic conditions.

Each of the analyzed silica-requiring species were diatoms (Supplementary Table 1), so this model only looked at the factors that determine the biomass of the most abundant diatoms. The biomass of silica-requiring species as a whole was expected to be affected by nutrients, silica, CO2 and depth. The model indicated that nutrients, temperature, zooplankton biomass, pCO2 and lake area played a significant role in explaining the biomass of Si-requiring phytoplankton species (Supplementary Table 23). Interestingly, nitrogen induces a decrease of the biomass of Si-requiring species (diatoms), while phosphorus is positively correlated. There is no similar result in

existing literature and a possible ecological explanation is lacking. Surprisingly, silica is not found to drive the abundance of Si-requiring species, but low concentrations of silica may already be sufficient to support diatom growth. Moreover, lake area and zooplankton biomass are strong positive indicators of diatom abundance, which is unique in literature. The decrease of diatom biomass with higher concentrations of CO2 is in agreement with the finding of Vogt et al. (2017). A possible explanation for this relationship is that diatoms lose their advantage of taking up CO2 under low CO2 concentrations (Low-Décarie, Fussmann, & Bell, 2014) and other species therefore can compete with diatoms if the pCO2 is higher.

Moreover, there is a constant decrease of Si-requiring and thus diatom biomass with increasing temperatures (Supplementary Table 7). A reason might be that warmer waters are less viscous than cold water, inducing that large diatoms sink more easily from the euphotic zone. This implicates that further global warming would lead to a decrease in diatom abundance and thus a shift in phytoplankton community.

Genus richness

Genus richness was mostly predicted by water temperature and chlorophyll a (Supplementary Table 27), which was hypothesized based on previous research (Stomp et al., 2011; Vogt et al., 2017). Furthermore, these studies also confirmed our finding that zooplankton biomass positively and phosphorus negatively influenced freshwater phytoplankton genus richness (Muylaert et al., 2010; Vogt et al., 2017). As temperature is linked to latitude, genus richness declines significantly over a latitudinal scale (Figure 5). This latitudinal scale of

freshwater phytoplankton genus richness is known and previously observed in similar research (Stomp et al., 2011).

Conclusion

In conclusion, there are large differences in the factors that drive different phytoplankton groups. Overall, light and nutrients were the main drivers of the abundances of taxonomic groups, while pCO2 and lake morphology were relatively important in explaining trait-based biomasses. Temperature and productivity of a lake largely influenced genus richness. Our results largely supported previous theories about the factors that influence different phytoplankton, but in addition found significant evidence for other possible drivers of phytoplankton communities (Table 6). This study was the first to examine the factors that shape phytoplankton communities on a continental scale and these conclusions add to the understanding of phytoplankton ecology and water management.

Table 6. The factors per phytoplankton group and genus richness that were hypothesized a prior to drive biomass and abundances, based on existing literature. Factors that are marked red were not found to be significant drivers. Green marked factors were factors that were confirmed to play a role in the biomass of that phytoplankton group. Factors that are indicated blue were drivers that are found in this study. * indicate the individual and interactive effect of 2 factors, while : only illustrates the interactive effect.

Group Factors

Total Nutrients Temperature Turbidity*Secchi DOC Depth

Diatoms Nutrients Temperature Zooplankton Silica Light Calcium Depth Cyanobacteria Nutrients Temperature pH Depth Light

Chlorophytes Nutrients Temperature Light CO2 Zoopla_nkton

Dinoflagellates Nutrients Temperature Light Calcium

Euglenophytes Nutrients Temperature Light O2 CO2

Cryptophytes Nutrients Light Zooplankton O2 Temperature DOC

Chrysophytes Nutrients Temperature pH Alkalinity Silica Zooplankton Light

N-fixation Nutrients Temperature Light Silica Depth CO2

Mixotrophy Nutrients Temperature Light CO2 Zooplankton DOM Area ANC

Si-requirement

Nutrients Silica CO2 Depth Zooplankton Temperature Area

Colony formation

Nutrients Temperature pH CO2 Zooplankton Depth Area Genus

richness

Literature

:

Anderson, D. M., Glibert, P. M., & Burkholder, J. M. (2002). Harmful algal blooms and eutrophication: nutrient sources, composition, and consequences. Estuaries, 25(4), 704-726.

Arenovski, A. L. (1994). Distribution, Abundance and Ecology of Mixotrophic Algae in Marine and Freshwater Plankton Communities (No. WHOI-94-22). Woods Hole Oceanographic Institution MA.

Barsanti, L., Vismara, R., Passarelli, V., & Gualtieri, P. (2001). Paramylon (β-1, 3-glucan) content in wild type and WZSL mutant of Euglena gracilis. Effects of growth conditions. Journal of applied phycology, 13(1), 59-65. Barton, K. (2009). MuMIn: multi-model inference. R package version 1. 0. 0. http://r-forge. r-project. org/projects/mumin/.

Beaulieu, M., Pick, F., & Gregory-Eaves, I. (2013). Nutrients and water temperature are significant predictors of cyanobacterial biomass in a 1147 lakes data set. Limnology and Oceanography, 58(5), 1736-1746.

Berman, T., & Rodhe, W. (1971). Distribution and migration of Peridinium in Lake Kinneret: With 7 figures and 1 table in the text. Internationale Vereinigung für Theoretische und Angewandte Limnologie: Mitteilungen, 19(1), 266-276.

Berman-Frank, I., Lundgren, P., & Falkowski, P. (2003). Nitrogen fixation and photosynthetic oxygen evolution in cyanobacteria. Research in microbiology, 154(3), 157-164.

Boynton, W. R., Kemp, W. M., & Keefe, C. W. (1982). A comparative analysis of nutrients and other factors influencing estuarine phytoplankton production. In Estuarine comparisons (pp. 69-90). Academic Press. Buetow, D. E. (1962). Differential effects of temperature on the growth of Euglena gracilis. Experimental cell research, 27(1), 137-142.

Buetow, D. E., & Schuit, K. E. (1968). Phosphorus and the growth of Euglena gracilis. The Journal of protozoology, 15(4), 770-773.

Canfield, D. E., & Hodgson, L. M. (1983). Prediction of Secchi disc depths in Florida lakes: impact of algal biomass and organic color. Hydrobiologia, 99(1), 51-60.

Carey, C. C., Ibelings, B. W., Hoffmann, E. P., Hamilton, D. P., & Brookes, J. D. (2012). Eco-physiological adaptations that favour freshwater cyanobacteria in a changing climate. Water research, 46(5), 1394-1407. Carpenter, S. R., Cole, J. J., Kitchell, J. F., & Pace, M. L. (1998). Impact of dissolved organic carbon, phosphorus, and grazing on phytoplankton biomass and production in experimental lakes. Limnology and Oceanography, 43(1), 73-80.

Castenholz, R. W. (1961). The effect of grazing on marine littoral diatom populations. Ecology, 42(4), 783-794. Chapman, A. D., & Pfester, L. A. (1995). The effects of temperature, irradiance, and nitrogen on the encystment and growth of the freshwater dinoflagellates Peridinium cinctum and P.willei in culture (Dinophyceae) 1. Journal of Phycology, 31(3), 355-359.

Cloern, J. E. (1987). Turbidity as a control on phytoplankton biomass and productivity in estuaries. Continental shelf research, 7(11-12), 1367-1381.

Codd, G. A. (1995). Cyanobacterial toxins: occurrence, properties and biological significance. Water Science and Technology, 32(4), 149-156.

Converti, A., Casazza, A. A., Ortiz, E. Y., Perego, P., & Del Borghi, M. (2009). Effect of temperature and nitrogen concentration on the growth and lipid content of Nannochloropsis oculata and Chlorella vulgaris for biodiesel production. Chemical Engineering and Processing: Process Intensification, 48(6), 1146-1151.

Cremona, F., Tuvikene, L., Haberman, J., Nõges, P., & Nõges, T. (2018). Factors controlling the three-decade long rise in cyanobacteria biomass in a eutrophic shallow lake. Science of the Total Environment, 621, 352-359. Dauta, A., Devaux, J., Piquemal, F., & Boumnich, L. (1990). Growth rate of four freshwater algae in relation to light and temperature. Hydrobiologia, 207(1), 221-226.

Diehl, S., Berger, S., Ptacnik, R., & Wild, A. (2002). Phytoplankton, light, and nutrients in a gradient of mixing depths: field experiments. Ecology, 83(2), 399-411.

Dokulil, M. T., & Teubner, K. (2000). Cyanobacterial dominance in lakes. Hydrobiologia, 438(1-3), 1-12.

Duke, C. S., Litaker, W., & Ramus, J. (1989). Effect of temperature on nitrogen-limited growth rate and chemical composition of Ulva curvata (Ulvales: Chlorophyta). Marine Biology, 100(2), 143-150.

Edmondson, W. T. (1972). Nutrients and phytoplankton in Lake Washington.

Egge, J. K. (1998). Are diatoms poor competitors at low phosphate concentrations?. Journal of Marine Systems, 16(3-4), 191-198.

Egge, J. K., & Aksnes, D. L. (1992). Silicate as regulating nutrient in phytoplankton competition. Marine ecology progress series. Oldendorf, 83(2), 281-289.

Elliott, J. A., Jones, I. D., & Thackeray, S. J. (2006). Testing the sensitivity of phytoplankton communities to changes in water temperature and nutrient load, in a temperate lake. Hydrobiologia, 559(1), 401-411. Falconer, I. R. (1999). An overview of problems caused by toxic blue–green algae (cyanobacteria) in drinking and recreational water. Environmental Toxicology: An International Journal, 14(1), 5-12.

Falkowski, P. G., Katz, M. E., Knoll, A. H., Quigg, A., Raven, J. A., Schofield, O., & Taylor, F. J. R. (2004). The evolution of modern eukaryotic phytoplankton. science, 305(5682), 354-360.

Fay, P. (1992). Oxygen relations of nitrogen fixation in cyanobacteria. Microbiology and Molecular Biology Reviews, 56(2), 340-373.

Feld, C. K., Segurado, P., & Gutiérrez-Cánovas, C. (2016). Analysing the impact of multiple stressors in aquatic biomonitoring data: A ‘cookbook’with applications in R. Science of the Total Environment, 573, 1320-1339. Feminella, J. W., & Hawkins, C. P. (1995). Interactions between stream herbivores and periphyton: a

quantitative analysis of past experiments. Journal of the North American Benthological Society, 14(4), 465-509. Fichez, R., Jickells, T. D., & Edmunds, H. M. (1992). Algal blooms in high turbidity, a result of the conflicting consequences of turbulence on nutrient cycling in a shallow water estuary. Estuarine, Coastal and Shelf Science, 35(6), 577-592.

Gasol, J. M., García-Cantizano, J., Massana, R., Guerrero, R., & Pedrós-Alió, C. (1993). Physiological ecology of a metalimnetic Cryptomonas population: relationships to light, sulfide and nutrients. Journal of Plankton Research, 15(3), 255-275.

Gervais, F. (1997). Light‐dependent growth, dark survival, and glucose uptake by cryptophytes isolated from a freshwater chemocline. Journal of Phycology, 33(1), 18-25.

Gervais, F. (1998). Ecology of cryptophytes coexisting near a freshwater chemocline. Freshwater Biology, 39(1), 61-78.

Granéli, E., Carlsson, P., & Legrand, C. (1999). The role of C, N and P in dissolved and particulate organic matter as a nutrient source for phytoplankton growth, including toxic species. Aquatic Ecology, 33(1), 17-27.

Guiry, M. D. (2012). How many species of algae are there?. Journal of phycology, 48(5), 1057-1063.

Guiry, M.D. & Guiry, G.M. (2020). AlgaeBase. World-wide electronic publication, National University of Ireland, Galway. https://www.algaebase.org; searched on 05 April 2020.

Hansson, T. H., Grossart, H. P., del Giorgio, P. A., St‐Gelais, N. F., & Beisner, B. E. (2019). Environmental drivers of mixotrophs in boreal lakes. Limnology and Oceanography, 64(4), 1688-1705.

Hibberd, D. J. (1976). The ultrastructure and taxonomy of the Chrysophyceae and Prymnesiophyceae

(Haptophyceae): a survey with some new observations on the ultrastructure of the Chrysophyceae. Botanical Journal of the Linnean Society, 72(2), 55–80. doi:10.1111/j.1095-8339.1976.tb01352.x

Hoek, C., Mann, D., Jahns, H. M., & Jahns, M. (1995). Algae: an introduction to phycology. Cambridge university press.

Huisman, J.M., Matthijs, H.C.P., and Visser, P.M. (2005) Harmful Cyanobacteria. Springer Aquatic Ecology Series 3. Dordrecht, the Netherlands: Springer.

Huisman, J., Sharples, J., Stroom, J. M., Visser, P. M., Kardinaal, W. E. A., Verspagen, J. M., & Sommeijer, B. (2004). Changes in turbulent mixing shift competition for light between phytoplankton species. Ecology, 85(11), 2960-2970.

Ishwaran, H. & Kogalur, U. B. (2016). RandomForestSRC: Random Forest for Survival, Regression and Classification (RF-SRC). R package version 2.2. 0.

Jäger, C. G., Diehl, S., & Schmidt, G. M. (2008). Influence of water‐column depth and mixing on phytoplankton biomass, community composition, and nutrients. Limnology and oceanography, 53(6), 2361-2373.

Jansson, M., Blomqvist, P., Jonsson, A., & Bergström, A. K. (1996). Nutrient limitation of bacterioplankton, autotrophic and mixotrophic phytoplankton, and heterotrophic nanoflagellates in Lake Örträsket. Limnology and Oceanography, 41(7), 1552-1559.

Järvenpää, M., & Lindström, K. (2004). Water turbidity by algal blooms causes mating system breakdown in a shallow-water fish, the sand goby Pomatoschistus minutus. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1555), 2361-2365.

Jensen, J. P., Jeppesen, E., Olrik, K., & Kristensen, P. (1994). Impact of nutrients and physical factors on the shift from cyanobacterial to chlorophyte dominance in shallow Danish lakes. Canadian Journal of fisheries and aquatic sciences, 51(8), 1692-1699.

Johansson, N., & Granéli, E. (1999). Cell density, chemical composition and toxicity of Chrysochromulina polylepis (Haptophyta) in relation to different N: P supply ratios. Marine Biology, 135(2), 209-217. John, D. M., Whitton, B. A., & Brook, A. J. (2011). Freshwater algal flora of the British Isles. Cambridge University Press.

Kitaya, Y., Azuma, H., & Kiyota, M. (2005). Effects of temperature, CO2/O2 concentrations and light intensity on cellular multiplication of microalgae, Euglena gracilis. Advances in Space Research, 35(9), 1584-1588.

Knisely, K., & Geller, W. (1986). Selective feeding of four zooplankton species on natural lake phytoplankton. Oecologia, 69(1), 86-94.

Kosten, S., Huszar, V. L., Bécares, E., Costa, L. S., van Donk, E., Hansson, L. A., ... & De Meester, L. (2012). Warmer climates boost cyanobacterial dominance in shallow lakes. Global Change Biology, 18(1), 118-126. Kristiansen, J. R. (1986). Silica-scale bearing chrysophytes as environmental indicators. British phycological journal, 21(4), 425-436.

Lange, K., Liess, A., Piggott, J. J., Townsend, C. R., & Matthaei, C. D. (2011). Light, nutrients and grazing interact to determine stream diatom community composition and functional group structure. Freshwater Biology, 56(2), 264-278.

Lee, Z., Shang, S., Du, K., & Wei, J. (2018). Resolving the long‐standing puzzles about the observed Secchi depth relationships. Limnology and Oceanography, 63(6), 2321-2336.

Leliaert, F., Smith, D. R., Moreau, H., Herron, M. D., Verbruggen, H., Delwiche, C. F., & De Clerck, O. (2012). Phylogeny and molecular evolution of the green algae. Critical reviews in plant sciences, 31(1), 1-46.

Lewitus, A. J., & Caron, D. A. (1990). Relative effects of nitrogen or phosphorus depletion and light intensity on the pigmentation, chemical composition, and volume of Pyrenomonas salina (Cryptophyceae). Marine Ecology Progress Series, 171-181.

Litchman, E., & Klausmeier, C. A. (2008). Trait-based community ecology of phytoplankton. Annual review of ecology, evolution, and systematics, 39, 615-639.

Low-Décarie, E., Fussmann, G. F., & Bell, G. (2014). Aquatic primary production in a high-CO2 world. Trends in Ecology & Evolution, 29(4), 223-232.

Ma, J., Brookes, J. D., Qin, B., Paerl, H. W., Gao, G., Wu, P., ... & Xu, H. (2014). Environmental factors controlling colony formation in blooms of the cyanobacteria Microcystis spp. in Lake Taihu, China. Harmful algae, 31, 136-142.

Mann, D. G. (1999). The species concept in diatoms. Phycologia, 38(6), 437-495.

Moisander, P. H., Steppe, T. F., Hall, N. S., Kuparinen, J., & Paerl, H. W. (2003). Variability in nitrogen and phosphorus limitation for Baltic Sea phytoplankton during nitrogen-fixing cyanobacterial blooms. Marine Ecology Progress Series, 262, 81-95.

Montagnes, D. J., & Franklin, M. (2001). Effect of temperature on diatom volume, growth rate, and carbon and nitrogen content: reconsidering some paradigms. Limnology and Oceanography, 46(8), 2008-2018.

Muylaert, K., Pérez-Martínez, C., Sánchez-Castillo, P., Lauridsen, T. L., Vanderstukken, M., Declerck, S. A., ... & Vyverman, W. (2010). Influence of nutrients, submerged macrophytes and zooplankton grazing on

phytoplankton biomass and diversity along a latitudinal gradient in Europe. In Fifty years after the ‘‘Homage to Santa Rosalia’’: Old and new paradigms on biodiversity in aquatic ecosystems (pp. 79-90). Springer, Dordrecht. Naimi, B. (2015). usdm: uncertainty analysis for species distribution models. R package version 1.1-15. R Documentation http://www.rdocu‑mentation. org/packages/usdm.

Naselli-Flores, L., & Barone, R. (2000). Phytoplankton dynamics and structure: a comparative analysis in natural and man-made water bodies of different trophic state. Hydrobiologia, 438(1-3), 65-74.

Nelson, D. M., Tréguer, P., Brzezinski, M. A., Leynaert, A., & Quéguiner, B. (1995). Production and dissolution of biogenic silica in the ocean: revised global estimates, comparison with regional data and relationship to biogenic sedimentation. Global Biogeochemical Cycles, 9(3), 359-372.

Olaveson, M. M., & Nalewajko, C. (2000). Effects of acidity on the growth of two Euglena species. Hydrobiologia, 433(1-3), 39-56.

O’neil, J. M., Davis, T. W., Burford, M. A., & Gobler, C. J. (2012). The rise of harmful cyanobacteria blooms: the potential roles of eutrophication and climate change. Harmful algae, 14, 313-334.

Paerl, H. W., & Huisman, J. (2009). Climate change: a catalyst for global expansion of harmful cyanobacterial blooms. Environmental microbiology reports, 1(1), 27-37.

Paerl, H. W., & Otten, T. G. (2013). Harmful cyanobacterial blooms: causes, consequences, and controls. Microbial ecology, 65(4), 995-1010.

Pålsson, C., & Granéli, W. (2004). Nutrient limitation of autotrophic and mixotrophic phytoplankton in a temperate and tropical humic lake gradient. Journal of Plankton Research, 26(9), 1005-1014.

Peterson, B. G., Carl, P., Boudt, K., Bennett, R., Ulrich, J., Zivot, E., ... & Wuertz, D. (2018). Package ‘PerformanceAnalytics’. R Team Cooperation.

Piorreck, M., Baasch, K. H., & Pohl, P. (1984). Biomass production, total protein, chlorophylls, lipids and fatty acids of freshwater green and blue-green algae under different nitrogen regimes. Phytochemistry, 23(2), 207-216.

Ptacnik, R., Solimini, A. G., Andersen, T., Tamminen, T., Brettum, P., Lepistö, L., ... & Rekolainen, S. (2008). Diversity predicts stability and resource use efficiency in natural phytoplankton communities. Proceedings of the National Academy of Sciences, 105(13), 5134-5138.

Reynolds, C. S. (1984). The ecology of freshwater phytoplankton. Cambridge University Press.

Rigosi, A., Carey, C. C., Ibelings, B. W., & Brookes, J. D. (2014). The interaction between climate warming and eutrophication to promote cyanobacteria is dependent on trophic state and varies among taxa. Limnology and Oceanography, 59(1), 99-114.

Roijackers, R. M. M., & Kessels, H. (1986). Ecological characteristics of scale‐bearing Chrysophyceae from The Netherlands. Nordic journal of botany, 6(3), 373-385.

Round, F. E., Crawford, R. M., & Mann, D. G. (2007). Diatoms: biology and morphology of the genera. Cambridge university press.

Sandgren, C. D. (Ed.). (1988). Growth and reproductive strategies of freshwater phytoplankton. Cambridge University Press.

Scheffer, M., Rinaldi, S., Gragnani, A., Mur, L. R., & van Nes, E. H. (1997). On the dominance of filamentous cyanobacteria in shallow, turbid lakes. Ecology, 78(1), 272-282.

Schelske, C. L., & Stoermer, E. F. (1971). Eutrophication, silica depletion, and predicted changes in algal quality in Lake Michigan. Science, 173(3995), 423-424.

Schindler, D. W., Hecky, R. E., Findlay, D. L., Stainton, M. P., Parker, B. R., Paterson, M. J., ... & Kasian, S. E. M. (2008). Eutrophication of lakes cannot be controlled by reducing nitrogen input: results of a 37-year whole-ecosystem experiment. Proceedings of the National Academy of Sciences, 105(32), 11254-11258

Schöllhorn, E., & Granéli, E. (1996). Influence of different nitrogen to silica ratios and artificial mixing on the structure of a summer phytoplankton community from the Swedish west coast (Gullmar Fjord). Journal of Sea Research, 35(1-3), 159-167.

Serruya, C., & Pollingher, U. (1977). Lowering of water level and algal biomass in Lake Kinneret. Hydrobiologia, 54(1), 73-80.

Shilo, M. (1971). Toxins of chrysophyceae. Microbial toxins, 3, 67-103.

Smetacek, V. (1999). Diatoms and the ocean carbon cycle. Protist, 150(1), 25-32.

Smith, V. H. (2003). Eutrophication of freshwater and coastal marine ecosystems a global problem. Environmental Science and Pollution Research, 10(2), 126-139.

Stomp, M., Huisman, J., Mittelbach, G. G., Litchman, E., & Klausmeier, C. A. (2011). Large‐scale biodiversity patterns in freshwater phytoplankton. Ecology, 92(11), 2096-2107.

Sukenik, A., Eshkol, R., Livne, A., Hadas, O., Rom, M., Tchernov, D., ... & Kaplan, A. (2002). Inhibition of growth and photosynthesis of the dinoflagellate Peridinium gatunense by Microcystis sp.(cyanobacteria): a novel allelopathic mechanism. Limnology and Oceanography, 47(6), 1656-1663.

Team, R. C. (2013). R: A language and environment for statistical computing.

Tomas, C. R. (1978). Olisthodiscus luteus (Chrysophyceae) I. Effects of salinity and temperature on growth. Motility and survivals 1, 2. Journal of Phycology, 14(3), 309-313.