Effects of Climate , Limnological Features and Watershed Clearcut Logging on Long-Term Variation in Zooplankton Communities of Boreal Shield Lakes

In Canada, climate change and forest harvesting may both threaten the ecological integrity of boreal lakes. To disentangle the effects of natural variation in climate and lake environments from those of logging, we evaluated long-term variation (1991–2003) in zooplankton communities of six boreal lakes in Ontario. We monitored concomitantly changes in zooplankton abundance and composition in three undisturbed and three harvested lakes, five years prior and eight years after watershed clearcut logging. We tested the hypothesis that long-term natural variation in climate and lake environments will be more important drivers of zooplankton community changes than short-term impacts of logging. We used space/time interaction tests and asymmetric eigenvector maps to model zooplankton responses to environmental changes and logging. Year-to-year variation in zooplankton abundance and composition were almost an order of magnitude whereas among-lake variation was stable through time. Breakpoints in time series of zooplankton in each lake were not directly related to logging. Climatic and limnological features were the most important drivers of long-term variation in the zooplankton community, shading the effect of logging. These results highlight the need to better understand the pressures exerted by climate change on boreal lake ecosystems in the context of anthropogenic pressure, such as logging.


Introduction
Lakes are the major landscape features of the boreal ecoregion in Canada.Climate change and watershed clearcut logging are two major disturbances in Canadian forests and represent important threats for the ecological integrity of boreal lakes [1,2].Climate warming in boreal regions (increasing temperature and irradiance) can strongly affect physical (ice cover, thermocline depth), chemical (water quality), and biological (plankton productivity) characteristics of lakes across large spatial scales [2,3].Although lakes are recognized as sentinels of climate change [4][5][6] and watershed land-use [7,8] at regional scale, their responses to climate fluctuations may be confounded by multiple watershed anthropogenic disturbances [9,10].Inversely, climatic conditions, especially year-to year differences in snow accumulation and rain precipitation, modulate the impacts of watershed disturbances by clearcut logging in boreal lakes [1,11,12].Complex interactions between natural variation in climate and anthropogenic watershed disturbances make it difficult to disentangle their specific effects on ecological communities and functions of boreal lakes.In particular, the scarcity of long-term, broad scale (watershed) studies has limited our understanding of the combined effect of climate change and watershed human disturbances on plankton communities and lake functioning [13][14][15].
Since 2000, comparative and experimental studies have documented the impacts of watershed clearcut logging on biogeochemical cycles and ecological processes of boreal lakes in Canada [1,11,12].Forest harvesting can cause minor to strong increases in nutrient and organic carbon inputs from catchments to lakes, which subsequently affect water quality and transparency, lake productivity and trophic food webs [1,16,17].An increase in dissolved organic carbon (DOC) changes the light regime in the euphotic zone and reduces water transparency [18,19].Higher inputs of phosphorus and nitrogen (TP and TN) increase algal biomass (Chlorophyll-a) and primary productivity [17,20].Generally, the impacts of watershed logging on water quality of boreal lakes are enhanced by high precipitation and runoff [12].However, most of these impacts are short-term (<3 years) and have subtle bottom-up effects on zooplankton [12,17,[21][22][23] and fish communities [24][25][26][27][28].A slight increase in crustacean zooplankton density has been observed in Finnish lakes after watershed logging [18,29].Both calanoid copepod and cladoceran biomass declined in Alberta lakes after partial logging of the watershed [11].In Québec lakes, rotifer biomass increased whereas calanoid biomass decreased at short term after clearcut logging [21,22].However, changes in zooplankton species richness, and composition were negligible and not clearly related to logging disturbance [23].Scaling up in food webs, post-harvest enhancements in nutrients, and dissolved organic carbon likely favor early growth of yellow perch due to better feeding conditions on zooplankton prey [28,30], only slightly affect littoral fish [25], and do not influence the habitat of the lake trout [24] in Boreal Shield lakes.Watershed logging increases the depletion of oxygen under ice cover in shallow lakes of the Boreal Plain and can enhance fish winterkill [27].
As zooplankton organisms respond both to bottom-up and top-down forces (multiple force hypothesis [31]), it is very difficult to disentangle the effects of natural variation in abiotic (climate, water quality, nutrients) and biotic (algal resources, grazing, and predation) factors from the direct and indirect effects of logging disturbances.These difficulties arise from the confounding effects of natural variation at regional scale in climate, limnological features, and trophic interactions.Most studies assessing the impacts of watershed forest harvesting on boreal lakes did not measure the interactions between natural environmental variation and harvesting disturbance prior and after logging.Thus, difficulties in the interpretation of results from non-replicated experimental design remain an ongoing issue [32].
To investigate long-term impacts of watershed forest harvesting on Boreal Shield lakes in regards to natural environmental variation in climate and limnological features, we monitored changes in zooplankton communities of three harvested lakes five years before and eight years after logging, concomitantly with three neighbouring undisturbed lakes.Long-term monitoring (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) of preand post-disturbance variations in climate, water quality, and zooplankton communities allowed us to contrast spatio-temporal variation in zooplankton under two scenarios: undisturbed watersheds vs partially logged watersheds.We investigated the stability of zooplankton community through time in undisturbed and harvested lakes to determine breakpoints in zooplankton space/time variation and to relate it to logging disturbance.We tested the hypothesis that long-term space/time natural variability in the climate and lake environments will have more influence on the variation of zooplankton abundance and composition than the short-term impacts of logging disturbance.A previous study based on a comprehensive lake ecosystem model (LakeWeb-model) indicated that it is often realistic to assume that logging activities do not strongly influence lake community structures [33].
Our study is the first to assess long-term variation (>10 years) in the zooplankton community in oligotrophic lakes subjected to logging disturbance in Canada's Boreal Shield.Other Canadian studies used short-term (<3 years) comparative [21][22][23] or experimental [11,17,30] designs.Since climate fluctuations and limnological conditions may induce considerable, but natural variation in zooplankton community [34], impact assessment of watershed disturbances must be replicated at the scale of the ecosystems [13].The final aim should be to test if the change in lake community structure and function is beyond what one may simply expect from natural variation in unperturbed ecosystems [14,35,36].C in January.Annual precipitation was on average 609 mm rain and 171 cm snow, and the snow cover generally extended from November to May (1990-2003 at Mine Centre station, station 6025203, Environment and Climate Change Canada).Lake watersheds were covered by shallow soils with abundant boulders and bedrock outcrops; they are located at an elevation of approximately 427 m, have moderate slopes (15-50%) and relatively flat topography [16].Forest cover consisted of jack pine (Pinus banksiana), trembling aspen (Populus tremuloides), and black spruce (Picea mariana), with some paper birch (Betula papyrifera), eastern white cedar (Thuja occidentalis), red pine (Pinus resinosa), and eastern white pine (Pinus strobus).

Materials and Methods
zooplankton community [34], impact assessment of watershed disturbances must be replicated at the scale of the ecosystems [13].The final aim should be to test if the change in lake community structure and function is beyond what one may simply expect from natural variation in unperturbed ecosystems [14,35,36].

Study Area, Multiple Before/After-Control-Impact (MBACI) Experimental Logging Design and Limnological Features
The study area was in the Coldwater Lakes Experimental Watersheds in the Boreal-Great Lakes transition forest on the Canadian Shield (49°07′16′′-40°08′05′′N; 92°12′22′′-92°13′34′′W), approximately 200 km northwest of Thunder Bay, Ontario, Canada (Figure 1).The climate was characteristic of the north-temperate boreal ecoregion, with mean air temperatures of 18 °C in July and −15 °C in January.Annual precipitation was on average 609 mm rain and 171 cm snow, and the snow cover generally extended from November to May (1990-2003 at Mine Centre station, station 6025203, Environment and Climate Change Canada).Lake watersheds were covered by shallow soils with abundant boulders and bedrock outcrops; they are located at an elevation of approximately 427 m, have moderate slopes (15-50%) and relatively flat topography [16].Forest cover consisted of jack pine (Pinus banksiana), trembling aspen (Populus tremuloides), and black spruce (Picea mariana), with some paper birch (Betula papyrifera), eastern white cedar (Thuja occidentalis), red pine (Pinus resinosa), and eastern white pine (Pinus strobus).The study design was analogous to a MBACI (multiple before/after-control-impact) design [37,38] including three undisturbed lakes (L20, L38, and L80) and three watershed-logged (harvested) lakes (L26, L39, and L42) monitored from 1991 to 2003 (1993-2003 for L38).Some of the study lakes were hydrologically connected (Figure 1).Lake L26 was a 1st order lake that flows in L20 (2nd order); Lake L42 was also a 1st order lake upstream of L39 (2nd order) which flowed into L38 (3rd order); L80 was a 1st order lake not hydrologically connected to the other studied lakes.Due to forest industry planning, the intensity of clearcut logging varied among lakes.Watersheds of L26, L39, and L42 were subjected to partial clearcutting (33%, 62%, and 71%, respectively) during summer 1996 The study design was analogous to a MBACI (multiple before/after-control-impact) design [37,38] including three undisturbed lakes (L20, L38, and L80) and three watershed-logged (harvested) lakes (L26, L39, and L42) monitored from 1991 to 2003 (1993-2003 for L38).Some of the study lakes were hydrologically connected (Figure 1).Lake L26 was a 1st order lake that flows in L20 (2nd order); Lake L42 was also a 1st order lake upstream of L39 (2nd order) which flowed into L38 (3rd order); L80 was a 1st order lake not hydrologically connected to the other studied lakes.Due to forest industry planning, the intensity of clearcut logging varied among lakes.Watersheds of L26, L39, and L42 were subjected to partial clearcutting (33%, 62%, and 71%, respectively) during summer 1996 after five years of pre-harvesting sampling; additional logging was carried out in summer 1998 on the watersheds of L26 and L39 (12% and 16%, respectively).Thus, logging disturbance was moderate in L26 catchment (in total 45%) and more intense in L39 (78%) and L42 (71%) catchments.Shoreline disturbances also varied among lakes; L26 had no shoreline disturbance with either no cutting at all or an undisturbed shoreline buffer strip, proportional to shoreline slope, of 30-90 m in width.L39 and L42 had more extensive shoreline disturbance (62% and 42%, respectively).The watersheds of the three other lakes (L20, L38, L80) remained unlogged but due to hydrological connectivity among lakes (Figure 1), two of them (L20, L38) may have been indirectly impacted by logging disturbance on the watersheds of headwater lakes (L42, L39, L26).Therefore, as the experimental design was not perfect, the six studied lakes exhibited a gradient in the intensity of logging disturbance, with lakes directly impacted (L26, L39, L42), lakes slightly impacted through their hydrologic connectivity (L20 and L38), and a lake (L80) not impacted.
Lakes showed some variation in morphometry and water quality conditions (Table 1).The undisturbed lakes tended to be shallower than the harvested lakes, but had similar mean depth.On average, lakes were of medium size, but the undisturbed lakes tended to be bigger than harvested lakes.Drainage basin areas of undisturbed lakes were also larger than these of harvested lakes.Average water renewal time (WRT) was seven years and ranged from 1 to 13 years.The undisturbed lakes had a shorter WRT than the harvested lakes (see also Supplementary Materials Table S1, for morphological features in each lake).Water quality was characteristic of dilute, slightly acidic, and clear waters of Boreal Shield lakes (Table 1).The undisturbed lakes were slightly more enriched in nutrients, dissolved organic carbon, and Chlorophyll-a (Chl-a), and had higher alkalinity.In contrast, the harvested lakes showed lower concentrations of nutrients, dissolved organic carbon, and Chlorophyll-a, as well as lower alkalinity (see also Supplementary Materials Table S2, for water quality in each lake).

Zooplankton Sampling and Analysis
Zooplankton was collected from 1991 to 2003 (1993-2003 for L38) using a Wisconsin plankton net (11.8 cm of diameter, 0.01 m 2 of area, 80 µm mesh size) vertically hauled from 1 m above sediments to the water surface at the deepest point of each lake.Generally, sampling was carried out every month during winter and every two weeks in May, June, July, August, September, and October.However, sampling intensity was irregular among years and lakes.Thus, for this study, we used only samples collected twice per month during the ice-free season (May to October).Collected zooplankton organisms were preserved in 4% formalin solution.
To determine spatio-temporal changes in zooplankton community composition among lakes and years, zooplankton species were aggregated into four taxonomic groups: Rotifera, Cladocera, Calanoida, and Cyclopoida.Abundance of taxa in each taxonomic group was estimated as densities per liter (Ind.•L−1 ) based on numbers of organisms counted, volume of the subsample, volume of the concentrated sample, and volume of water filtered at each sampling date and in each lake.Densities estimated during the ice-free season were averaged for each lake and year.

Natural Environmental Heterogeneity and Logging Disturbance Impact
To evaluate long-term natural environmental variation at regional scale, we selected a group of variables encompassing space/time variation in climatic and limnological features, able to drive changes in zooplankton communities.We considered year-to-year (time) variation in three major climatic variables, i.e., mean May-October air temperature ( • C), cumulative snow (cm) in the previous winter, and cumulative rain (mm) between May-October standardized to unit of watershed area.We considered among-lakes (space) variation in morphometry and hydrology based on three important variables: surface, maximum depth, and water residence time.To quantify among-lakes and among-years variation (space/time) in water quality, we used several physico-chemical variables (alkalinity, K, Mg, Ca, Chl-a, DOC, pH, TON, TP, and SO 4 , see Table 1).We used among-years variation in the percentages of clearcut logging on watershed area in the undisturbed and harvested lakes to quantify variation in the intensity of logging disturbance.

All Lakes Grouped
To assess long-term natural variation in climatic and limnological features, we performed principal component analyses (PCA) to quantify temporal variation in regional climate, spatial variation in morphometry of lakes, and water quality, using the selected variables (listed previously).We scaled variables to unit variance to allow comparison of factors [53].Additionally, we investigated the stability of climate throughout the study period with a time-constrained cluster analysis based on a multivariate regression tree analysis [53].This method allows the identification of one or several breakpoints in a data series [54].The optimal size of the tree was determined by minimising the cross-validation relative error [55].
Two Groups: Undisturbed and Harvested Lakes Because our experimental study did not follow a perfect MBACI design, the impact of clearcut logging on zooplankton abundance and community composition could not have been analysed by a simple three-way ANOVA.Instead, we performed two-way ANOVA without replication [56] on the harvested (L26, L39, and L42) and the undisturbed (L20, L38, and L80) groups of lakes separately.With this procedure, we tried to detect an interaction between lakes (space) and years (time).A significant interaction in the harvested group of lakes would be indicative that spatial pattern (among-lakes) in zooplankton community has changed through time (years) and may suggest an impact of logging disturbances, or conversely that the temporal variations differed among lakes.No significant interaction was expected in the group of undisturbed lakes.Briefly, this method consists of representing space and time by spatial and temporal eigenfunctions in two-way ANOVA to reduce the degree of freedom required [54].Specifically, we used Moran's eigenvector maps (MEM) to model among-years variation and Helmert contrast to represent the among-lakes variation since virtually no degree of freedom gain was possible because of the small number of lakes per group (n = 2-3).

Influence of Natural Environmental Variation and Clearcut Logging on Zooplankton Community
To test the effects of temporal (among-years: AEM, see below) and spatial (among-lakes: Lake, see below) patterns of variation in environmental factors, including logging disturbance, on zooplankton community, we used redundancy analysis (RDA) and variation partitioning [57] on the harvested and undisturbed groups of lakes separately.Zooplankton abundance data were transformed using the Hellinger transformation (i) to prevent similarity to reflect double zeros; (ii) to eliminate differences in total zooplankton abundance among years and, because it includes a square-root transformation; and (iii) to reduce the importance of the most abundant zooplankton groups, as recommended for RDA [58].Temporal structure was modelled with asymmetric eigenvector maps (AEM) [59].AEM produce a group of orthogonal temporal variables from the directional connexion between years (time) [54].Forward selection was applied separately to the environmental and AEM variables to reduce the number of variables before RDA modelling.

Analysis of Zooplankton Stability through Time in Individual Lakes
Since among-lake variation (Lake) was significant for both the harvested and undisturbed groups of lakes, we investigated the stability through time of total zooplankton abundance in each lake by detecting breakpoints with a time constrained cluster multivariate regression tree [54].We also evaluated how temporal pattern of variation in environmental variables and logging disturbance influenced zooplankton community composition for individual lake using RDA.Temporal structure was modeled at two levels: (i) linearly using a regression with the sampling years and (ii) non-linearly with asymmetric eigenvector maps (AEM).
All statistical analyses were made using the statistical software R [60].PCA analyses were done with FactoMineR library [61].Multivariate regression tree to detect temporal clusters were implemented with mvpart library.ANOVA for space-time interaction in absence of replication was implemented with the STI library.The vegan package [62] was used for RDA and variation partitioning.Forward selection was done with packfor library [59].We used library AEM [63] to create AEM variables.

Space/Time Variation in Climate and Limnological Features
Temporal variation in regional climate was investigated for the 13 years of the study.The first two axes of the PCA explained 70.2% of the total variance in climatic conditions (Figure 2a-left panel).Air temperature was opposed to cumulative rain on the first PCA axis (38.8%).This gradient reflected higher temperature and lower rain precipitation during summers of 1991 and 1995 before logging, and in 1998 and 2003 after logging.In contrast, years 1992 and 1993 before logging were colder and rainier.An increase in cumulative snow during winter was associated with the second PCA axis (31.4%) but not related to changes in air temperature and rain precipitation.Year-to-year climatic variation was generally important except between the year 1999 and 2000 which were similar in terms of climate.A time cluster analysis revealed two climatic periods: 1991-1998 and 1999-2003 (Figure 2a-right panel).The difference between these two periods is mainly associated with higher accumulation of snow during winter before 1998.Winter preceding the first clearcut disturbance in 1996 was the snowiest of the survey with total snow precipitation of ≈200 cm while winter of 1998, before the second wave of clearcut logging, was the least snowy with only ≈70 cm of cumulative snow.Summers of 1997 and 1998 were the driest with less than 300 mm of cumulative rain.Mean summer air temperature was at the lowest in 1992 (≈14  S1, for year-to-year variations in climatic variables).

Space/Time Variation in Climate and Limnological Features
Temporal variation in regional climate was investigated for the 13 years of the study.The first two axes of the PCA explained 70.2% of the total variance in climatic conditions (Figure 2a-left panel).Air temperature was opposed to cumulative rain on the first PCA axis (38.8%).This gradient reflected higher temperature and lower rain precipitation during summers of 1991 and 1995 before logging, and in 1998 and 2003 after logging.In contrast, years 1992 and 1993 before logging were colder and rainier.An increase in cumulative snow during winter was associated with the second PCA axis (31.4%) but not related to changes in air temperature and rain precipitation.Year-to-year climatic variation was generally important except between the year 1999 and 2000 which were similar in terms of climate.A time cluster analysis revealed two climatic periods: 1991-1998 and 1999-2003 (Figure 2a-right panel).The difference between these two periods is mainly associated with higher accumulation of snow during winter before 1998.Winter preceding the first clearcut disturbance in 1996 was the snowiest of the survey with total snow precipitation of ≈200 cm while winter of 1998, before the second wave of clearcut logging, was the least snowy with only ≈70 cm of cumulative snow.Summers of 1997 and 1998 were the driest with less than 300 mm of cumulative rain.Mean summer air temperature was at the lowest in 1992 (≈14 °C) and at the highest in 1991, 1995, 1998, 2001, and 2003   Spatial variation in limnological features (morphology, hydrology) among lakes and long-term spatio-temporal variation in water quality were also investigated.Three major variables accounted for 97.7% of spatial variation in limnological features among lakes (Figure 2b-left panel).The first PCA axis (60.6%) indicated an increase in maximum depth and water residence time (WRT), opposing shallow lakes (L80, L38) with short WRT (<2 years) to deeper lakes (L20, L26, L39, L42) with long WRT (>6 years) (Figure 2b-right panel).The second PCA axis (37.1%) was associated to an increase in lake area, opposing large lakes (L80, L20: 55 and 64 ha) to smaller lakes (L26, L38, L39, L42: 18-39 ha).Spatio-temporal variation in water quality was explained (70.4%) by an increasing gradient in alkalinity and associated ions (Ca, Mg, K), nutrients (TP, TON) and DOC concentrations, and algal biomass (CHL) along the first axis (59.2%), in opposition to a decreasing gradient in SO4 (Figure 2c-left panel).The second PCA axis (11.2%) was associated to higher water pH.These gradients in nutrients, chlorophyll and alkalinity were induced by Lake L80 which was mesotrophic, Spatial variation in limnological features (morphology, hydrology) among lakes and long-term spatio-temporal variation in water quality were also investigated.Three major variables accounted for 97.7% of spatial variation in limnological features among lakes (Figure 2b-left panel).The first PCA axis (60.6%) indicated an increase in maximum depth and water residence time (WRT), opposing shallow lakes (L80, L38) with short WRT (<2 years) to deeper lakes (L20, L26, L39, L42) with long WRT (>6 years) (Figure 2b-right panel).The second PCA axis (37.1%) was associated to an increase in lake area, opposing large lakes (L80, L20: 55 and 64 ha) to smaller lakes (L26, L38, L39, L42: 18-39 ha).Spatio-temporal variation in water quality was explained (70.4%) by an increasing gradient in alkalinity and associated ions (Ca, Mg, K), nutrients (TP, TON) and DOC concentrations, and algal biomass (Chlorophyll-a) along the first axis (59.2%), in opposition to a decreasing gradient in SO 4 (Figure 2c-left panel).The second PCA axis (11.2%) was associated to higher water pH.These gradients in nutrients, Chlorophyll and alkalinity were induced by Lake L80 which was mesotrophic, more alkaline, nutrient-rich, and productive than the five other lakes (Figure 2c-right panel).L80 was removed from the PCA analysis to represent more clearly the differences in water quality among interconnected undisturbed (L20, L38) and disturbed (L26, L42, L39) lakes.Once L80 was removed, lakes distributed along the first PCA axis according to an increasing gradient in alkalinity and associated ions (Ca, Mg, K), opposed to a decreasing gradient in nutrients (TP, TON), DOC, and Chl-a (Figure 2d-left panel).Water quality was similar in four lakes (L20, L26, L38, and L42) while L26 was slightly different with higher alkalinity but lower concentrations of nutrients and DOC, and algal biomass (Figure 2d-right panel).
Water 2017, 9, 733 9 of 18 more alkaline, nutrient-rich, and productive than the five other lakes (Figure 2c-right panel).L80 was removed from the PCA analysis to represent more clearly the differences in water quality among interconnected undisturbed (L20, L38) and disturbed (L26, L42, L39) lakes.Once L80 was removed, lakes distributed along the first PCA axis according to an increasing gradient in alkalinity and associated ions (Ca, Mg, K), opposed to a decreasing gradient in nutrients (TP, TON), DOC, and CHL (Figure 2d-left panel).Water quality was similar in four lakes (L20, L26, L38, and L42) while L26 was slightly different with higher alkalinity but lower concentrations of nutrients and DOC, and algal biomass (Figure 2d-right panel).

Space-Time Variation in Zooplankton Community
Over the long-term survey, we observed important variation in zooplankton abundance and composition among lakes (Figure 3).On average, zooplankton was less abundant in L20 and L26 (8 and 11 ind.•L−1 , respectively) than in L39 and L42 (20 and 30 ind.L −1 , respectively), while the highest abundances were recorded in L38 and L80 (92 and 43 ind.•L−1 , respectively) (Supplementary Materials Table S3).Dominance patterns also differed among lakes (Figure 3).In lakes L20, L26, L39, and L42, zooplankton was dominated by calanoid and cyclopoid copepods (75-81%).Rotifers represented a minor proportion of the zooplankton community in L20, L26, and L39 (<18%), and a moderate proportion in L42 (29%).In contrast, calanoid and cyclopoid copepods were less abundant (47-49%) and rotifers more important in L38 and L80 (~45%).Generally, cladocerans accounted for a small percentage of zooplankton abundances (≤7%) in all lakes.Year-to-year variations in zooplankton abundance were important, reaching almost an order of magnitude (Figure 3).The highest abundance of zooplankton occurred before logging in 1991 in all the lakes and also between 1993 and 1995; lower abundances of zooplankton were observed after logging in all lakes except L80 (Figure 3).For lakes, either directly (L26, L42, L39) or indirectly (L20 and L38), impacted by clearcut logging, zooplankton abundance declined during 1996, the first year after logging (Figure 3).Furthermore, all those lakes supported significantly more zooplankton before logging (1991)(1992)(1993)(1994)(1995) than after logging (1996-2003) (Table 2).Post-logging declines in abundance were observed in harvested lakes for several groups of zooplankton: rotifers and cyclopoid copepods in L26, rotifers, calanoid, and cyclopoid copepods in L39, and calanoid copepods in L42.In lakes indirectly impacted by logging, post-logging decline also occurred for the rotifers and the calanoid and cyclopoid copepods in L20 as well as for the cladocerans in L38.In contrast, the unperturbed L80 showed higher abundance of zooplankton in 1996 and no significant change before and after logging disturbance (Figure 3, Table 2).However, year-to-year variation in zooplankton abundance in this unperturbed lake was of the same magnitude (ten-fold) as in the other lakes directly or indirectly affected by clearcut logging (Figure 3).Table 2. Mean abundance (Ind.•L−1 ± SD) of total zooplankton and taxonomic groups in the two groups of lakes (harvested vs undisturbed) before (1991-1995) and after (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) the 1996 initial clearcut logging.Values in bold indicated significant differences in zooplankton abundances before and after logging, both in harvested and undisturbed lakes (Mann-Whitney U test, p < 0.05).

Space-Time Interaction in Harvested and Undisturbed Lakes
Space/time interaction was not detected in harvested lakes (p = 0.527) nor in undisturbed lakes (p = 0.205) (Table 3).This lack of space/time interaction indicated that the temporal pattern of variation in the zooplankton community was common to all sites and that the spatial pattern of variation was stable through time.For the harvested lakes (L26, L39, L42), the main time effect (Year) was significant (R 2 = 0.42; p = 0.003) as was the main spatial effect (Lake) (R 2 = 0.23; p = 0.001).This was also the case for the undisturbed lakes (L20, L38, L80) with a significant main effect for both Year (R 2 = 0.17; p = 0.002) and Lake (R 2 = 0.55; p = 0.001).In the undisturbed lakes, spatial variation in zooplankton abundance among lakes was much more important than temporal variation among years.Conversely, in the harvested lakes, temporal variation among years was more important than the spatial variation among lakes.

Influence of Environment and Logging Disturbance on Zooplankton Space-Time Variation
We tried to detect the influence of environmental drivers on zooplankton abundance while considering temporal (AEM) and spatial (Lake) patterns (Figure 4).For the harvested lakes (L26, L39, and L42), the combined influence of environmental conditions and spatial (Lake) and temporal (AEM) structures explained only 24.3% (R 2 adj %) of zooplankton abundance variation over the 1991-2003 period (Figure 4-top panel).Combined effect of environmental and spatial variables explained 17.8% of the variation in zooplankton abundances.The pure spatial component (Lake) explained the largest proportion (7.5%; p = 0.02) of the total variance of zooplankton abundance among the harvested lakes, whereas the pure temporal component (AEM) explained 4.7% of zooplankton abundance variation among years (p = 0.047).The pure environmental component (climate, limnological features) also explained 4.7% of the variation in zooplankton abundance and was associated to changes in air temperature and DOC concentrations.However, the effects of these environmental factors were not significant taken alone (p = 0.07).Finally, the percentage of clearcut logging in the watershed did not significantly affect zooplankton abundances in the harvested lakes.

Influence of Environment and Logging Disturbance on Zooplankton Space-Time Variation
We tried to detect the influence of environmental drivers on zooplankton abundance while considering temporal (AEM) and spatial (Lake) patterns (Figure 4).For the harvested lakes (L26, L39, and L42), the combined influence of environmental conditions and spatial (Lake) and temporal (AEM) structures explained only 24.3% (R 2 adj %) of zooplankton abundance variation over the 1991-2003 period (Figure 4-top panel).Combined effect of environmental and spatial variables explained 17.8% of the variation in zooplankton abundances.The pure spatial component (Lake) explained the largest proportion (7.5%; p = 0.02) of the total variance of zooplankton abundance among the harvested lakes, whereas the pure temporal component (AEM) explained 4.7% of zooplankton abundance variation among years (p = 0.047).The pure environmental component (climate, limnological features) also explained 4.7% of the variation in zooplankton abundance and was associated to changes in air temperature and DOC concentrations.However, the effects of these environmental factors were not significant taken alone (p = 0.07).Finally, the percentage of clearcut logging in the watershed did not significantly affect zooplankton abundances in the harvested lakes.For the undisturbed lakes, variance partitioning was conducted only for L20 and L38 because RDA was over-influenced by the difference in limnological features of L80 (See Figure 2b,c).The RDA explained 53.3% of the variance in zooplankton abundance from 1991 to 2003 (Figure 4-bottom panel).The combined influence of environmental factors (as expressed by changes in Chl-a, DOC, and rain precipitation, and spatial component (Lake) explained the largest proportion (18.1%) of total variance of zooplankton abundance in these two undisturbed lakes, indirectly impacted.Generally, zooplankton abundance was lower during years with higher cumulative rain, and DOC and Chl-a concentrations.Indeed, DOC seemed to increase during the post-logging period from 1997 to 2003 in the harvested lakes (Supplementary Materials Figure S2a) coinciding with lower abundances of zooplankton (Figure 3).Pure environmental and pure spatial effect (Lake) were both significant (p = 0.046 and 0.043, respectively) while pure time effect (AEM) was not significant (p = 0.075) (Figure 4-bottom panel).The common fraction between environment, temporal and spatial variables was particularly important explaining 16.2% of the variation, indicating more spatio-temporal structure in the environmental factors in the undisturbed lakes than in the harvested lakes (only 0.9%).

Temporal Patterns of Variation in Zooplankton Abundance in Each Lake
Since spatial variation (Lake) was significant for both groups of lakes, we completed the analysis of temporal year-to-year variation in zooplankton abundance for each lake.First, we tried to model the temporal structure as either a linear (regression) or non-linear (AEM) trend (Table 4).Trends were detected for only two of the three harvested lakes.A decreasing linear model explained 23.5% of the long-term variation in zooplankton abundance in L26 and a non-linear model explained 46.8% of the long-term variation in zooplankton abundance in L42 (Table 4).We did not detect significant temporal variation in zooplankton abundance in the most impacted harvested lake (L39) and in the undisturbed lakes (L20, L38, and L80).Secondly, we tried to detect breakpoints in the time series of zooplankton abundance.We found breaks in all directly impacted harvested lakes (L26, L39, and L42) and in one indirectly impacted undisturbed lake (L38) between the years 1995 and 1996, corresponding to the first clearcut logging event (Figure 5).Another break was detected between 1998 and 1999 in two harvested lakes (L39 and L42) after additional clearcut logging in their catchments during summer of 1998.Those two breaks (1995-1996 and 1998-1999) defined a cluster that corresponds to logging disturbances occurring between 1996 and 1998 for L39 and L42.However, breaks also occurred independently of clearcut logging in the harvested lakes between the years 1992-1993 (L42) and 1993-1994 (L26, L39).The undisturbed lake L38 also depicted a break between the years 1992-1993 and the pre-post logging 1995-1996 years.L20 also showed a break between the years 1994-1995.Another important break between 2000 and 2001 was revealed for all harvested lakes (L26, L39, L42), and for an indirectly impacted lake (L20) (Figure 5).Unperturbed L80 did not show breakpoints in the zooplankton abundance time series.

Discussion
Forecasting the state of aquatic systems across the boreal zone is challenging because of the increasing development of natural resources and its interactions with climate changes.Our study of two-decade time series of climate, water quality and zooplankton community of six boreal lakes in Ontario, some of them subjected to direct or indirect impacts of logging, revealed important spatio-temporal variation in abiotic and biotic conditions.There were also significant long-term changes in zooplankton abundance which were not, however, consistently associated with the occurrence and intensity of logging disturbances, thus confirming our hypothesis that variations in limnological features and climate exert a stronger effect on zooplankton community changes than short-term impacts of logging disturbance.There is a considerable potential for interactions between climate change and human perturbations affecting boreal lakes, including overharvesting by clearcut logging of watersheds.Natural variation in regional climate and limnological features may confound the responses of lakes to watershed disturbances at the scale of the ecosystems, even in a replicated experimental design as in our study, because lakes are not identical ecosystems.Furthermore, hydrological connectedness among some of the disturbed and undisturbed lakes, as frequently observed in the boreal Canadian Shield, may have partially compromised our MBACI type design, and confounded the responses of lakes, either directly (L26, L39, L42) or indirectly (L20, L38) affected by clearcut logging disturbances.Indeed, the only lake showing no temporal changes in zooplankton community is headwater L80, which has no hydrologic connection to the other study lakes.Several studies conducted at the watershed scale showed that spatial variation in limnological features confounded the responses of boreal lakes both to climate changes (Schindler et al., 1996) [64] and logging disturbances [1,10].
In our study, year-to-year variation in regional climate was generally important and not in synchrony with pre-and post-logging periods.We detected two climatic periods (before and after 1998) that were characterized by changes in snow accumulation during winter, air temperature, and rain accumulation during summer.On the one hand, the winter preceding clearcut logging of 1996 was the snowiest of the study, which may have increased the inputs of nutrients and DOC from deforested catchments in harvested lakes as observed after 1996.On the other hand, lower snow accumulation during winter after 1998 could have attenuated inputs of nutrients and DOC from catchments to the harvested lakes.Thus pre-and post-logging (wet vs dry) changes in climate were counteracting the potential effects of logging disturbances.Finally, spatio-temporal variation in morphometry, water quality, and trophic conditions among lakes were important in terms of size and depth, drainage basin, water renewal time, and trophic status, inducing differential responses of lakes to long-term climatic changes and short-term logging disturbances.For instance, higher depth, smaller drainage basin, and longer water renewal time in the harvested lakes (L26, L39, and L42) may

Discussion
Forecasting the state of aquatic systems across the boreal zone is challenging because of the increasing development of natural resources and its interactions with climate changes.Our study of two-decade time series of climate, water quality and zooplankton community of six boreal lakes in Ontario, some of them subjected to direct or indirect impacts of logging, revealed important spatio-temporal variation in abiotic and biotic conditions.There were also significant long-term changes in zooplankton abundance which were not, however, consistently associated with the occurrence and intensity of logging disturbances, thus confirming our hypothesis that variations in limnological features and climate exert a stronger effect on zooplankton community changes than short-term impacts of logging disturbance.There is a considerable potential for interactions between climate change and human perturbations affecting boreal lakes, including overharvesting by clearcut logging of watersheds.Natural variation in regional climate and limnological features may confound the responses of lakes to watershed disturbances at the scale of the ecosystems, even in a replicated experimental design as in our study, because lakes are not identical ecosystems.Furthermore, hydrological connectedness among some of the disturbed and undisturbed lakes, as frequently observed in the boreal Canadian Shield, may have partially compromised our MBACI type design, and confounded the responses of lakes, either directly (L26, L39, L42) or indirectly (L20, L38) affected by clearcut logging disturbances.Indeed, the only lake showing no temporal changes in zooplankton community is headwater L80, which has no hydrologic connection to the other study lakes.Several studies conducted at the watershed scale showed that spatial variation in limnological features confounded the responses of boreal lakes both to climate changes (Schindler et al., 1996) [64] and logging disturbances [1,10].
In our study, year-to-year variation in regional climate was generally important and not in synchrony with pre-and post-logging periods.We detected two climatic periods (before and after 1998) that were characterized by changes in snow accumulation during winter, air temperature, and rain accumulation during summer.On the one hand, the winter preceding clearcut logging of 1996 was the snowiest of the study, which may have increased the inputs of nutrients and DOC from deforested catchments in harvested lakes as observed after 1996.On the other hand, lower snow accumulation during winter after 1998 could have attenuated inputs of nutrients and DOC from catchments to the harvested lakes.Thus pre-and post-logging (wet vs dry) changes in climate were counteracting the potential effects of logging disturbances.Finally, spatio-temporal variation in morphometry, water quality, and trophic conditions among lakes were important in terms of size and depth, drainage basin, water renewal time, and trophic status, inducing differential responses of lakes to long-term climatic changes and short-term logging disturbances.For instance, higher depth, smaller drainage basin, and longer water renewal time in the harvested lakes (L26, L39, and L42) may have limited their trophic upsurge response after logging due to a better retention of nutrients and DOC following snowy winters or during rainy summers.Inversely shallow depth, larger drainage ratio and shorter water renewal time in undisturbed lakes increased the inputs of nutrients and organic matter from catchments to undisturbed lakes.
Total abundance of zooplankton was highly variable among lakes and years (1-175 Ind.•L −1 ).Over the long-term survey, mean zooplankton abundance ranged from 8 to 92 Ind.•L −1 among lakes.Zooplankton abundance was higher in undisturbed lakes than in harvested lakes, due to more nutrients and a higher productivity in L38 and L80.Our study confirms that even at small regional scale, more nutrient-enriched and shallow lakes support greater zooplankton abundance than less-enriched and deeper lakes, as observed along a large latitudinal gradient (Pinto-Coelho et al., 2005) [65].Calanoid and cyclopoid copepods accounted for most of the abundance of zooplankton in the deeper lakes with long water residence time, either harvested (L26, L39, L42) or undisturbed (L20), where conditions were more oligotrophic and stable.In opposition, rotifer represented a larger portion of the zooplankton abundance in unstable and oligo-mesotrophic lakes with shallow depth, nutrient-enriched water, and rapid water renewal (L38, L80).Overall, zooplankton dominance and variation patterns in the studied lakes were typical of Precambrian Shield boreal lakes of north-eastern Canada [21,65,66].
In our long-term but small-scale study, both space (among lakes) and time (among years) effects were significant.The lack of space/time interaction indicated that the temporal pattern of variation in the zooplankton community was common to all lakes and that the spatial pattern of variation among lakes was stable through time over the long-term survey.The limited number of lakes and the high among-year variation might have precluded us to demonstrate a space/time interaction for the harvested lake group.Notably, we found differences in the variation partitioning of space and time between harvested and undisturbed lakes.The temporal sources of variation in zooplankton abundance were more important than the spatial sources in harvested lakes, indicating a probable short-term and weak influence of logging disturbance.In contrast, for undisturbed lakes, the space represented a bigger proportion of the variation explained than the temporal sources which indicated a more stable spatial pattern among lakes through time.Our results reflected some of the spatio/temporal variation patterns in zooplankton community observed at larger scales in northeastern America.Spatial variation exceeded temporal variation in abundances of both functional groups and large taxonomic groups of zooplankton in unmanipulated lakes of eastern North America [9].In a large-scale study of temperate lakes of North America, most of the variation in zooplankton abundance was spatial and apportioned primarily among lakes and regions rather than among years [34].However, significant interactions between years and lakes within regions indicated an additional source of variation and demonstrated considerable dependence between space and time variations in zooplankton abundance at large scale.
The composition and abundance of zooplankton communities can help interpret the condition of lakes with respect to natural environmental factors and watershed disturbances [9].Sensitivity of zooplankton indicators for regional lake monitoring has been tested in northeastern United States [7].In general, sensitivity of zooplankton groups such as rotifers, copepod calanoids and cyclopoids, and large and small cladocerans was low across regions.Calanoid copepods (diaptomids) were the only indicators sufficiently sensitive to be used in lake monitoring programs [7].Calanoid copepods were also the most sensitive indicators of logging disturbances in eastern Québec Lakes, declining by 43% in biomass the first year after logging [21].In our study, mean abundances of total zooplankton and some taxonomic groups (rotifers calanoid and cyclopoid copepods) were lower after logging in all three harvested lakes.However, it is difficult to relate these declines to a direct effect of logging disturbances as analogous zooplankton declines were also observed in the undisturbed lakes.Furthermore, we found either a linear or non-linear decreasing trend in zooplankton abundance over the long-term survey in only two out of the three harvested lakes (L26, L42), and not in the most impacted lake (L39).Finally, the results of breakpoints analysis in zooplankton time series did not confirm a direct effect of logging disturbances on zooplankton abundance.Breakpoints occurred after the first logging in 1996 in the three harvested lakes, but also before (1992)(1993) and after (2001) logging.Breakpoints in zooplankton time-series were also observed in two undisturbed lakes (L38, L20) indirectly impacted by logging through their connectivity with harvested lakes.L20 did not show a breakpoint during the clear-cut logging year (1996) but between 1994 and 1995.L38 did show a breakpoint after logging in 1996 and also before in 1992-1993.Finally, L80, which was isolated and mesotrophic, did not show any significant changes in zooplankton abundance during the survey and no breakpoint in the time series.
Variance partitioning of the effects of environmental conditions, space (Lake) and years (AEM) explained higher variance in zooplankton communities of two (L20, L38) of the undisturbed lakes (53%) than in the three harvested lakes (34%).Logging intensity did not have a significant effect on zooplankton community in harvested lakes.Two climatic factors (air temperature and precipitation) and two limnological factors (DOC, Chl-a) were retained in the models.However, their effects were weak or not significant when temporal and spatial effects were taken out.Lower temperature during rainier summers may increase the inputs of nutrients and DOC, enhancing the potential impact of watershed logging [12].In boreal lakes of eastern Quebec, an increase of DOC, by limiting light penetration, dampened harvesting effects due to higher inputs of nutrients, and the trophic upsurge of both phytoplankton and zooplankton after logging [20,21].In our study, DOC seemed to increase during the post-logging period from 1997 to 2003 in the harvested lakes coinciding with lower abundances of zooplankton both in the harvested and undisturbed lakes.Lake models suggested that harvesting operations must be extensive and last a long period to induce enough change in DOC and light penetration to significantly influence planktonic food web [33].Previous studies in the Boreal ecozone of Canada suggested that significant changes in zooplankton community would be detected more easily in lakes with a large drainage ratio (watershed/lake volume ratio >4) and important clearcut logging (>40% of the watershed) [12,67].Finally, we were not expecting strong cascading effects of the littoral fish, mainly cyprinids, because logging impacts in cyprinids of the harvested lakes were small [25].
In summary, our long-term study reports few, transient, or no measurable effects of logging disturbances on zooplankton community beyond natural variation induced by climate changes and limnological variation at the regional scale of the Coldwater watersheds as previously reported in short-term studies of the impacts of forest harvesting on zooplankton communities of boreal lakes in the Canadian Shield ecoregion [17,[21][22][23].Zooplankton community structure showed high natural variability among lakes and years that seemed to be important in harvested as well as in undisturbed lakes.Therefore, natural variability in climate and limnological features seemed to be more important in determining the ecological patterns within boreal lakes of eastern Canada than short-term impacts of logging.Considering the climatic change projection in the boreal shield, where temperature is projected to increase by 8.0 • C during winter and an overall increase in precipitation of 11% distributed more in winter and spring [68], these findings bear important meanings since longer growing seasons and frost-free periods will drastically modify the limnological features and aquatic biodiversity of north-temperate lakes [6].We expect less zooplankton abundance in boreal lakes with warmer and wetter climate, particularly for sensitive groups as the Copepoda Calanoida which would decline when higher rain accumulation and runoff will increase DOC and nutrients and decrease water transparency [21].

Acknowledgments:
The Ontario Ministry of Natural Resources funded this research.We acknowledge M. Bozek for supervising the Coldwater Lakes field program from 1991 to 1993, and the CNFER technical staff, particularly Pat Furlong, Mike Friday, and Blake Laporte for sampling and fieldwork, and Mike Carneiro for data management.Special thanks to Emma Mangas who helped to separate and organise the zooplankton samples before taxonomic analysis.The postdoctoral fellow (D.Lévesque) was supported by Discovery grants from the National Science and Engineering Research Council (NSERC) to B.P.A.

2. 1 .
Study Area, Multiple Before/After-Control-Impact (MBACI) Experimental Logging Design and Limnological Features The study area was in the Coldwater Lakes Experimental Watersheds in the Boreal-Great Lakes transition forest on the Canadian Shield (49 • 07 16 -40 • 08 05 N; 92 • 12 22 -92 • 13 34 W), approximately 200 km northwest of Thunder Bay, Ontario, Canada (Figure 1).The climate was characteristic of the north-temperate boreal ecoregion, with mean air temperatures of

Figure 2 .
Figure 2. Biplot of the principal component analysis (PCA) based on correlations: (a) among climate variables measured from 1991 to 2003 (n = 13), (b) among morphometric variables of the six monitored lakes (n = 6), (c) among physical and chemical variables of the six monitored lakes from 1991 to 1999 (n = 76), and (d) among physical and chemical variables in the monitored lakes (except L80) from 1991 to 1999 (n = 63).L80 was deleted to represent more clearly the differences in water quality among interconnected undisturbed and disturbed lakes.For all left panels, arrows represent environmental variables within the correlation circle.Sites are represented on the right panel.Additionally, sites for climatic variables (panel a) are represented following a time constrained cluster analysis based on a multivariate regression tree analysis of the three climatic variables where two groups were found, as indicated by a solid line (1991-1998) and a dotted line (1999-2003).Confidence ellipses based on Lake are shown for physical and chemical variables in panel c.

Figure 2 .
Figure 2. Biplot of the principal component analysis (PCA) based on correlations: (a) among climate variables measured from 1991 to 2003 (n = 13), (b) among morphometric variables of the six monitored lakes (n = 6), (c) among physical and chemical variables of the six monitored lakes from 1991 to 1999 (n = 76), and (d) among physical and chemical variables in the monitored lakes (except L80) from 1991 to 1999 (n = 63).L80 was deleted to represent more clearly the differences in water quality among interconnected undisturbed and disturbed lakes.For all left panels, arrows represent environmental variables within the correlation circle.Sites are represented on the right panel.Additionally, sites for climatic variables (panel a) are represented following a time constrained cluster analysis based on a multivariate regression tree analysis of the three climatic variables where two groups were found, as indicated by a solid line (1991-1998) and a dotted line (1999-2003).Confidence ellipses based on Lake are shown for physical and chemical variables in panel c.

Figure 3 .
Figure 3. Stacked barplot of the abundance of total zooplankton (Ind.•L−1 ) and the relative abundances of taxonomic groups in the six monitored lakes from 1991 to 2003 (black: Cyclopoida; pale grey: Calanoida; dark grey: Rotifera; white: Cladocera).Arrows represent connectivity between lakes.Percentages of the watershed deforested by clearcut logging in 1996 and 1998 are also shown for harvested lakes.

Figure 3 .
Figure 3. Stacked barplot of the abundance of total zooplankton (Ind.•L−1 ) and the relative abundances of taxonomic groups in the six monitored lakes from 1991 to 2003 (black: Cyclopoida; pale grey: Calanoida; dark grey: Rotifera; white: Cladocera).Arrows represent connectivity between lakes.Percentages of the watershed deforested by clearcut logging in 1996 and 1998 are also shown for harvested lakes.

Figure 4 .
Figure 4. Percentage (R 2 adj %) of variation in zooplankton abundance explained by ENVIRONMENTAL (df = 2 and df = 3, respectively); SPATIAL (df = 2 and df = 1, respectively), and TEMPORAL (df = 1 and df = 2, respectively) variables for harvested (top panel: n =39) and undisturbed (bottom panel: n = 24) lakes.A forward selection was applied on environmental and temporal variables.Values < 0 are not shown.Based on PCA results, L80 was excluded from this analysis.

Figure 4 .
Figure 4. Percentage (R 2 adj %) of variation in zooplankton abundance explained by ENVIRONMENTAL (df = 2 and df = 3, respectively); SPATIAL (df = 2 and df = 1, respectively), and TEMPORAL (df = 1 and df = 2, respectively) variables for harvested (top panel: n =39) and undisturbed (bottom panel: n = 24) lakes.A forward selection was applied on environmental and temporal variables.Values < 0 are not shown.Based on PCA results, L80 was excluded from this analysis.

Figure 5 .
Figure 5.Time clusters on the variation in zooplankton abundance in each monitored lake based on multiple regression tree analysis.Full arrow shows the first year of clearcut logging in 1996 and the dotted arrow represents the second partial logging in 1998 for L39 and L42.

Figure 5 .
Figure 5.Time clusters on the variation in zooplankton abundance in each monitored lake based on multiple regression tree analysis.Full arrow shows the first year of clearcut logging in 1996 and the dotted arrow represents the second partial logging in 1998 for L39 and L42.

Table 3 .
Effect of time, space and their interaction on zooplankton abundance in the two groups of lakes (harvested vs undisturbed).Values are significant at p ≤ 0.05.All tests were performed with 999 permutations.* significant effect.

Table 4 .
Modelling of the temporal patterns for the zooplankton groups (Cyclopoida, Calanoida, Rotifera, and Cladocera) for each harvested and undisturbed lake from 1991 to 2003.Bold values indicated significant long-term temporal pattern of variation.AEM: Assymetric Eigenvector Map.