Phytoplankton and littoral invertebrate assemblages in 4 boreal lakes recovering from acidification and 4 minimally disturbed reference lakes studied over 2 decades were used to determine the pathways and trajectories of change under the influence of climatic variability. Assemblage composition (species presence–absence data) but not dominance patterns (invertebrate abundance/phytoplankton biovolume) of acidified lakes became more similar to those of reference lakes (distance decreased with time), indicating that detection of recovery varies as a function of chosen metrics. Acidified lakes had more pronounced shifts in assemblage composition than did reference lakes. The most marked differences were noted for phytoplankton assemblages. Assemblages in acidified lakes had mean between-year Euclidean distances almost 2× greater than those of assemblages in reference lakes. Trends in water chemistry showed unequivocal recovery, but responses of phytoplankton and invertebrate assemblages, measured as between-year shifts in assemblage composition, were correlated with interannual variability in climate (e.g., North Atlantic Oscillation, water temperature) in addition to decreased acidity. The finding that recovery pathways and trajectories of individual acidified lakes and the environmental drivers explaining these changes differed among assemblages shows that biological recovery is complex and the influence of climatic variability is poorly understood.
Acidification of inland waters has been recognized as a major environmental problem for >3 decades. Awareness that emissions from burning of fossil fuels were causing biodiversity loss in lakes and streams led to international action plans to protect and restore natural resources (Stoddard et al. 1999). Ensuing reductions in the emissions of sulfur dioxide and NOx compounds resulted in improved air quality, and subsequent increases in surface-water pH and alkalinity have been attributed to the decreased deposition of these acidifying compounds (e.g., Stoddard et al. 1999, Skjelkvåle et al. 2005). However, whereas a growing body of literature supports chemical recovery of surface waters, empirical evidence of biological recovery is scarce and findings are equivocal (Skjelkvåle et al. 2003, Stendera and Johnson 2008, Ormerod and Durance 2009). Furthermore, because biological recovery is often the ultimate goal of legislative action, failure to achieve biological objectives means that acidification is still considered a foremost problem affecting the biodiversity of inland surface waters in northern Europe (Johnson et al. 2003) and elsewhere (e.g., Monteith et al. 2005, Kowalik et al. 2007, Burns et al. 2008).
Acidification, a prolonged press type of disturbance (sensu Lake 2000), often results in distinctive and persistent changes in biological assemblages (Økland and Økland 1986). Response to disturbance often is stabilized by ecological feedback mechanisms resulting in resilience (e.g., Blindow et al. 2002), and growing evidence indicates that perturbed systems resist returning to the predisturbed state. For example, acidified sites can be maintained in the impaired state by inadequate physicochemical properties (e.g., acid episodes can retard recovery, Kowalik et al. 2007) and biological interactions (Ledger and Hildrew 2005, Monteith et al. 2005). In addition to resistance to change, ecosystems can undergo transitions between states (regime shifts) induced by interactions between internal ecosystem processes and external fluctuations (Scheffer and Carpenter 2003). Scheffer and Carpenter (2003) used the expression catastrophic regime shifts for irreversible ecosystem changes. Hence, a number of internal and external factors, such as habitat quality, connectivity, dispersal abilities, and species interactions, can individually or collectively affect the rates and trajectories of biological recovery (Yan et al. 2003).
Climatic variability also influences the recovery of stream invertebrate assemblages. Wet winters increase acidity in Welsh (UK) moorland and forest streams and reset recovery (Ormerod and Durance 2009). Positive (i.e., warmer, wetter) North Atlantic Oscillation phases were accompanied by reduced interannual stability ( = similarity) in invertebrate assemblages (Durance and Ormerod 2007). Circumneutral streams were more sensitive to climatic variability than were acidified streams, a result suggesting that acidification stress could override or interact complexly with the stress arising from global warming (Christensen et al. 2006, Durance and Ormerod 2007). Thus, understanding the factors important for ecosystem recovery is difficult. Long-term data can be valuable for disentangling many of the factors driving or impeding recovery, and these types of data are of particular interest for determining the pathways and trajectories of systems undergoing recovery under the constraints set by climatic variability.
We used water-chemistry, phytoplankton, and littoral benthic invertebrate assemblage data from 4 minimally disturbed (reference) and 4 acidified lakes studied during the last 2 decades to track recovery pathways from acidification. We tested 4 specific hypotheses: 1) acidified lakes show chemical and biological recovery as indicated by increased pH and diversity (e.g., biological assemblages of acidified lakes are more dissimilar to those of reference sites at the start than at the end of the study period, and dissimilarity decreases with time), 2) trajectories of biological recovery vary between assemblages, 3) long-term assemblage trends, indicated by between-year shifts in assemblage composition (measured as Euclidean distance), are related to decreased acidity, and 4) climatic variability subsumes the recovery signal from acidification.
Methods
Study area
In the late 1980s, Sweden initiated long-term monitoring of multiple habitats and trophic levels of lakes to follow the effects of acidification and recovery of regionally representative lake ecosystems (Johnson 1999). Two lake categories were monitored: 1) minimally disturbed reference lakes (sensu Stoddard et al. 2006) to determine natural, interannual changes and 2) acidified lakes (defined as exceedence of critical load of S > 0) to assess the degree of degradation and subsequent recovery.
Samples for analysis of water chemistry and phytoplankton and littoral invertebrate assemblages have been collected during the last 21 y (1988–2008) from 8 lakes in the boreonemoral or mixed forest ecoregion of southern Sweden (Fig. 1). Four lakes with a mean annual minimum pH > 6.2 were considered minimally disturbed reference lakes, and 4 lakes with mean annual maximum pH < 5.8 were judged as being acidified based on estimates of exchange of preindustrial alkalinity for anthropogenic SO42− (Persson 1996). All 8 lakes are relatively small (mean lake surface area = 0.47 ± 0.52 km2), shallow (mean depth = 5.8 ± 2.5 m), and nutrient poor (mean annual PO4-P concentration <6.5 µg/L) (Appendix 1). Catchment land use/cover also is similar between reference and acidified lakes (% forest: mean = 66%, range = 54–85%; % mire mean = 11%, range = 4–20%). More information on these lakes is available from: http://www.ma.slu.se.
Figure 1
Map of Sweden showing the 8 study lakes and 3 main ecoregions: 14 (Central Plains), 22 (Scandinavian Shield), and 20 (Boreal Highlands) according to Illies (1966).
Sampling
Surface-water samples (0.5 m) were collected 6 to 8 times a year (monthly during the open-water phase and once in February) at a mid-lake station in each lake. Water was collected with a Plexiglas® sampler and kept cool during transport to the laboratory. Samples were analyzed for variables indicative of acidity (pH and SO42− concentration), nutrients (PO4-P), and water color (absorbance at 420 nm of filtered water). All physicochemical analyses were done at the Department of Aquatic Sciences and Assessment following international (ISO) or European (EN) standards when available (Wilander et al. 2003). Annual mean values were used in statistical analyses to down-weight seasonal differences.
Phytoplankton was sampled in August of each year by taking a water sample from the epilimnion (0–4 m) with a Plexiglas tube sampler (diameter = 3 cm). In lakes with a surface area >1 km2 a single mid-lake site was used for sampling. In lakes with a surface area <1 km2, 5 random epilimnetic water samples were taken and mixed to form a composite sample from which a subsample was taken and preserved with Lugol's iodine solution (2 g potassium iodide and 1 g iodide in 100 mL distilled water) supplemented with acetic acid. Phytoplankton counts were made using an inverted light microscope and the modified Utermöhl technique commonly used in the Nordic countries (Olrik et al. 1989). Taxa were identified to the lowest taxonomic unit possible (usually species), and species-specific biovolume measures were calculated from geometric formulae. In 1992, enumeration of phytoplankton changed from counts of predominating taxa to counts of all species in a sample. Therefore, phytoplankton counts made before 1992 are not included here.
Benthic invertebrates were collected from wind-exposed, vegetation-free littoral habitats in late autumn (October–November) each year. Five replicate samples were taken by standardized kick sampling with a hand net (0.5-mm mesh size). Each sample was taken by disturbing the bottom substratum for 20 s along a 1-m-long stretch of the littoral region at a depth of ∼0.5 m. Samples were preserved in 70% ethanol in the field and processed in the laboratory by sorting against a white background with 10× magnification. Invertebrates were identified to the lowest taxonomic unit possible (generally to species level) and counted with dissecting and light microscopes.
Climatic data
The effect of interannual variability in climate on phytoplankton and invertebrate assemblages was tested using the North Atlantic Oscillation (NAO) winter index (NAOwinter; December–March, provided by Climate Analysis Section, National Center for Atmospheric Research, Boulder, Colorado; http://www.cgd.ucar.edu/cas/jhurrell/nao.stat.winter.html; Hurrell 1995). The NAOwinter index is calculated from the difference in sea surface pressure between the Azores and Iceland. Positive values are associated with mild, wet winters in northwestern Europe and negative values are associated with cold, dry winters (Hurrell 1995). During the period 1988 to 1995, the NAOwinter index averaged +2.97 (min: +0.72 for 1988, max: +5.08 for 1989), compared to +0.32 for 1996 to 2008 (min = −3.78 for 1996, max = +2.80 for 2007). For southern Sweden, the NAOwinter index shows a switch from warm to cold winters between 1995 and 1996 (NAOwinter decreased from +3.96 to −3.78, respectively).
Data analysis
Biological response variables
Four measures of diversity and abundance, derived from the species × site data sets, were calculated using the PAST software (PAST - PAlaeontological STatistics, version 1.73; available from http://folk.uio.no/ohammer/past/index.html) to determine if phytoplankton and invertebrate assemblages are recovering from acidification: 1) total number of individuals (invertebrates) or biovolume (phytoplankton) in a sample, 2) total number of taxa in a sample (taxon richness), 3) Simpson diversity (Simpson 1949), and 4) evenness, calculated as the logarithm of Shannon diversity (H; Shannon and Weaver 1949) divided by the number of taxa (eH/S). In addition, 2 measures of assemblage composition were calculated using nonmetric multidimensional scaling (NMDS) on presence–absence and abundance (invertebrate) or biovolume (phytoplankton) data in PAST 1.73. NMDS was carried out on a similarity matrix built on presence–absence data (Sørensen similarity) and a dissimilarity matrix built on log(x)-transformed biovolume (phytoplankton) and abundance (invertebrate) (Bray–Curtis dissimilarity). Measures of both similarity/dissimilarity are recommended because they ignore joint absences, a combination that predominates in site × species matrices (Faith et al. 1987). Moreover, presence–absence data gives more weight to species with a reduced distribution, whereas abundance/biovolume data is more recommended for use where many taxa are widely distributed (Rabeni et al. 2002) (i.e., presence–absence data distinguished differences in species composition without considering the abundance/biovolume of each species).
Functional groups also were determined. Phytoplankton was divided into autotrophic, mixotrophic, and heterotrophic biomass groups based on the classification scheme of Jansson et al. (1996). Bacillariophyceae, Conjugatophyceae, Cryptophyceae, Cyanophyceae, Loxophyceae, Prasinophyceae, Xanthophyceae, and Chlorophyceae were considered to be functionally autotrophic, except the green algae Polytoma and Polytomella (mixotrophic taxa). Chrysophyceae, Craspedophyceae, Dinophyceae, Euglenophyceae, Haptophyceae, and Raphidophyceae were regarded as mixotrophic. Taxa were considered heterotrophic when they contained no photosynthetic apparatus (e.g., euglenozoan flagellates). Functional feeding groups for invertebrates (gatherers/collectors, parasites, xylophagous taxa, predators, miners, active and passive filter feeders, shredders and grazers/scrapers) were assigned with Asterics 3.1.1. (IRV Software, Vienna, Austria; http://www.fliessgewaesserbewertung.de/en/download/berechnung/).
Analysis of trends
Water quality and temperature
Regression was used to test the conjecture that water-chemistry variables are responding to decreased emissions of acidifying compounds (restoration from acidification). Preliminary analyses showed that several chemical variables were autocorrelated with time (year of sampling). Therefore, relationships between response variables and time were tested using Kendall's τ rank correlation coefficient (Kendall 1938), a nonparametric test of concordance (agreement between 2 rankings). The same test was used to evaluate the temporal trends of water temperature in the lakes. The Wilcoxon signed-rank test was used to detect differences in physicochemical variables between acidified and reference lakes (Wilcoxon 1945). Pearson product–moment correlation analyses were used to determine the relationship between the NAOwinter index and the water-quality variables.
Biological variables
Structural and functional metrics of phytoplankton and invertebrates were log(x)-transformed when necessary to improve normality and homoscedasticity before analyses with a mixed-model, repeated-measures analysis of variance (ANOVA) in Statistica (version 5; Statsoft, Inc., Tulsa, Oklahoma). The model contained treatment (acidified lakes vs reference lakes; fixed factor), year (sampling years between 1988 and 2008; fixed factor), and lakes nested within treatment (lakes[treatment]; random factor). Variance components were calculated for each term in the model to provide information beyond the significance level about the variation contributed by each term. The statistical design for ANOVA showing factors, multipliers, and formulas to determine estimated mean squares, variance components, and F-ratio tests are shown in Appendix 2.
The key term in this ANOVA was the interaction treatment × year which allowed us to test whether all acidified lakes converged with control lakes over the study period, a result that would indicate consistent recovery patterns. To avoid detection of false positives, significance levels were corrected using the procedure by Holm (1979), which is less conservative than the sequential Bonferroni correction.
Interannual changes in assemblage composition
NMDS was used to condense the species × site data sets into 2 measures (1st and 2nd axis scores) representing lake-specific between-year and among-lake similarities/dissimilarities in assemblage composition. NMDS axis scores for phytoplankton and invertebrate assemblages were calculated using Sørensen similarity and Bray–Curtis dissimilarity as described above. A distance matrix based on the individual lake–year NMDS results was calculated as the Euclidean distance (√[(xji – xi1)2 + (xi2 – xj2)2]) between the NMDS axis 1 and 2 scores. Between-year similarity/dissimilarity for pairs of successive years for each lake were extracted from the Euclidean distance matrix and used to quantify the degree of dissimilarity of assemblage composition between successive sampling years for the individual lakes and for the lake groups (acidified and reference lakes).
Recovery—movement towards reference conditions
The mean coordinates of individual lake–year combinations of the reference sites was used as a measure of the reference condition and to test the assumption that biological assemblages of acidified lakes became less dissimilar to those of reference lakes with time. The reference condition for phytoplankton and invertebrate assemblages was estimated as the mean of the NMDS axis 1 and 2 scores for all 4 reference lakes. The Euclidean distance between individual lake–year combinations of acidified and reference lakes and the mean of the reference lakes was extracted from a distance matrix of NMDS axis 1 and 2 scores calculated with Sørensen similarity and Bray–Curtis dissimilarity as described above. A Wilcoxon signed-rank test was used to test if between-year distance was lower for acidified than reference lakes. Kendall concordance of Euclidean distance of lake–year combinations to the reference population was used to determine if distance decreased in acidified lakes over time.
Environmental correlates of between-year shifts in assemblage composition
The Euclidean distance between individual lake–year combinations of acidified and reference lakes was extracted from a distance matrix of NMDS axis 1 and 2 scores calculated with Sørensen similarity and Bray–Curtis dissimilarity as described above. Spearman rank correlation analysis was used to determine the associations between shifts (between-year Euclidean distance) in phytoplankton and littoral invertebrate assemblages and water-chemistry and the NAOwinter index across lakes.
Results
Trends in surface-water chemistry
Acidified lakes had lower pH than reference lakes (Wilcoxon signed-rank test, p < 0.001). Mean annual pH was 5.3 ± 0.41 for the 4 acidified compared to 6.6 ± 0.17 for the 4 reference lakes (Appendix 1). Mean surface-water pH increased significantly with time in 5 of the 8 lakes during the 21-y period (Appendix 1). However, acidic episodes occurred during the recovery period. Annual minimum pH was <5.0 (mean: 4.4 ± 0.41) in the acidified lakes compared to >6.0 (mean: 6.2 ± 0.17) for the reference lakes. One lake (Härsvatten) has had pH < 5.0 during the 21-y of study, but the other 3 acidified lakes also have experienced acidic episodes. Annual mean pH < 5.5 was recorded in 8 (Övre Skärsjön) to 15 (Rotehogstjärnen) of the 21 study years. Significant decreases in SO4 concentration were noted for all 8 of the lakes. Likewise, increasing trends in lake water color were noted for all 8 of the lakes, with acidified lakes having higher water color. No temporal trends or differences were noted between acidified and reference lakes for PO4-P concentration.
Surface water temperatures increased significantly only in the reference lakes Fräcksjön and Stora Skärsjön (Appendix 1). A common pattern arising from the Pearson correlation analysis was that trends in water temperature were independent of the NAOwinter index (p > 0.05). However, correlations of the NAO with physicochemical variables were lake specific, and no common patterns were revealed for either reference or acidified lakes (Appendix 3).
Patterns of biological recovery
Recovery patterns based on both structural and functional metrics of phytoplankton (Fig. 2A–I) and invertebrates (Fig. 3A–J) were highly variable among acidified lakes.
Figure 2
Temporal trends of phytoplankton taxon richness (A), biovolume (B), Simpson Index (C), evenness (D), nonmetric multidimensional scaling (NMDS) axis 1 (Sørensen's similarity) (E), NMDS axis 1 (Bray–Curtis dissimilarity) (F), proportional biomass of autotrophs (G), heterotrophs (H), and mixotrophs (I) in reference lakes (represented by the mean and 95% confidence interval of 4 lakes: Allgjuttern, Fräcksjön, Stora Envättern, Stora Skärsjön) and individual acidified lakes (Brunnsjön, Härsvatten, Övre Skärsjön, Rotehogstjärnen) between 1992 and 2008.
Figure 3
Temporal trends of littoral invertebrate taxon richness (A), abundance (B), Simpson Index (C), evenness (D), nonmetric multidimensional scaling (NMDS) axis 1 (Sørensen's similarity) (E), NMDS axis 1 (Bray–Curtis dissimilarity) (F), and % abundance of grazers (G), shredders (H), collectors (I), and predators (J) in reference lakes (represented by the mean and 95% confidence interval of 4 lakes: Allgjuttern, Fräcksjön, Stora Envättern, Stora Skärsjön), and individual acidified lakes (Brunnsjön, Härsvatten, Övre Skärsjön, Rotehogstjärnen) between 1988 and 2008.
Phytoplankton
Species evenness (Fig. 2D) and biovolume-based assemblage composition (NMDS 1, Bray–Curtis; Fig. 2F) showed the most consistent recovery for Brunnsjön and Övre Skärsjön. These metrics were either within or close to the 95% confidence interval reflecting target conditions in neutral lakes. In contrast, phytoplankton in lake Rotehogstjärnen was most resistant to recovery. Species richness (Fig. 2A), biovolume (Fig. 2B), Simpson diversity (Fig. 2C), and species composition (NMDS, Sørensen; Fig. 2E) showed different degrees of variability approaching reference conditions. Recovery trends of functional groups were highly variable and were characterized by periods where these groups showed intermittent approaches to reference conditions (Fig. 2G–I). As a result of this strong variability between and idiosyncratic responses within lakes, most comparisons in the ANOVA were not significant, except the significant temporal variation detected for biovolume-based assemblage composition (NMDS 1, Bray Curtis; Table 1). None of the metrics had a significant treatment × year interaction, i.e., no consistent recovery occurred across all acidified lakes (Table 1).
Table 1
Analysis of variance results showing mean squares (MS), sources in variation (F ratios), and variance components (VC; %) for structural and functional metrics of phytoplankton assemblages. Significant terms are highlighted in bold. The remaining terms were not significant after Holm correction of p levels. NMDS = nonmetric multidimensional scaling. * = p < 0.05, ** = p < 0.01, *** = p < 0.001.
Invertebrates
Overall, the invertebrates showed a better recovery response than the phytoplankton. Most structural metrics indicated that acidified lakes were within or close to the 95% confidence interval representing reference conditions (Fig. 3A–D). Brunnsjön and Övre Skärsjön seemed to have a better recovery response, whereas Härsvatten had a poorer response. Both composition- and abundance-based assemblage similarity showed variable degrees of approaching assemblage structure in reference lakes (Fig. 3E, F). Recovery responses among functional feeding groups were variable among lakes and over time. Predators generally showed the weakest response, and only Brunnsjön or Rotehogstjärnen were consistently within or close to reference conditions (Fig. 3J). These patterns also were reflected by the collectors, although Härsvatten and Övre Skärsjön were close to reference conditions (Fig. 3I). Shredders and grazers showed the most consistent recovery patterns, reflected by an approach to reference conditions during the 2nd ½ of the study period (significant year term in the ANOVA; Fig. 3G, H, Table 2). Only shredders showed a consistent recovery across acidified lakes, indicated by the significant interaction term, which explained 82% of the variation in the model.
Table 2
Analysis of variance results showing mean squares (MS), sources in variation (F ratios), and variance components (VC; %) for structural and functional metrics of littoral invertebrate assemblages. Significant terms are highlighted in bold. The remaining terms were not significant after Holm correction of p levels. NMDS = nonmetric multidimensional scaling. * = p < 0.05, ** = p < 0.01, *** = p < 0.001.
Between-year shifts in assemblage composition
Mean between-year shifts (Euclidean distance) of assemblage composition supported the idea that changes were stronger in acidified than in reference lakes (Appendix 4). Changes in phytoplankton and littoral invertebrate assemblage composition were significant when measured with Sørensen similarity (presence–absence data) (Wilcoxon signed-rank test, p < 0.05), but differences were not significant when biovolume or abundance data were used (Bray–Curtis dissimilarity). The strongest differences were noted for phytoplankton assemblages. Acidified lakes had a mean Euclidean distance of almost 2× that of reference lakes (acidified: 0.037 ± 0.039, reference: 0.019 ± 0.016; Sørensen similarity). Between-year shifts in littoral invertebrate assemblages were 0.035 ± 0.024 for acidified and 0.027 ± 0.025 for reference lakes (Sørensen similarity).
Distance to reference conditions
Distance to the mean of reference lakes was on average 2.8× (invertebrates, Bray–Curtis) to 4.0× (phytoplankton, Sørensen) greater for individual acidified lakes than for individual reference lakes. On average, phytoplankton assemblages in acidified lakes were 0.134 ± 0.044 (Sørensen similarity) to 0.113 ± 0.038 (Bray–Curtis dissimilarity) Euclidean distance units from the mean of the 4 reference lakes (NMDS axis 1 and axis 2 scores) (Table 3). Similar mean distances were found for littoral invertebrate assemblages (0.100 ± 0.054 for Sørensen and 0.101 ± 0.060 for Bray–Curtis). Härsvatten was the furthest from the reference lakes (Euclidean distance >0.152). Between-year distance of individual lakes to mean coordinates of the 4 reference lakes decreased significantly with time in all 4 acidified lakes, indicating recovery in terms of assemblage composition (Kendall concordance, p < 0.05) (Fig. 4A–H). However, interpretation of recovery was dependent on the taxonomic group evaluated. For example, phytoplankton assemblages of Härsvatten and Rotehogstjärnen showed marked changes over time (Härsvatten: Sørensen τ = −0.441, p = 0.014, Bray–Curtis τ = −0.706, p < 0.0001, Fig. 4F; Rotehogstjärnen: Sørensen τ = −0.103, p = 0.564, Bray–Curtis τ = −0.529, p = 0.003, Fig. 4H), whereas invertebrate assemblage dissimilarity did not decrease with time (Fig. 4B, D). In contrast, invertebrate assemblages of 2 lakes (Brunnsjön and Övre Skärsjön) had significant negative slopes during the 21-y study period (Brunnsjön: Sørensen τ = −0.162, p = 0.305, Bray–Curtis τ = −0.457, p = 0.004, Fig. 4A; Övre Skärsjön: Sørensen τ = −0.514, p = 0.001, Bray–Curtis τ = −0.571, p = 0.0003, Fig. 4C).
Figure 4
Plots of Euclidean distance of acidified lakes Brunnsjön (A, E), Härsvatten (B, F), Övre Skärsjön (C, G), and Rotehogstjärnen (D, H) to the mean of the reference condition of 4 minimally disturbed lakes against year. Euclidean distance was based on nonmetric multidimensional scaling (NMDS) axes 1 and 2 coordinates for Bray–Curtis dissimilarity and Sørensen similarity of presence–absence of invertebrate (A, B, C, D) and phytoplankton (E, F, G, H) assemblages. Statistics are shown for significant (p < 0.05) regressions.
Table 3
Mean (±1 SD) Euclidean distance of individual lake–year combinations to mean coordinates of nonmetric multidimensional scaling (NMDS) axis 1 and axis 2 scores of 4 reference lakes. Euclidean distance was calculated using NMDS on Sørensen similarity and Bray–Curtis dissimilarity for phytoplankton (1992–2008) and invertebrate (1988–2008) assemblages of 4 reference and 4 acidified lakes. Different letters in parentheses show significant differences (p < 0.05) between reference and acidified lakes using a Wilcoxon signed-rank test.
Environmental drivers of between-year shifts in assemblage composition
Between-year distance of phytoplankton assemblages was negatively correlated with lake surface-water temperature (−0.336) and pH (−0.318) (Sørensen similarity) and water color (−0.307; Bray–Curtis dissimilarity) (Table 4). Between-year distance of invertebrate assemblages was negatively correlated with temperature (−0.249), pH (−0.281) (Sørensen similarity), and N concentration (−0.167; Bray–Curtis dissimilarity), but positively related to NAOwinter index (0.249; Sørensen similarity).
Table 4
Spearman correlation (ρ) of between-year Euclidean distance of phytoplankton and littoral invertebrate assemblages in 4 acidified and 4 reference lakes and selected water chemistry variables (annual means between 1988 and 2008) and the North Atlantic Oscillation winter index (NAOwinter). Sørensen similarity and Bray–Curtis dissimilarity were used in a nonmetric multidimensional scaling ordination. * = p < 0.05, ** = p < 0.01, *** = p < 0.001
Discussion
Decreasing trends in surface-water acidity are consistent with results from several previous studies showing chemical recovery from acidification across Europe and North America (Stoddard et al. 1999, Eshleman et al. 2008). However, despite widespread evidence of chemical recovery, support for biological recovery has been patchy, and the environmental variables driving or delaying recovery have been difficult to ascertain. These results imply that response of acid-sensitive species (and assemblages) to decreased acidity is more complex in the recovery phase than in the impairment phase. A number of hypotheses have been raised to explain discrepancies between the impairment and recovery phases of acidification (Yan et al. 2003, Ledger and Hildrew 2005, Monteith et al. 2005). Briefly, recovery can be prevented or impeded by impoverished water quality (Lepori et al. 2003, Kowalik et al. 2007), biotic interactions (Ledger and Hildrew 2005), or scarcity of acid-sensitive colonists (Monteith et al. 2005). Choice of response variable also can affect ability to detect recovery. Clear shifts in lake assemblages of phytoplankton (Findlay 2003), zooplankton (Frost et al. 2006), and epilithic algae (Battarbee et al. 1988, but see Vinebrooke et al. 2003) have been related to decreased acidity, whereas results for other groups, such as fish and benthic invertebrate assemblages, have been more equivocal (Burns et al. 2008, Stendera and Johnson 2008). Our study further highlights the complexity of biological recovery responses to acidification. This complexity is reflected in strongly idiosyncratic changes over time as a function of the response variable chosen (functional vs structural, univariate vs multivariate) to track change for both assemblages.
Quantitative and qualitative aspects of recovery
Our results were based on a complementary analysis approach that allowed identification of trends of convergence of acidified lakes with reference lakes (quantitative recovery aspect) and of temporal variability shown by individual lakes during recovery (qualitative aspect). Our study showed that phytoplankton and invertebrate assemblage composition (Sørensen similarity) tended to converge quantitatively with those in reference lakes. The rate of recovery varied among lakes (indicated by the negative slopes in Fig. 4), but the result suggests that acidified lakes tend to regain sets of species present in reference lakes. This result supports the findings of other studies that have documented recovery of phytoplankton (Findlay 2003, Stendera and Johnson 2008) and invertebrate assemblages from acidification (Raddum et al. 2001, Lento et al. 2008, Stendera and Johnson 2008). However, convergence with reference conditions was equivocal with abundance-based measures of similarity (Bray–Curtis dissimilarity). Invertebrate assemblages in Härsvatten and Rotehogstjärnen and phytoplankton assemblages in Övre Skärsjön tended to diverge from reference conditions during recovery. This result suggests that recovery patterns of assemblage composition and dominance patterns are not harmonic and that internal abiotic and biotic lake characteristics (Yan et al. 2003) individually or collectively influence population frequencies during recovery.
Recovery of phytoplankton and invertebrate assemblage composition (Sørensen similarity) rather than abundance-based assemblage structure (Bray–Curtis dissimilarity) underlines the importance of considering species dominance patterns for understanding ecological patterns and processes (Hillebrand et al. 2008). Species distributions tend to become more homogenized regionally as a result of recolonization of acid-sensitive biota, a pattern consistent with landscape-level recovery from acidification (DGA, RKJ, and Willem Goedkoop, Swedish University of Agricultural Sciences, unpublished data). However, local lake conditions contribute to species sorting and mediate the degree to which recolonizing species can become numerically dominant in the assemblages. Consequently, the recovery signal from patterns of species distribution can be obscured when assessing evenness, which can be mediated through biological interactions (Ledger and Hildrew 2005). Thus, our results provide mixed support of our 1st hypothesis (decreased dissimilarity between acidified and reference lakes over time). Compositional similarity data support it, whereas abundance-based dissimilarity data do not. However, the results do support our 2nd hypothesis (recovery trajectories differ between assemblages as a function of the response variable used to detect change). This finding agrees with those of earlier studies (e.g., Stendera and Johnson 2008) and shows that inference about biological recovery can be biased if results are based on single-assemblage studies. Comparative assessment approaches using multiple assemblages are needed to evaluate integral ecological responses to environmental stress (Johnson and Hering 2009).
We acknowledge that quantitative assessment of assemblage recovery, expressed as the Euclidean distances of acidified lakes to the mean of reference lakes between sampling years is a static approach that underestimates temporal variability in the reference lakes. Reference lakes also can undergo environmental change over time, particularly when broad-scale environmental drivers (e.g., climate) cause coherent changes across aquatic ecosystems at macroecological scales (Li et al. 2006, Kent et al. 2007, Fisher et al. 2008). Accordingly, defining convergence of impacted sites on a temporal mean of reference sites could either over- or underestimate recovery when the reference sites also show directional temporal change. Assemblage structure in our reference lakes showed no marked directional change despite variability over the study period (Figs 2A–I, 3A–J). Thus, our quantitative recovery estimates can be viewed with confidence.
Our results provide quantitative evidence of recovery and evidence for qualitative differences between lakes and assemblages. Thus, they highlight the importance of addressing short- and long-term temporal dynamics when assessing recovery. Acidified lakes had greater temporal variability than reference lakes, as has been observed in other stressed environments (Angeler and Moreno 2007), but variability in temporal patterns was assemblage and lake specific. As a result, a consistent convergence of structural and functional metrics in all acidified lakes with reference conditions was not observed. Instead, each acidified lake showed specific periods of approximation to and convergence with reference conditions, punctuated by periods during which stressed lakes deviated from these desired states. This pronounced variability can be partially explained by the environmental variables that were correlated with this temporal variability.
Environmental correlates of recovery
Our results provide mixed support for our 3rd hypothesis (biological recovery, measured as decreases in between-year Euclidean distance, is a result of decreased acidity). The effects of decreased acidity were evident in species composition of assemblages (Sørensen similarity) but not in species abundances (Bray–Curtis dissimilarity) in both phytoplankton and invertebrate assemblages. Stendera and Johnson (2008) showed that changes in phytoplankton taxon richness in these boreal lakes were related not only to increased pH but also to other environmental variables, such as increased water color (a proxy for dissolved organic C [DOC]) and nutrient concentrations. Our study showed that high between-year shifts in phytoplankton assemblage composition were negatively correlated with temperature, water color, and pH. Increased pH should permit colonization of acid-sensitive species, but concurrent increases in water color (DOC concentrations) and temperature also could have influenced assemblage structure (Wetzel 1995). Interannual shifts in invertebrate assemblage composition also were negatively correlated with temperature and pH. In addition, NO2 + NO3 concentration and the NAOwinter index were significantly correlated with temporal assemblage turnover.
The imprint of climate on recovery
The NAO is an important regional-level driver of climate in southern Sweden (Weyhenmeyer 2004). Positive NAO values signify higher precipitation and, thus, more variable hydraulic and chemical conditions. Invertebrates seemed to respond more readily to the NAOwinter index than did phytoplankton assemblages. Life-history characteristics mediate the strength of biological responses to climate forcing (Adrian et al. 2006). Much shorter generation times of phytoplankton relative to invertebrates could uncouple phytoplankton responses from climate signals and could make these responses visible only via indirect or integrated effects (Ottersen et al. 2001). These effects could be manifested through responses to abiotic or biotic variables that are more strongly affected by the NAO. Rusak et al. (2008) showed climate imprints on zooplankton dynamics, and these effects could have at least some influence on phytoplankton. We acknowledge that using single annual summer samples in our analyses could have caused us to miss potential NAO effects on spring phytoplankton development (Weyhenmeyer et al. 1999), but our results suggest that this potential NAO effect does not influence phytoplankton during later successional stages. This conclusion is in agreement with the results of a recent meta-analysis of data from European lakes, which showed no correlations between phytoplankton assemblages in all seasons with the NAO (Blenckner et al. 2007).
Between-year Euclidean distance of invertebrate assemblage composition was correlated positively with the NAOwinter index across lakes. Durance and Ormerod (2007) found that the imprints of the NAO on stream invertebrate assemblages were particularly strong in circumneutral streams. They suggested that acidification stress in streams could override a climate signal, possibly because of the complex interaction of the natural and anthropogenic disturbance regimes to which stream organisms are exposed under the constraints set by climatic conditions. The difference between the study of Durance and Ormerod (2007) and ours could be because lake ecosystems are much less affected by changes in surface-water hydrology than are stream ecosystems, and hence, stream invertebrate assemblages respond more readily to climatic variability regardless of acidification status.
The significant correlations of both assemblages with temperature also highlight the importance of climatic variables in shaping among-year shifts in assemblage composition. Correlation of between-year Euclidean distance of invertebrate assemblages and NAOwinter index values were significant for acidified lakes, and greater between-year distance of invertebrate assemblages was associated with warmer winters (i.e., positive NAO values). Several studies have shown how interannual variability in NAO affects lake assemblages by affecting temperature, ice cover, and timing of spring algae blooms (Weyhenmeyer et al. 1999). In the boreal lakes of our study, the effect of broad-scale climatic variability might have context-dependent effects on water quality which, in turn, can mediate shifts in invertebrate and phytoplankton assemblages. Thus, recovery might be offset by effects of interannual, climate-driven changes on local lake variables (but see Evans et al. 2006). This supports our 4th hypothesis (climatic effects subsume the signal of biological recovery). In addition to climate-mediated imprints on recovery, a further plausible explanation for the marked interannual variability in assemblage structure is that acidic episodes can inflict recurrent stress on lake assemblages. This conjecture is supported by data showing that annual mean pH was <5.5 during several of the 21 y of study in the acidified lakes.
Predisturbance conditions
An apparent weakness of our study is that predisturbance measures of water chemistry and phytoplankton and littoral invertebrate assemblages are not available. However, palaeoreconstruction of preindustrial pH indicates that all 4 lakes were relatively well buffered prior to acidification (diatom-inferred pH ranged from 6.0 [Rotehogstjärnen] to 6.7 [Övre Skärsjön]) (Guhrén et al. 2003, Erlandsson et al. 2008). These preindustrial estimates of pH are close to the range of present-day measures of pH in the 4 reference lakes (mean pH ranged from 6.4 to 6.8) and indicate that water quality of the 4 reference lakes is a reasonable estimate of the predisturbed state. Furthermore, the finding that phytoplankton and littoral invertebrate assemblage composition are moving toward the reference population lends support to the conjecture that the acidified lakes show signs of approaching their predisturbed abiotic state.
Conclusion
Our results are of general interest concerning recovery of ecosystems from disturbance, but are particularly relevant to lake management. Phytoplankton and littoral invertebrate assemblages responded differently to improved water quality among acidified lakes. Thus, investigators should consider using multiple groups of organisms when a priori knowledge of response (recovery) signatures is poor. Use of phytoplankton or invertebrate assemblages alone would have led us to a biased conclusion regarding biological recovery. Moreover, the choice of the response variable influenced our ability to detect recovery. Assemblage composition (NMDS ordination) was powerful for detecting changes, but the dissimilarity measure used (whether based on presence–absence or biovolume/abundance) influenced the results. One of our most interesting findings was the congruence of trends in water chemistry and biology in acidified and reference lakes. Our results show the influence of factors other than decreased acidity, such as the underlying importance of climate, as drivers of recovery trajectories and pathways.
Acknowledgments
We thank the many people involved in collecting and processing the physicochemical, phytoplankton, and invertebrate samples. In particular, we thank Björn Wiklund and Putte Olsson for sorting and Lars Eriksson for identifying the invertebrates, and Eva Herlitz, Ann-Marie Wiederholm, and Isabel Quintana for identifying phytoplankton. Financial support for this project was partly provided by the EU funded projects EURO-LIMPACS (contract no. GOEC-CT-2003-5055), WISER (contract no. 226273), and REFRESH (contract no. 244121). The Swedish Environmental Protection Agency is acknowledged for making these data available. We thank S. Ormerod, M. Barbour, and an anonymous referee for helpful comments on a previous version of the manuscript.
Literature Cited
Appendices
Appendix 1.
Mean (±1 SD) values of water-chemistry and temperature variables and Kendall τ rank correlation coefficients of water-chemistry and temperature variables against year for 4 reference and 4 acidified lakes sampled from 1988 to 2008. Different letters in parentheses shows significant differences (p < 0.05) between reference and acidified lakes. Only significant regressions are shown. * = p < 0.05, ** = p < 0.01, *** = p < 0.001.
Appendix 2.
Statistical design for analysis of variance showing factors, degrees of freedom, multipliers, and formulas to determine estimated mean squares (E[MS]), variance components, and F-ratio tests. Factors and levels: Treatment (T) has i = …a levels, Lake(treatment) (L[T]) was replicated n times (r = …n) and served as the error term for the spatial variation, Year (Y) has j = …b levels, all combinations of the temporal variation effect were replicated N times (s = …N), and the term Y × L(T) is the error term for the temporal variation. Two factors are orthogonal and 1 is nested: T and Y are fixed, and L nested in T (L[T]) is random. Multipliers and formulas to determine E(MS) and variance components, denominators for F-ratio tests (den.; numbers correspond to numbers in the sources of variation column), and degrees of freedom (df) are shown. σe2, error term for between-sample variation; σE2, error term for within-sample (repeated measures) variation.Models: between samples: Yir = μ + Ti + L(T)r(i), within samples: Yijr = μ + Yj + TYij + YL(T)jr(i)
Appendix 3.
Correlation coefficients from Pearson correlation relating the North Atlantic Oscillation winter [NAOwinter] index to water-quality variables. Correlations that were not significant (p > 0.05) are not shown.* = p < 0.05, ** = p < 0.01.
Appendix 4.
Mean (±1 SD) between-year shifts (Euclidean distance) of phytoplankton and invertebrate assemblages in 4 reference and 4 acidified lakes. Euclidean distance was calculated using NMDS on Sørensen similarity and Bray–Curtis dissimilarity for phytoplankton (1992–2008) and invertebrate (1988–2008) assemblages. Different letters in parentheses show significant differences (p < 0.05) between reference and acidified lakes using a Wilcoxon signed-rank test.
