Lakes can be regarded as islands within a matrix of terrestrial land cover (Dodson, 1992) and numerous observational lake gradient studies support species area theory: aquatic bacteria (Bell et al., 2005; Reche et al., 2005), phytoplankton (Søndergaard et al., 2005), zooplankton (Dodson, 1992; Keller and Conlon, 1994), fish (Tonn and Magnuson, 1982; Post et al., 2000; Ingram et al., 2012), submerged vegetation (Rørslett, 1991) and birds (Hoyer and Canfield, 1994). Despite the consistency of results across taxonomic groups, species-area relationship has traditional been tested in a two-dimensional habitats, which do not hold in lenthic ecosystems as they have an additional vertical dimension (water column), making diversity predictions more complex (Drakare et al., 2006). The species sorting paradigm of metacommunity ecology stresses the dominant role of local habitat factors on community composition and biodiversity (Leibold et al., 2004; Elliott et al., 2007; Soininen, 2014). For instance, in lakes an assessment of size effects on species richness can be confounded by the interactions of lake volume with other variables (acidity and nutrient status) (Sprules, 1977; Confer et al., 1983; Søndergaard et al., 2005). The importance of many of these in-lake variables (e.g., trophic condition) can be magnified through global change stressors (e.g., eutrophication) that can result in dramatic changes in the abiotic and biotic environment of lakes (Hansen et al., 2001). The resulting changes may lead to interactions with ecosystem size that could affect community structure (Angeler and Alvarez-Cobelas, 2005). Cause-effect relationships between biodiversity, natural and anthropogenic factors cannot be explored by observational studies alone (Frost et al., 1988). Experimental approaches can discern whether the role of global stress impacts is dependent on ecosystem size by controlling for potential confounding effects of spatial context and local variability in the abiotic and biotic environment.
In this experimental study, we manipulated water volume in mesocosms, to assess how volume (for simplicity referred to as size), affect phytoplankton communities that dwell in three dimensional habitats (Marquet et al., 1993; Carr et al., 2003). We also manipulated nutrient levels to mimic mesotrophic and hypertrophic conditions in lakes, to investigate how the interaction between ecosystem size and anthropogenic stressors can influence phytoplankton diversity. Shallow lakes are the predominant lake type worldwide (Wetzel, 2001), and eutrophication impacts are predicted to increase due to changes in temperature and precipitation regimes associated with climate change, which is influenced by the transport of nutrients from the land to lakes and prolonged water residence time, changing external and internal nutrient loading patterns (Meyer et al., 1999; Ipcc, 2007; Giorgi and Lionello, 2008; Wilhelm and Adrian, 2008; Jeppesen et al., 2009). Thus, the importance of understanding how the interactions between cultural eutrophication and lake size can influence diversity is relevant for society in terms of ecosystem service provisioning (Scheffer et al., 2001). Cultural eutrophication is well-known to induce regime shifts from clear to turbid water in lakes, resulting in a drastic change in ecosystem services such as fisheries, recreation and fresh water quality (Pope et al., 2014).
Our goal was to investigate how the interaction between a stressor (eutrophication) and lake size can influence phytoplankton community composition and diversity using common metrics of species richness, diversity, evenness and total biovolume. Phytoplankton communities respond quickly to environmental change and are indicators of eutrophication (Scheffer et al., 1993; Carpenter et al., 1999; Carvalho et al., 2013). They also show different community dynamics in contrasting ecosystem states (Angeler et al., 2013; Baho et al., 2014), where high nutrient levels generally favour certain species, for instance cyanobacteria. However, the role of ecosystem size in community assembly under the influence of anthropogenic stressors is not comprehensively explored. Here we test the hypothesis that the influence of eutrophication on phytoplankton is dependent on lake size and expect differences in the dynamics of key species and community composition.
The mesocosm experiment was conducted in Lake Erken, Sweden for approximately six months (May to early November 2011) corresponding to the ice-free period. Lake Erken is a mesotrophic lake (annual average of 25 µg TP L–1) located southeast of Sweden (59°49’59” N, 18°33’55” E) and has a surface area of 24 km2, an average depth of 9 m and a maximum depth of 21 m (Elliott et al., 2007). Ice cover normally occurs between December and April and the lake has an average water retention time of seven years (Elliott et al., 2007).
A thorough description of the experimental design can be found in Landkildehus et al. (2014). Briefly, the experiment consisted of two nutrient levels (low, high) crossed with two water depths (shallow, deep). Each treatment [shallow - high nutrients (SH); shallow - low nutrients (SL); deep - high nutrients (DH); deep - low nutrients (DL)] was replicated four times, resulting in a total of 16 mesocosms. The mesocosms consisted of cylindrical impermeable fiberglass tanks (diameter 1.2 m and a thickness of 4 mm) with two different depths 1.2 and 2.2 m corresponding to the shallow and deep treatments, respectively. A combination of sand and mud from a nearby mesotrophic lake was added as bottom substratum (ca. 10 cm) before being filled with filtered (500 µm) lake water to 1020 and 2150 litres, corresponding to two contrasting ecosystem sizes (shallow and deep water levels, respectively). The mesocosms received inocula of plankton (phytoplankton and zooplankton) from five nearby lakes, to enable the development of a diverse flora and fauna, thus addressing one of the major shortcomings of mesocosm experiments regarding the limited representation of species (Stewart et al., 2013). Plankton inocula were collected from five nearby lakes by performing five vertical hauls over the entire water column, using plankton nets of mesh size of 50 µm. The contents were pooled, mixed and a 1-litre subsample was added to each mesocosm. Similarly, five litres of sediments were collected from the same five lakes and carefully sieved (through 10 mm) to avoid introducing fish eggs or large mussels. Submerged vegetation, Eurasian water-milfoil (Myriophyllum spicatum), was collected from the host lake and planted in each mesocosm. Planktivorous fish (juvenile roach (Rutilus rutilus) originating from Lake Erken) were stocked in each mesocosm based on fish biomass (ca. 10 grams per enclosure) and dead fish were replaced when necessary (see Supplementary Table S1-S3 for information about the length and weight of fish individuals). Nets were installed on top of each enclosure to deter birds.
Phosphorus concentrations were adjusted to 25 µg TP L–1 and 200 µg TP L–1, for the low and high nutrients mimicking mesotophic and hypertrophic conditions, respectively. Nitrogen was added to achieve phosphorus to nitrogen ratio of 1:20 (by molecular weight). At the beginning of the experiment, nutrient additions were only required for the high nutrient treatments. Thereafter, monthly nutrient additions were needed to maintain the two nutrient levels throughout the duration of the experiment. Nutrient loss was compensated by monthly nutrient (total phosphorus and total nitrogen) dosing (see Supplementary Table S4 and Figure S1 for observed nutrients values). Nutrient compensation was determined from previous experiments (González Sagrario et al., 2005; Jeppesen et al., 2007b). The water in each mesocosm was continuously mixed using standard 5 watt aquaria pumps.
Sampling of phytoplankton was initiated on 16th of May 2011. Thereafter samples were taken twice per month resulting in 13 sampling points. Water samples representative of the entire water column were taken with the aid of a Plexiglas® tube (length 60 cm and internal diameter of 9.5 cm) sampler and pooled. The pooled samples were homogenized and transferred to 50 ml glass bottles and 500 ml plastic bottles for phytoplankton and water chemistry analyses, respectively. The phytoplankton samples were preserved with 0.5 ml of acid Lugol’s solutions (Olrik et al., 1989). Phytoplankton was counted using an inverted light microscope and a modified Utermöhl technique commonly used in Scandinavia (Olrik et al., 1989). Taxa were identified to the lowest possible taxonomic unit (generally species). Biovolumes (mm3 L–1) were calculated from geometric conversions following protocols developed by Blomqvist and Herlitz (1998). All water chemistry analyses were conducted by SWEDAC certified laboratories (Swedish Board for Accreditation and Conformity Assessment, SWEDAC; http://swedac.se/en/) at the Erken laboratory, using International (ISO) and European standards when available (http://search.swedac.se/en/accreditations/1239/a000043-001).
All statistical analyses were conducted in R 3.0.2 (R Development Core Team, 2012). Some mesocosms (one shallow - low nutrient and one deep - low nutrient) were excluded from the statistical analyses due to birds that used the pontoon as resting ground disturbing the experimental nutrient balance with their excrements. In one tank (deep and low nutrient), a non-stocked fish (an adult perch 150 g) was discovered and therefore this tank was excluded from analyses. The influence of the birds (despite using nets as a deterrent) and the non-stocked fish could be detected by means of non-metric multidimensional scaling ordination of water quality variables. As a solution, we randomly removed one replicate from each treatment to achieve a balanced experimental design.
Phytoplankton community structure across treatments was characterized with commonly used metrics (richness, diversity, evenness and total biovolume) following the recommendations of Jost (2007) and Tuomisto (2010) to obtain statistically unbiased measures. The Shannon index was transformed into species equivalents by exponentiation to make species richness and diversity trends comparable (Jost, 2007). The exponentiated Shannon index was considered as a measure of ‘diversity’ reflecting a different ecological process from species richness (‘richness’) (Tuomisto, 2010). Evenness was obtained dividing ‘diversity’ by ‘richness’, thereby expressing evenness independently from richness (Tuomisto, 2012). We also assessed temporal community turnover across treatments using Bray-Curtis and Sørensen dissimilarities.
Repeated measures analysis of variance (rm-ANOVA) was used to test for statistical differences in the phytoplankton community metrics between treatments. The treatment effects, time effect and their interaction (treatment x time) were investigated. Log transformation was carried out when necessary to fulfil the assumptions of parametric tests. The assumptions of sphericity were similarly considered and when violated Huynh Feldt corrections were implemented (Huynh and Feldt, 1976). When a significant main effect was detected, a posteriori comparisons were made using Tukey HSD test.
The rm-ANOVA was complemented with permutational multivariate analysis of variance (PERMANOVA) followed by non-metric multidimensional scaling (NMDS) ordination to contrast phytoplankton communities across treatments and time. PERMANOVAs were based on Bray-Curtis and Sørensen similarity matrices of square-root transformed species biovolume matrices with 9999 unrestricted permutations of raw data. Significant differences were inferred at an α-level of 0.05. For both univariate and multivariate comparisons, the interaction terms were considered crucial for inference because they indicate that phytoplankton structure and their temporal patterns during the experiment in low and high nutrient mesocosms vary depending on the size of the mesocosms.
The trajectory of the phytoplankton communities’ development across treatments (averaged across replicates) and time were visually assessed using a NMDS ordination based on Bray-Curtis distance matrix. This ordination technique starts by scaling objects in full-dimensional space then reduces these dimensions to a few while maintaining the distance relationship among the objects (Legendre and Legendre, 1998). Spearman’s rank correlation analysis was used to investigate the relationship between the individual phytoplankton species (raw biovolume data of individual species) and the univariate community metrics (species richness, diversity, evenness and total biovolume) with NMDS dimensions 1 and 2.
Univariate analyses showed significant treatment effects for all metrics (Tab. 1, Fig. 1): species richness, diversity, evenness and total biovolume. Post-hoc analysis revealed that the species richness, diversity and evenness were higher in the low nutrient treatments compared to the high nutrient treatments irrespective of water depth (Tukey’s HSD test: DL>DH, SL>DH, SL>SH, DL>SH, SH=DH, SL=DL; P<0.05), while the reverse trend was found for total biovolume (Tukey’s HSD test: DL<DH, SL<DH, SL<SH, DL<SH, SH=DH, SL=DL; P<0.05) (Fig. 2). Moreover, the effect of time was significant for all community metrics, whereas the interaction term (treatment x time) was significant for species richness and diversity.
Cyanobacteria were the predominant group in the deep and high nutrient mesocosms (Fig. 3), while Charophyta were generally more abundant in the shallow and high nutrient mesocosms. Cyanobacteria peaked almost simultaneously in both high nutrient treatments and reached higher biovolumes in the deeper mesocosms. In contrast, the two low nutrient treatments had comparable phytoplankton development with the exception of Chlorophyta that had higher biovolume in the shallow compared to the deep mesocosms.
Permutational multivariate analysis of variance based on both Bray-Curtis (considering both presence/absence and biovolume of taxa) and Sørensen dissimilarity matrices (emphasizing presence/absence only) detected significant treatment and time effects, with a significant interaction (Tab. 2). Phytoplankton communities differed among the four treatments (PERMANOVA analysis, P<0.05). Effects of nutrient and depth manipulations were visually apparent in non-metric multidimensional scaling (NMDS) ordination when tracking phytoplankton communities’ development over time (Fig. 4), which corroborate the results of the PERMAVONA analyses. Phytoplankton communities, similar at the start of the experiment, showed different degrees of divergence over time depending on the treatment (Fig. 4). Nutrient effects were visible along the second NMDS axis, where the low nutrient treatments tended to aggregate towards the top of the figure, whilst the high nutrient treatment clustered in lower right corner (Fig. 4). The Spearman rank correlation analyses showed that 10 species correlated significantly with the first NMDS axis whereas 20 taxa and total biovolume with the second NMDS axis (Tab. 3).
Generally, community analyses indicated declining turnover rates through time. By contrast, high nutrient treatments exhibited slower turnover rates compared to the low nutrient treatments. The two dissimilarity matrices indicated that treatment effects were only significant for Sørensen (Tukey’s HSD test: DL>DH, SL>DH; P<0.05), whereas the time effect and the interaction term were significant for both matrices (Fig. 5, Tab. 4).
Nutrient enrichment decreased species richness, diversity and evenness but increased the biovolume of phytoplankton, which corroborates with previous studies (Interlandi and Kilham, 2001; Smith, 2003; Smith and Schindler, 2009). Our findings suggest that common metrics of phytoplankton biodiversity (e.g., species richness and diversity) were similar between size treatments with the same nutrient status. This deviates from the results of previous lake studies that supported the species area theory (Tonn and Magnuson, 1982; Rørslett, 1991; Dodson, 1992; Hoyer and Canfield, 1994; Keller and Conlon, 1994; Post et al., 2000; Bell et al., 2005; Reche et al., 2005; Søndergaard et al., 2005). Despite these findings, phytoplankton species composition markedly differed across size treatments. This does not only support previous studies that inference varies with metrics chosen (univariate vs. multivariate) (Angeler and Goedkoop, 2010), but also with the volume of the water column in mesocosms. These findings complement and refine insight from previous lakes studies that focused, like terrestrial studies, on the two-dimensional aspect of surface area (Marquet et al., 1993; Carr et al., 2003).
Our results suggest that further cross-system comparisons in the future (e.g., Drakare, Lennon & Hillebrand, 2006) can benefit from incorporating ecosystems with three-dimensional habitat configuration, specifically marine and freshwater environments but also terrestrial habitats where such dimensionality can vary substantially (e.g., height of grasslands vs forests). At least in a planktonic environment, the three-dimensional habitat configuration mediates important abiotic influences on phytoplankton; that is, light, pH, temperature, oxygen, nutrients and stratification, can show important gradients in the water column (Kalff, 2002; Angeler and Rodrigo, 2004; Angeler et al., 2004; Ingram et al., 2012), which enhances habitat heterogeneity and presumably beta-diversity (Veech and Crist, 2007). In this study, we had identical mixing regimes induced by pumps, and therefore can discard stratification as a variable influencing community dynamics. Instead, light regime in the water column comparable to natural conditions might have been important for phytoplankton distribution and community composition in the water column (Havens et al., 1998; Diehl, 2002).
The differences in phytoplankton community composition between treatments (both univariate and multivariate analyses), indicate interactions between habitat size (volume) and nutrient enrichment. Nevertheless, phytoplankton community development over time appeared to follow the nutrient model proposed by Reynolds (1980, 1984) starting with species generally associated with spring bloom, followed by green algae and cyanobacteria. Generally, the effects of nutrient enrichment differed with depth of the water column. Cyanobacteria (mainly Microcystis aeruginosa) were predominant in deep mesocosms, whereas in the shallow nutrient-enriched mesocosms species belonging to the phylum Charophyta (strictly planktonic algae such as species from the genera Euastrum, Elakatothrix and Staurodesmus) reached higher densities than cyanobacteria, presumably due to their preference for higher light intensities (Kruk and Segura, 2012). The inability of cyanobacteria to prevail in the shallow eutrophic mesocosms relative to the deep ones might be attributed to photoinhibition mechanisms, where high light intensities reduce their ability to photosynthesize (Scheffer et al., 1997; Havens et al., 1998). However, in favourable conditions (high nutrient and low light intensity), cyanobacteria can to some extent control light regime through self-shading mechanisms by reaching high densities and forming large colonies, where self-shading does not impact their photosynthetic efficiencies as shown by Agusti and Phlips (1992). This may underline a positive feedback that suppresses photoinhibition and create optimum conditions for cyanobacteria to thrive, which has been previously been demonstrated by Veraart et al. (2012).
Anthropogenic disturbances can potentially homogenize communities by favouring competitively superior species (Devictor et al., 2008; Macdougall et al., 2013). Our findings that nutrient enriched mesocosms had lower evenness and community turnover compared to the low nutrient treatments support this conjecture. Despite this, phytoplankton community structure differed markedly between shallow and deep mesocosms with high nutrient status. However, differences were less pronounced between low-nutrient mesocosms of different size, indicating a strong effect of ecosystem size in mediating eutrophication effects in shallow lakes. The decrease in community turnover could be related to seasonal changes. Although a longer study period might provide more insights to the effect of anthropogenic stressors on temporal community turnover, our results partially suggest that nutrient enrichment can interfere and interrupt the gradual seasonal phytoplankton community dynamics by enabling certain species to maintain prolonged dominance (Scheffer and Nes, 2007).
The scale of perception of organisms, which is correlated to their body size (Peters, 1986) might have some implications for species-area relationship concept (Azovsky, 2002). For instance, the study performed by Azovsky (2002) found that smaller organisms (ciliates and diatoms) showed a weaker and less steep increase of species richness with area than larger ones. Thus, for large species groups such as fish, area might be more appropriate whereas volume, better at expressing habitat heterogeneity across the water column, might be more relevant for smaller planktonic organisms. Developing a framework to evaluate habitat heterogeneity in lentic ecosystems, despite not being the scope this study, could contribute to the ongoing discussion on whether heterogeneity is the most important factor that shape species-area relationship (Macarthur and Wilson, 1967; Tews et al., 2004; Whittaker and Fernández-Palacios, 2007; Báldi, 2008).
In summary, we found that the effects of nutrient enrichment on phytoplankton diversity and communities were dependent on ecosystem size. The dominance of cyanobacteria, often considered a nuisance algal group by water quality managers (Scheffer et al., 1997) due to their wide-ranging impacts, including toxic effects, shading of submerged plants and alteration of food webs (Havens, 2008), is not only dependent on nutrient status but rather on the combination of nutrient and ecosystem size as observed in our experiment. Such interactions increase the uncertainty of global change outcomes related to continued changes in nutrient regimes of lakes (Anderson et al., 2005; Jeppesen et al., 2007a; Jeppesen et al., 2007b) and should therefore be given greater consideration in risk/impact assessment and management. Furthermore, identifying interactions between stressors and ecosystem size can help managers and decision maker to decide when and where to act, as resources (financial and time) are limited (Wilson et al., 2006; Brown et al., 2013; Ghedini et al., 2013; Côté et al., 2016). Priority action must be identified aiming at achieving the highest benefits to with respect to biodiversity and ecosystem functions that are highly valued by society (Wilson et al., 2006; Brown et al., 2013; Ghedini et al., 2013; Côté et al., 2016).