SITE DESCRIPTION AND CORE COLLECTION

During the Last Glacial Maximum, Hudson Bay was the growth center for the Laurentide Ice Sheet and thus the site of the thickest ice cover (Sella et al., 2007), depressing the Earth's crust well below present-day levels (Stewart and Barber, 2010). The low-lying coastal plains of the HBL are a legacy of this last ice age and are currently ∼120 m a.s.l. at their highest point. Adjustment of the Earth's crust to the large change in surface load as the ice sheet receded is why the HBL experiences one of the highest rates of isostatic rebound in the world, between 1.2 and 1.3 m vertical increase per century (Webber et al., 1970; Stewart and Barber, 2010). The emergent marine sediment is largely composed of fine silts and clays with low permeability (Martini, 1981). The low permeability of underlying glacial deposits and marine silty clays, as well as the dominance of permafrost and low topographic relief, makes the HBL an ideal landscape for waterlogging and wetland development (Martini, 2006; Kuhry and Turunen, 2006; Sannel and Kuhry, 2011). Most of the lakes in the Sutton River region of the HBL are likely of glacial origin (unlike thermokarst lakes that are more common in the Churchill area), and it is possible that some of the lakes near the coast are the product of recently emerged shorelines from the high rates of isostatic rebound in this region.

The HBL is the most southern subarctic region in the world due to the cooling effect of the world's largest polar inland sea, Hudson Bay (Fig. 1). The Arctic influence of this cold sea on the surrounding landscape is apparent in numerous ways. For example, the steep southward dip of the northern limit of the treeline and the sharp transition in permafrost development to the west of HB closely tracks the mean July position of the Arctic Front (Bryson, 1966), which is largely a response to cold air masses generated over HB (Rouse, 1991). Despite its southern latitude (e.g., similar latitude to Glasgow, Scotland), the Sutton River region of the HBL has a distinctly “northern” character with a polar (or nordicity) value (French and Slaymaker, 2012) exceeding that of Yellowknife, Northwest Territories, ∼8° to the north. The climate of the HBL region is not only influenced by Arctic air masses, but also by the circulation of HB sea water, as cold Arctic waters enter HB in the northwest corner of Foxe Basin and Hudson Strait, generating currents that cool the western coast of HB (Martini 2006). The Arctic influence of the cold waters of HB penetrates far inland, resulting in a cold region dominated by permafrost (Rouse, 1991; McKendry and Roulet, 1994). Based on 1971–2000 averages for this part of the HBL, the mean annual temperature is -4.5 °C, the mean annual maximum temperature is 17.3 °C, the mean annual minimum temperature is -29.1 °C, and the average annual precipitation is 518.7 mm (McKenney et al., 2010).

Sediment cores were collected from 13 lakes in late July to early August of 2009 or 2010 (55°6′48.6″N, 84°2′30.12″W; 53°54′29.7″N, 84°8′42.72″W) (Fig. 1). These lakes were sampled to cover a geographical area of northern Ontario from the HB coastline to ∼120 km inland, near the Sutton River. Seven of the lakes are named (Hawley, North Raft, Spruce, Aquatuk, Kinushseo, Opinnagau East, and Raft Lakes), and the remaining six lakes were given unofficial names (Billbear, Stuart, Julison, Wolfgang, Sam, and Cassie Lakes) (Fig. 1). Sediment cores were retrieved from approximately the deepest part of each lake (usually at the center) by lowering a Glew (1989) gravity corer (inside core tube diameter of 7.5 cm) off the pontoon of a float plane. The surface sediments (uppermost 0.5 cm interval) and the bottom sediments (bottom-most undisturbed 0.5 cm interval: mean = 21.75 cm) were extruded in the field using a Glew (1988) vertical extruder and placed into Whirlpak^® bags.

For each lake, limnological data were collected between 2009 and 2011, with some lakes monitored each field season, some for two seasons, and some for only one year (Table 1). For water chemistry analyses, a composite water sample of the lake water column was collected at approximately the deepest part of the lake from the Secchi depth to the surface. The Secchi depth of some lakes extended to the lake bottom, and these were sampled from 0.5 m above the lake bottom to the surface. All chemical analyses were performed at the Ontario Ministry of the Environment's (MOE) Dorset Environmental Science Centre, using standard MOE protocols (Ontario Ministry of the Environment, 1983). In addition, profiles of water temperature and oxygen levels were recorded at each site. Paterson et al. (this issue) provides a detailed account of the limnological data and collection procedures.

The transect of 13 study lakes crosses two zones of permafrost development. Five of the study lakes (Billbear, Stuart, Julison, Sam, and Wolfgang Lakes) lie within the zone of continuous permafrost, a thin, sparsely vegetated band from the HB shoreline to approximately 60 km inland (Fig. 1). The remaining eight lakes are located in the zone of discontinuous permafrost, where needle-leaf trees occur. Four of the 13 study lakes are relatively deep (maximum depths ∼11 to 35 m) and located at a relatively high elevation (120 m above the surrounding lowlands) near the Sutton Ridges, an outcrop of Precambrian Shield bedrock (Gunn and Snucins, 2010) within an otherwise limestone-dominated region (Fig. 1; Table 1). These four deep lakes are somewhat of an anomaly in this region, whereas the remaining nine study lakes (0.7 to 2.7 m deep) are more typical of the shallow ponds that characterize the HBL (Fig. 1; Table 1). For more details regarding lake selection, see Paterson et al. (this issue).

FIGURE 1.

Location of the 13 study lakes in the Hudson Bay Lowlands, Ontario, Canada.

The study lakes range in nutrient concentrations from oligotrophic to mesotrophic (total phosphorus (TP) range = 4.3 to 28.8 µg/L, average = 11.7 µ/L) with some of the highest values recorded in the shallowest lakes located close to the HB coastline (e.g., Billbear, Julison, Stuart, and Wolfgang Lakes) (Table 1). All lakes were slightly alkaline (pH range = 7.6 to 8.1, average = 7.9) and are strongly phosphorus limited (Table 1). Water column profiles of temperature and oxygen measured indicate that all lakes were well oxygenated throughout their water columns and that the deeper lakes at the time of sampling were weakly to moderately stratified, although strong thermal stratification was observed at Hawley Lake in 2001 (Gunn and Snucins, 2010). A more detailed description of the water chemistry of our sites is provided by Paterson et al. (this issue).

THE TOP-BOTTOM PALEOLIMNOLOGICAL APPROACH

To examine whether modern-day algal (diatoms) and invertebrate (chironomids and cladocerans) sedimentary assemblages differ from past assemblages, we use a paleolimnological approach that compares two discrete time intervals across a spatial gradient (the “top-bottom” or “before and after” approach) (Smol, 2008). These types of regional assessments can efficiently address broad-scale questions about environmental change. The bottom samples provide an integrated sample of the biological communities present in the pre-impact environment (e.g., baseline or natural conditions). In contrast, the top surface sediment samples provide an integrated sample of the biological communities present in modern-day aquatic habitats (post-impact). Given the shift in the climate system of the HB region from a relatively cool and stable state to one of intense warming in the mid-1990s (Gagnon and Gough, 2005; Hochheim and Barber, 2010; Bhiry et al., 2011), and the reported stability among biological assemblages of lakes in the HB region over at least the past ∼200 years (Laing et al., 2002; Pienitz et al., 2004; Smol et al., 2005) and likely the past ∼1500 years (e.g., Ponader et al., 2002; Fallu et al., 2005), our objective was to determine whether recent biological assemblages (i.e., post-1990s warming) are different from those assemblages deposited prior to the onset of 1990s warming. Although we cannot provide an absolute age for each bottom sample, we are confident that bottom samples (average bottom interval = 21.4 cm, minimum bottom interval = 9.5 cm, maximum bottom interval = 30 cm; Appendix Table Al) were deposited prior to the onset of regional warming.

TABLE 1

A selection of measured limnological variables for the 13 study lakes. Data presented are 3-year means (2009 to 2011).

Comparison of biological assemblages preserved in the top and bottom samples allows for a rapid assessment (relative to detailed core analyses) of lakes across a spatial gradient, providing important insights into regional changes. This is a particularly effective approach for exploratory research in remote locations where little direct environmental monitoring has been performed and baseline data are lacking (Smol, 2008). Although top-bottom paleolimnological analyses have advantages over detailed sedimentary analyses (i.e., relatively cost- and time-effective), it is important to understand the limitations and assumptions of this approach (Smol, 2008). As only two discrete sedimentary intervals are examined from each lake, the aim is to assess whether the modern samples are different from pre-impact samples, as the specific timing of any identified changes cannot be determined (Smol, 2008). It is also important to note that a given sedimentary slice may represent a different number of years of accumulation across lakes (e.g., the top 0.5 cm in one lake may represent ∼3 years of accumulation, whereas the bottom 0.5 cm may represent ∼8 years of accumulation). Furthermore, the bottom sedimentary intervals across the study lakes may represent different points in time, and it is therefore important to determine that the bottom sample is truly representative of pre-disturbance conditions. In this study, 6 of the 13 lakes have detailed cores available that have been radioisotopically dated (see below) providing us with an approximation of the average age represented by the bottom samples and an estimate of the approximate number of years represented by top and bottom samples. Given the available evidence that the climate of the Hudson Bay region has remained relatively unchanged over the past ∼200 years (e.g., Laing et al., 2002; Pienitz et al., 2004; Smol et al., 2005), and likely longer, we are confident that our bottom samples (although undoubtedly different in age) represent a relatively stable climatic environment prior to the onset of substantial warming over the past few decades and that biological community changes are a response to recent warming trends.

RADIOISOTOPIC ANALYSES

To ensure that the bottom sediment samples were deposited prior to the onset of warming in the HBL, ²¹⁰Pb analysis using gamma spectrometry was performed following Appleby (2001). Chronologies from 6 of the 13 lakes were determined via radiometric dating (Hawley, North Raft, Spruce, and Aquatuk Lakes—Rühland et al. (2013); Sam and Wolfgang Lakes— unpublished data) using the constant rate of supply (CRS) model (Appleby, 2001). ²¹⁰Pb activities from these four lakes together with ²¹⁰Pb activities from the top and bottom sedimentary intervals of 8 of the 9 remaining lakes (Sam, Opinnagau East, Billbear, Stuart, Wolfgang, Julison, Cassie, and Kinushseo Lakes) were analyzed (there was insufficient sediment from Raft Lake for this analysis). A comparison of the decay of the ²¹⁰Pb activities in the top and bottom intervals to ²¹⁴Bi (a proxy for background or supported ²¹⁰Pb activity levels) allowed for a determination of whether the bottom samples represent background conditions.

SEDIMENT SAMPLE ANALYSES: LOSS ON IGNITION, DIATOMS, CHIRONOMIDS, CLADOCERANS

NUMERICAL ANALYSES

Changes in the biological communities between top and bottom sediment samples (i.e., differences between present-day and pre-disturbance conditions) from each lake were summarized through principal components analysis (PCA), using the default options available in the program CANOWIN, version 4.5 (ter Braak and Šmilauer, 2002) as well as the vegan package (Oksanen et al., 2010) for the R software environment (R Development Core Team, 2010). All percent relative abundance species data were square-root transformed prior to analyses to equalize the variance among taxa. This unconstrained ordination procedure was used to examine how the lake samples varied with respect to species composition (i.e., samples plotted in species space) and is referred to herein as PCA_spp . For each biological proxy, pre-disturbance (i.e., bottom) samples were run passively in the ordinations, and all modern (i.e., top) samples were run actively. The results of the PCA ordinations were used in two different ways. First, sample scores for PCA axis 1 were plotted against PCA_spp axis 2 for each biological proxy. Plotting the modern and pre-disturbance assemblages in the same ordination space enabled a comparison of the relative changes in biological composition through time (modern and pre-disturbance), highlighting the similarities and differences among sites in terms of both species composition and magnitude of change (i.e., trajectory through time). Second, a PCA was used to examine how the lake samples varied with respect to measured environmental variables (i.e., samples plotted in environmental space) and is referred to as PCA_env . Environmental variables were transformed prior to analyses (see Paterson et al., this issue), and a Pearson correlation matrix was used to determine correlations among the measured variables. This allowed us to better understand some of the underlying patterns in our species distribu- tions and make comparisons among the similarities and differences of the study lakes by comparine PCA_spp and PCA_envp plots.

We examined the possibility of performing a redundancy analysis (RDA) using the program CANOWIN, version 4.5 (ter Braak and Šmilauer, 2002), to examine relationships between the distribution of modern biological indicators and the measured environmental variables from the 13 lakes. An RDA with forward selection only identified one variable (lake depth) as explaining a significant amount of the variation in both the chironomid and cladoceran assemblages and was therefore determined not to be an informative exercise. This same procedure resulted in four significant variables (lake depth, conductivity, TP, and color) for the diatom assemblages. However, given the lack of dimensionality for the chironomid and cladoceran indicator groups, RDA was not pursued further in this study.

To test whether there were significant differences in species' assemblages through time and across all study lakes, an analysis of similarity (ANOSIM) was performed using the vegan package (Oksanen et al., 2010) for the R software environment (R Development Core Team, 2010). To compare the magnitude of change among the three indicator groups, species relative abundance data were transformed into a dissimilarity matrix, and the Bray-Curtis dissimilarity coefficient (DC) (Clarke et al., 2006) was calculated between each top and bottom interval for each lake. For each lake, diatom DCs were compared to corresponding cladoceran and chironomid DCs, and cladoceran DCs were compared to corresponding chironomid DCs (e.g., Sweetman et al., 2008) to determine which of the indicator groups has experienced the greatest amount of change through time.

To determine whether there was a global (i.e., across all lakes) trend at the species level, a comparison was made of the percent relative abundances between top and bottom samples for select species and genera from each indicator group. In addition, differences in species diversity between modern and pre-disturbance samples were calculated using Hill's N2, where N2 represents the number of very abundant taxa in a given sample (Hill, 1973). Prior to calculating species diversity, all raw species data were rarefied to a common base sum to offset the differences in the numbers of individuals counted for each sample (Birks and Line, 1992; Birks and Birks, 2008). All species were included in this test (i.e., no cutoff criterion was applied). To test whether the species abundances and species diversity (N2) in the modern samples were significantly different from species abundances and diversity in the pre-disturbance samples, we used a nonparametric Wilcoxon Signed-Rank test, using the vegan package (Oksanen et al., 2010) for the R software environment (R Development Core Team, 2010). This measure should be viewed with some caution as sediment intervals at the bottom of the core generally represent a longer time period than modern intervals, potentially resulting in higher diversity values in the bottom samples (Smol, 1981).

RADIOISOTOPIC TRENDS

Radioactive decay of ²¹⁰ Pb was at background levels (i.e., approached the ²¹⁴Bi curve) in all eight of the bottom samples analyzed with gamma spectrometry as well as the six lakes where detailed dating analysis was undertaken (Hawley, North Raft, Spruce, and Aquatuk Lakes—data not shown for Wolfgang and Sam Lakes) (Appendix Fig. A1). Comparisons of radioactive activity between top and bottom samples (Fig. A1) suggest that our top intervals represent the last few years of sediment accumulation, and the majority of our bottom intervals represent background (i.e., ∼150 years). Considering that this region has only experienced substantial increases in temperature in the past few decades, all of our bottom lake sediment samples undoubtedly represent pre-warming conditions in the HBL.

GENERAL TRENDS IN PALEOINDICATORS

For each indicator group (diatom, chironomid, and cladoceran), it was often challenging to attain adequate counts for useful environmental interpretations in each lake. As a result, only subsets of the original 13 top and bottom samples (i.e., 26 sedimentary intervals) were available for analyses, and lake samples that lacked subfossil remains differed among the indicator groups (Appendix Table Al). The percent organic matter content (%OMC) measured through LOI at 550 °C varied considerably across lakes (Table A1). However, the modern intervals in all lakes had higher %OMC than the bottom intervals.

In the four shallowest lakes (Billbear, Julison, Cassie, and Stuart Lakes—all <1.0 m depth), which were located close (i.e., between ∼11 and ∼70 km) to the HB coastline, diatom valves were too scarce for a meaningful count (often less than 10 valves in an entire microscope slide) in the pre-disturbance (i.e., bottom) samples. In these lakes, other autochthonous siliceous indicators (e.g., chrysophyte cysts and scales) were equally scarce, whereas siliciclastic materials were extremely high. Due to the almost complete lack of diatom valves in these four samples, sodium polytungstate (SPT) treatment (Tapia and Harwood, 2002) was not deemed useful. The bottom samples of Hawley, North Raft, Spruce, and Aquatuk Lakes were treated with SPT, resulting in concentrated diatom samples that were adequately low in siliciclastic materials to allow counts of satisfactory quality (i.e., between 200 to 300 valves) for the bottom intervals (Table A1). In all samples, diatom valves and the thinly silicified spines of certain chrysophyte cysts that were encountered (even when their numbers were very low) were remarkably well preserved, suggesting that dissolution was not the reason for the scarcity of siliceous indicators in many of the “bottom” samples. The resulting lake set available for the diatom analyses consisted of 13 modern samples and 9 pre-industrial samples.

The ability to attain adequate counts for chironomid analysis varied across lakes. Although considerable effort was dedicated to attain the standard minimum of 50 head capsules for each sample, this was not possible for several lakes due to the scarcity of chironomids in the sediment and/or inadequate amounts of sediment available (Table A1). Nevertheless, the average number of head capsules counted across all samples was 68 (min = 32, max =112). Achieving adequate counts for the bottom sediment samples of Julison (38), the top of Sam (36), and both the top (40) and bottom (32) of Stuart Lakes was not possible. Achieving higher taxonomic resolution was likewise difficult in these samples due to the poor preservation of diagnostic characteristics (e.g., premandibles and mandibles). Furthermore, within the pre-disturbance samples of the shallowest lakes near to the HB coast (Billbear, Julison, Stuart, and Cassie Lakes), achieving an adequate number of chironomid subfossils was challenging due to the apparent poor preservation of the chitinous head capsules.

The density of cladoceran remains was low in all samples, and they were virtually absent in the pre-disturbance samples from Cassie, Julison, Spruce, and Stuart Lakes, as well as in the modern sediments of the latter. Due to this low density of cladoceran remains, the recommended number of individual clad-ocerans for a representative estimation of the sedimentary cladoceran assemblage (i.e., 70; Kurek et al., 2010) was not achieved for all samples (Table A1). However, when overall taxonomic richness is low, smaller counts of individuals can be sufficient to characterize cladoceran assemblages (e.g., Nevalainen, 2010). Therefore, we are confident that the lower than ideal counts are representative of the cladoceran assemblage composition. The resulting lake set available for the cladoceran analyses consisted of 12 modern samples and 9 pre-disturbance samples. Due to difficulties associated with the degree of fragmentation and debris in the samples, we were unable to distinguish the specific members of Bosmina spp.

CHANGES IN BIOLOGICAL ASSEMBLAGE COMPOSITION

For both algal and invertebrate groups, no significant difference across all 13 lakes was found in a comparison of the present-day assemblages with the historical assemblages (one-way ANOSIM, diatoms: R = 0.02, p = 0.3; chironomids: R = -0.05, p = 0.87; cladocerans: R = -0.13, p = 0.98; n = 999 permutations). Although ANOSIM did not identify a significant change for any indicator group (i.e., at the assemblage level) across lakes, there were nonetheless significant changes in some of the lakes, as well as significant and ecologically interesting changes at the species level between top and bottom samples (Fig. 2). For example, the modern intervals of all lakes had significantly higher percent relative abundances of planktonic diatoms (Cyclotella taxa, p = 0.01; Fragilaria tenera, p = 0.01) and significantly lower abundances of small benthic fragilarioid taxa (Fragilaria construens complex, p = 0.02) (Fig. 2, part A). The only chironomid species showing significantly higher abundances in the recent sediments was the warm-water, littoral Cladotanytarsus mancus-type (p = 0.012). Corynocera ambigua also showed slight increases in modern sediments (p = 0.06), whereas the grouped species of littoral Cricotopus had significantly higher relative abundances in the pre-disturbance samples in almost all lakes (p = 0.007) (Fig. 2, part B). Species from the littoral genus Psectrocladius decreased in several shallow lakes in the modern samples, but the difference was not significant across all lakes (p = 0.15; Fig. 2, part B). Although there was similarly little change since the 19th century within the cladoceran assemblages, two cladoceran taxa had significant differences in relative abundance between the top and bottom samples (Fig. 2, part C): the planktonic taxa Bosmina spp. (higher abundances in the present-day sediment samples; p = 0.03), and the littoral taxon Alona guttata (lower abundances in present-day samples; p = 0.04). Although Chydorus brevilabris was dominant in the study lakes, differences between top and bottom relative abundances were not significant (p = 0.36) (Fig. 2, part C).

There was a clear, significant increase in diatom species diversity (N2) (p = 0.01) in all but one lake (Fig. 2, part A), but no trend was observed for cladoceran diversity (p = 0.91) (Fig. 2, part B), or for chironomid diversity (p = 0.74) (Fig. 2, part C). A comparison of the Bray-Curtis dissimilarity coefficients (DCs) indicated different magnitudes of change among the three indicator groups (Fig. 3, parts A-C), with the highest change recorded in the diatom assemblages (mean DC = 0.40, n = 9) closely followed by the chironomids (mean DC = 0.36, n = 13) and with the least change recorded in the cladocerans (mean DC = 0.26, n = 9). Wilcoxon Signed-Rank tests deemed that differences in DCs between diatoms and chironomids were not significant (p = 0.82, n = 9), with five lakes showing greater change in the diatoms and four lakes showing higher change in the chironomids (Fig. 3, part A). Differences in DCs were significantly higher in diatoms than in cladocerans in all but one lake (p = 0.05, n = 9) (Fig. 3, part B) and were significantly higher in chironomids than in cladocerans in all lakes (p = 0.03, n = 9) (Fig. 3, part C).

SPECIES AND ENVIRONMENTAL ORDINATIONS

Examination of PCA biplots of the top (modern) sediment assemblages with the bottom (pre-disturbance) assemblages run passively in the same ordination space showed that PCA axis 1 explained 51% of the diatom variance, 39% of the cladoceran variance, and 21% of the chironomid variance in the species data (Fig. 4, parts A-C). PCA axis 2 explained 11%, 17%, and 17% in species variance for diatoms, cladoceran, and chironomid, respectively. PCA trajectories through time (i.e., comparisons between top and bottom assemblages for each lake) displayed highly variable magnitudes of change within and among indicator groups. PCA_env plots of the samples and the environmental variables (Fig. 4, part D) showed that axis 1 represented a gradient of lake depth and nutrients, whereas PCA_env axis 2 represented a mixed gradient of conductivity (and related variables), dissolved organic carbon (DOC), and color. These environmental gradients in the PCA_env are consistent with the distribution of species in the PCA_env for all indicator groups, with PCA_spp axis 1 clearly representing a gradient of lake depth with the deeper lakes (∼ 11.0 to 35 m) plotting in the opposite quadrant to the shallow lakes (0.7 to 2.7 m) in the ordinations. Surface area did not show a clear separation among sites in the ordination, with two of the shallowest and smallest lakes (e.g., Sam, 32 ha, and Cassie, 70 ha) plotting close to the larger and deeper lakes (e.g., Aquatuk, 828 ha, and North Raft, 668 ha) (Fig. 4, part D). Furthermore, the two largest lakes include the deepest (Hawley Lake, 35 m, 1235 ha) and one of the shallowest (Opinnagau East Lake, 1.55 m, 2723 ha) in the data set, and they plotted on opposite sides of the ordination along PCA axis 2 (Fig. 4, part D).

Diatom assemblages examined in the 13 HBL lakes varied in both species composition and the magnitude of change from pre-disturbance to modern times, with a clear separation of sites along a lake depth gradient on PCA axis 1 (Fig. 4, part A). The shallower lakes on the left side of the ordination contained a variety of benthic and periphytic diatoms, particularly from the genera Achnanthes, Cymbella, and Brachysira, and from the larger taxa of Navicula (e.g., N. cryptotenella, N. cryptocephala, and N. radiosa), and differences between top and bottom assemblages largely consisted of changes among these groups. In contrast with the diverse benthic assemblages of the shallower lakes, diatom assemblages in the pre-disturbance samples of the deeper lakes were dominated almost exclusively by small epilithic fragilarioid taxa including Fragilaria brevistriata varieties, F. construens varieties, F. leptostauron varieties, and F. pinnata varieties (Fig. 2, part A). Planktonic diatom taxa (e.g., Fragilaria tenera, Asterionella formosa, Cyclotella stelligera, and C. michiganiana) were almost exclusively observed in the modern assemblages of most of the study lakes, with generally higher abundances in the deeper lakes (Fig. 2, part A).

Although there was a clear separation of shallow sites from deep sites in the modern diatom assemblages along PCA axis 1, there are two notable exceptions. First, sample scores for the modern diatom assemblages from Billbear, one of the shallowest (0.8 m) and smallest (77 ha) lakes in the data set, had the greatest separation from the rest of the 13 lakes and plotted in the upper quadrants along PCA_spp axis 2 (Fig. 4, part A). Billbear Lake also had the highest total phosphorus concentrations (TP = 28.8 µ/L), well above the mean (TP = 12.1 µg/L: Table 1) and plotted on the far left of the PCA_env (Fig. 4, part D). Second, both modern and pre-disturbance sample scores for the largest site by surface area, Opinnagau East Lake (1.70 m deep, 2723 ha) plotted closer to PCA_spp sample scores for deeper lakes rather than shallow lakes (Fig. 4, part A), indicating that diatom assemblage composition was distinct from the generally similar assemblages found in the other shallow study lakes. Much like the assemblages found in the deeper lakes, Opinnagau East Lake consisted of high relative abundances of small, benthic Fragilaria pinnata taxa and a notable lack of the diverse assemblage of benthic diatoms such as Brachysira and Cymbella species and large Navicula cryptotenella and N. cryptocephala that were common to the shallow sites of this study (plotting on the left within the PCA_spp plot).

FIGURE 2.

Scatter plots showing the relationship between species diversity (Hill's N2) and the percent relative abundances of a selection of ecologically significant taxa and taxonomic groups recorded in the bottom and top sedimentary intervals of the 13 Hudson Bay Lowlands (HBL) study sites for (A) diatoms, including planktonic Cyclotella taxa (Dìscostella pseudostellìgera, Cyclotella michiganiana), the benthic Fragilaria construens complex (F. construens, F. construens var. venter), and planktonic Fragilaria tenera; (B) chironomids, including Cladotanytarsus mancus-type, Cricotopus taxa (C. intersectus-type, C. cylindraceus-type, C. type C, C. type P), and Psectrocladius taxa (P. calcaratus-type, P. septentrionalis-type, P. sordidellus-type, P. barbatipes-type); and (C) cladocerans. Samples that plot above the 1:1 line indicate that higher relative abundances of that taxon occurred in the top (modern) samples than in the bottom (pre-disturbance) samples. Significant differences between all tops and all bottoms were determined using a Wilcoxon-signed-rank test

Wolfgang Lake (1.2 m deep) displayed the largest diatom change through time, with the pre-disturbance assemblage having high relative abundances of benthic fragiliarioid taxa and small Navicula minima, similar to the bottom sediments of the deeper study lakes (Fig. 4, part A). Meanwhile, the modern diatom assemblages from Wolfgang Lake consisted of taxa common to the other shallow lakes of this study (e.g., high relative abundances of more diverse and complex benthic diatoms such as Brachysira vitrea, Cymbella taxa, and Achnanthes taxa) (Fig. 4, part A). Modest trajectories of species changes through time were observed for Raft, North Raft, Aquatuk, and Opinnagau East Lakes, with the least amount of change observed in Kinushseo, Sam, Hawley, and Spruce Lakes. Due to the paucity of diatoms in the pre-disturbance sediments of the shallow coastal sites, estimates of floristic changes in Billbear, Julison, Stuart, and Cassie Lakes could not be determined. Clearly, though, the concentration of diatoms increased greatly in the recent sediments.

FIGURE 3.

Scatter plots comparing the amount of change observed among the three biological indicators using a Bray-Curtis dissimilarity coefficient (DC) for (A) diatoms versus chironomids, (B) diatoms versus cladocerans, and (C) chironomids versus cladocerans. A similar number of lakes plotted above and below the 1:1 line in (A) and a Wilcoxon-signed-rank test established that DCs for both diatoms and chironomids were not significantly different, suggesting that both biological indicators underwent similar magnitudes of change between top and bottom assemblages. The majority of lakes plotted above the 1:1 lines for both (B) and (C) and a Wilcoxon-signed-rank test indicated that changes in both diatoms and chironomids were significantly greater than changes recorded between tops and bottoms for cladocerans.

Similar to the diatoms, there was a clear separation in the chironomid assemblages of shallow study sites from deep sites represented along PCA axis 1 (Fig. 4, part B). One exception to this was the modern sample score for the shallow Opinnagau East Lake, which plotted closer to the sample scores of the deeper lakes. Similar to the diatom ordination, PCA_spp axis 2 for chironomids distinguished the shallow, coastal Billbear Lake from all other shallow lakes. In addition, deeper North Raft Lake clearly plots along PCA axis 2 and is distinct from the rest of the deeper study lakes within the ordination (Fig. 4, part B).

The largest change in chironomid assemblages from pre-disturbance to modern conditions was observed in the shallow Opinnagau East Lake, which has an assemblage shift from predominantly Tanytarsus spp. and Orthocladiinae littoral species to Chironomini littoral species and a few taxa also found in the deeper lakes (Fig. 4, part B). The remaining shallow lakes did not record a unidirectional change through time, compared to the deeper lakes (with the exception of North Raft Lake), which all appear to converge toward similar modern sample scores within the PCA_spp ordination space. Of the deeper lakes, Aquatuk Lake's pre-disturbance chironomid assemblage was most similar to the chironomid assemblages of North Raft Lake (modern and pre-disturbance) and displayed the greatest change through time, to become more similar in species assemblage to those found in Spruce and Hawley Lakes (Fig. 4, part B).

Shallow lakes plotted in the left quadrants of the PCA_spp chironomid ordination diagram and are largely composed of littoral genera associated with macrophytes, including Limnophyes-Paralimnophyes, Psectrocladius spp. (mainly P. sordidellus-type), Cladotanytarsus mancus-type, Orthocladius and Cricotopus spp., Paratanytarsus, Tanytarsus, and Microtendipes (Fig. 4, part B; Larocque et al., 2006; Brooks et al., 2007). Modern PCA sample scores for Opinnagau East Lake are also represented by chironomid taxa living in littoral, macrophytic habitats, but were distinguished from the other shallow lakes by a composition of genera mainly from the Chironomini tribe, primarily the warm-water species Cladopelma lateralis -type but also Dicrotendipes, Endochironomus, Polypedilum, and Glyptotendipes (Bennike et al., 2004). The modern PCA_spp sample score for Opinnagau East Lake also plotted closer to the deeper lakes of this data set due to the large increase in Corynocera ambigua (31%) from pre-industrial to modern samples. The deeper Hawley, Spruce, and Aquatuk Lakes plot in the right quadrants of the PCA ordination and are composed of mainly profondai taxa of the Chironomini tribe (Chironomus spp., Stictochironomus, Sergentia coracina-type) as well as cold-stenotherm taxa (Micropsectra, Corynoneura arctica-type, Corynocera oliveri-type) (Larocque et al., 2006). The more generalist taxa, Constempellina and Stempellina, were found in all deep lakes. North Raft Lake is distinct from the other three deep lakes based on the dominance of the cold-stenotherm species Heterotrissocladius marcidus-type, Micropsectra insignilobus-type, Tanytarsus chinyensiss-type, as well as the carnivorous Procladius.

The PCA_spp of cladoceran assemblages from the nine lakes containing cladoceran remains in both the top and bottom sediments also exhibited a clear distinction between the deep and shallow lakes along PCA_spp axis l (Fig. 4, partC). Unsurprisingly, the cladoceran taxa most strongly associated with the deeper lakes were the planktonic taxa Bosmina spp. and the Daphnia longispina complex (Walseng et al., 2003); however, changes in the assemblages through time (albeit subtle) were most strongly associated with PCA_spp axis 2 and littoral taxa such as Alonella excisa, Alona guttata, and Chydorus brevilabris (Walseng et al., 2003) in both shallow (Kinushseo and Opinnagau) and deep (Aquatuk and North Raft) lakes.

FIGURE 4.

Principal components analysis (PCA) for the first two ordination axes showing differences in species assemblages between top and bottom sedimentary intervals for the 13 study sites for (A) diatoms, (B) chironomids, and (C) cladocerans. The length and direction of the arrows between top and bottom samples indicate the magnitude and direction of change in the biological indicators for each site. Numbers in italics represent the position that each taxon plots within the species ordination plots and correspond to taxonomic names listed in Appendix Table A2. For ease of interpreting the PCA species figures, ellipses were drawn around distinct groupings of taxa for (A) diatoms and (B) chironomids. Pre-disturbance samples were plotted passively in the ordination. (D) PCA of measured environmental variables for the 13 study sites.

ASSEMBLAGE CHANGES THROUGH TIME

In this exploratory paleolimnological survey of HBL lakes, the three indicator groups we studied experienced differing degrees of directional change since pre-disturbance times. The greatest amount of change was observed in the most diverse and species-rich indicator (diatoms), closely followed by chironomids, while the least change was observed in the least diverse, taxon-poor indicator (cladocerans). The biological assemblages and measured environmental variables both indicated that lakes were most strongly separated along a depth gradient for each indicator group (Fig. 4, parts A-D). This is not surprising, given the large range in depth of our relatively small lake set. However, the nature of the differences between shallow and deep lakes suggest that habitat structure and availability, rather than depth per se, more readily explained differences between shallow and deep lakes in terms of species composition. Although the differences in biological assemblages were often greater between the shallow and deep lakes than between tops and bottoms (with a few exceptions), lake depth does not appear to be directly important for determining the magnitude of biological change over time (Fig. 4, parts A-C). Lake surface area did not appear to be linked to lake depth, as the deepest lake (Hawley Lake) and one of the shallowest lakes (Opinnagau East Lake) were the two lakes with the largest surface areas (Table 1). The shallow lakes close to the HB coast (Billbear, Stuart, Julison, Wolfgang, Sam, and Cassie) were also generally the smallest lakes in the study set (32-357 ha). In general, there did not appear to be a directional pattern in the amount of biological change and lake surface area.

Aside from lake depth and surface area, other measured environmental variables, such as specific conductance (and related variables), DOC and related variables, and nutrients, as well as distance to shoreline, were also important in explaining the distribution of taxa among these lakes (i.e., along a spatial gradient). However, no directional patterns were observed in terms of the magnitude of biological change through time (comparison of DCs between tops and bottoms), with lakes arranged along these measured gradients such as lake depth (data not shown). This suggests that our measured environmental variables may not have a large influence on the magnitude of biological change through time. The lack of clear trends through time along the environmental gradients may be due to a number of factors, including the small sample size of lakes, the relatively short measured gradients (with the exception of lake depth), little change in the measured variables over time, and/or that the changes in biological assemblages are in response to factors that were not measured in this survey, such as warming-induced changes in habitat availability and lake water properties.

Biological indicators in the bottom sediments of several HBL study lakes (Table A1) were scarce, and sufficient counts could not always be obtained. Dissolution does not appear to be a plausible explanation for the lack of diatoms, because the few valves that were encountered were well-preserved as were the fine spines from chrsyophycean stomatocysts (also scarce in these samples). In addition, the amount of siliciclastic material was substantially higher in the bottom samples, suggesting that, in the past, larger amounts of sedimentary material (or correspondingly less organic matter during colder periods) were delivered from the surrounding catchment, and therefore there were very few diatom and cladoceran remains relative to the siliciclastic material. Adequate chironomid remains were attainable in almost all samples but nevertheless required substantial amounts of material (often four times the standard counting effort) to obtain a minimum number of individuals required for analyses, and these were often poorly preserved. Nevertheless, a change from virtually no diatom remains preserved in the pre-impact samples to modern assemblages that contained an ample variety of diatom species, is ecologically significant, and consistent with decreasing ice cover (Smol and Douglas, 2007).

A possible alternative explanation for the scarcity of biological indicators in some of our HBL lakes may be related to distance to the HB coast and the high rates of isostatic rebound in this region. For example, the most shallow study lakes (<1.0 m), near the coast (e.g., Billbear, Stuart, and Julison), are located at the lowest elevation, have small surface areas (32 to 357 ha), are on continuous permafrost, have longer periods of ice cover in summer, and would likely freeze to the bottom in winter. Due to its recent emergence from the sea, the landscape near the HB coast is of a young age relative to sites farther inland, and therefore soil and peat development would be shallower (Rouse, 1991). Collectively these characteristics may play a role in explaining the scarcity of biological remains in the pre-warming sediments of these shallow HBL lakes, as they were likely ice-covered for most of the summer. Indeed, diatoms were lacking in the bottom samples of the four shallowest lakes Billbear, Julison, Stuart, and Cassie Lakes (the latter is located on discontinuous permafrost). Cladocerans were lacking in Julison, Stuart, and Cassie Lakes but were countable (although difficult to attain adequate numbers) in Billbear Lake (Table A1). However, whether or not a site contains adequate biological indicators in the bottom sediments appears to be site-specific, as nearby shallow, coastal Wolfgang and Sam Lakes contained ample numbers of all three biological indicators in both top and bottom samples, in contrast to other shallow, coastal sites. Wolfgang and Sam Lakes had the highest organic matter content (top and bottom intervals) in this lake set (organic matter = 30.9% to 73.6%, and 54% to 77.9%, respectively; Table A2); however, there was insufficient sediment for %OMC analysis for the remaining shallow, coastal sites. At present, it is not possible to provide a conclusive explanation for the lack of biological indicators in the bottom sediments of some of our HBL lakes. This is indeed a topic of interest that warrants further studies requiring a more extensive data set.

In contrast to the bottom samples, the majority of modern samples contained all three biological indicators in adequate numbers for analyses (with the exception of cladocerans in Stuart Lake), and %OMC was higher than in bottom samples in all lakes where LOI was measured. Although there was insufficient material to undertake sedimentary chlorophyll a analysis (e.g., Michelutti et al., 2005) for the top and bottom samples, detailed analyses of the dated sediment cores from the four deepest lakes indicated a distinct directional increase in whole-lake primary production (inferred from sedimentary chlorophyll a and its diagenetic products) from very low levels in the past, to moderately higher levels over the past few decades (Rühland et al., 2013). Additionally, declines in carbon/nitrogen ratios as well as γ¹³C analysis in the upper intervals of sediment cores from these same lakes also suggest recent increases in lake autochthonous production (Brazeau et al., 2013). Recognizing that %OMC is only a very approximate measure for estimating lake production (Smol, 2008), an increasing trend in all lakes may suggest that the study lakes are currently more productive than in the past. An increase in primary production is consistent with the onset of warmer conditions, longer open-water periods, and longer growing seasons for aquatic biota (Smol and Douglas, 2007; Vincent et al., 2013).

BIOLOGICAL COMMUNITY RESPONSES TO WARMING: ESTABLISHMENT OF PLANKTONIC TAXA, INCREASE IN LITTORAL HABITATS, AND CHANGES IN SPECIES DIVERSITY

Despite the particularly muted response of cladocerans as a group, a common trend observed in both the diatoms and the cladocerans was an increase in the abundance of planktonic/pelagic biota relative to benthic/littoral taxa in the post-warming samples. As expected, chironomids did not track changes in lake planktonic habitat directly, principally due to their predominantly sedentary larval life strategies on benthic habitats (e.g., rocks, lake bottoms, submerged wood, aquatic macrophytes, etc). However, they are known to closely track air and water temperature and associated water column changes (Eggermont and Heiri, 2012) as discussed below.

The nature of the biological changes observed in our HBL lakes is consistent with changes in water column properties and habitat diversification that are indirectly related to warmer air temperatures such as a longer open-water season (i.e., reduction in ice cover) and associated effects on light, lake thermal structure, mixing strength and depth, and nutrient cycling (Smol and Douglas, 2007; Rühland et al, 2008; Nevalainen and Luoto, 2012; Winder et al., 2012). An extension of the growing season, as well as changes in the thermal properties of the four deep lakes (10 to 35 m), is a likely explanation for this shift (Rühland et al., 2013). Increases in planktonic diatoms in these deeper lakes are consistent with findings from a top-bottom analysis of 50 lakes from the central Canadian subarctic (Rühland et al., 2003), as well as from 40 lakes in northwestern Ontario (Enache et al., 2011), where deeper (>∼6 to 10 m) lakes generally recorded the largest increases in small, Cyclotella taxa and other planktonic diatoms in the modern sediments. In high-latitude regions, even modest warming during a longer open-water season can induce a cold monomictic lake (i.e., typically having water temperature near the maximum density of 4 °C) to thermally stratify as water warms up above this critical density threshold (Vincent et al., 2013).

Consistent with our understanding of ice conditions and thermal stratification patterns (Smol and Douglas, 2007), Hawley Lake (the deepest study lake) exhibited only modest biological changes through time, particularly among the diatoms and the cladocerans (Fig. 4, parts A and C). This is a very long, large lake, formed by deepening and widening of the Sutton River and is therefore river-like in many respects. Although strong thermal stratification was recorded in 2001 (Gunn and Snucins, 2010), it is likely Hawley Lake still mixes frequently during the summer, but as warming continues and exposure of the water column to sunlight lengthens, the frequency and duration of stratification events may be expected to increase. This has been found in other regions of the Arctic, where deep, wind-swept lakes that commonly experience short periods of weak stratification during the summer will mix less frequently during warmer and longer open-water periods and may even shift from polymixis to dimixis (Milot-Roy and Vincent, 1994; Vincent et al., 2013). Given its morphometry and that it is the headwater lake of the Sutton River and therefore closely linked to the river, it is not surprising that Hawley Lake would take longer to exhibit a change in thermal stability and thus does not exhibit the magnitude of biological change recorded in other deep lakes. It may be that the modest changes observed in Hawley Lake today represent the very first signs of responses to warming.

Aside from the four deepest lakes (10.9 to 34.6 m deep), the remainder of the study lakes were shallow (0.7 to 2.5 m deep), and therefore a strong signal among the planktonic taxa was not expected. Nevertheless, even these shallow lakes registered increases in modern abundances (albeit low) of a variety of planktonic diatoms (Fig. 2, part A) as well as in the planktonic cladoceran Bosmina (Fig. 2, part C). In the shallow study lakes, it is likely that with a longer open-water period, the water column experiences a longer period of radiative heating, changes in light habitat, and associated changes that ultimately favored the onset of minor planktonic development for both diatoms and cladocerans.

In addition to the relative increase in planktonic/pelagic taxa, warming and decreased periods of ice cover will also lead to an increase in the availability of benthic/littoral habitats favoring the development of more complex and diverse biological assemblages among both algal and invertebrate groups (Quinlan et al., 2005; Douglas and Smol, 2010). In particular, the distribution of chironomids is closely linked to climate (Eggermont and Heiri, 2012), and thus inferences of regional warming in the HBL may be deduced from the chironomid assemblages. Unlike the planktonic shift among diatom and cladoceran taxa observed in the HBL study lakes, chironomid assemblage shifts occurred mainly within littoral taxa. In general, there were higher modern abundances of Cladotanytarsus mancus (p = 0.01), a chironomid often associated with warmer conditions (Barley et al., 2006) and lower modern abundances of the littoral Orthocladiinae genera Cricotopus (p = 0.007), and in some lakes, Psectrocladius (mainly P. sordidellus-type) (p = 0.15), which tend to be more common in cold-water environments (Eggermont and Heiri, 2012), especially in higher latitude lakes located along the eastern HBL (Larocque et al., 2006). Also, there were higher abundances in the modern samples of Corynocera ambigua (p = 0.06), a chironomid taxon whose ecology is poorly understood (see Barley et al., 2006, and Medeiros and Quinlan, 2011). This chironomid taxon may be partially driving the large PCA_spp sample score change between pre-disturbance and modern conditions in Opinnagau East Lake. We interpret the other chironomid assemblage shifts in Opinnagau East Lake from a pre-industrial assemblage of Orthocladiinae and Tanytarsus to a modern assemblage composed of littoral taxa from the Chironomini tribe and Tanypodinae (Fig. 4, part B), generally found at intermediate to warm temperatures and less abundant in cool lakes (Medeiros and Quinlan, 2011; Eggermont and Heiri, 2012), as a response to regional climate warming.

In many of the HBL lakes, including the four deepest lakes (Hawley, North Raft, Spruce, and Aquatuk), as well as the shallow Wolfgang and Opinnagau East Lakes, pre-disturbance diatom assemblages were commonly dominated (up to 90%) by small benthic fragilarioid taxa (Fig. 4, part A). These epilithic (attached to rocks) diatoms are known to be opportunists and can compete well under a variety of challenging environmental conditions, including extended periods of ice cover and light limitation (Smol, 1988; Lotter and Bigler, 2000; Douglas and Smol, 2010). Although fragilarioid taxa often remained dominant, modern benthic assemblages were distinctly more species-rich and included the appearance of a variety of periphytic forms (e.g., Cymbella, Brachysira, Achnanthes, and Navicula species). This shift toward a diversified and complex benthic community at all sites is consistent with an increase in littoral habitat availability and variety as would be expected with warming (Douglas et al., 1994; Smol et al., 2005). An increase in species diversity and richness was particularly striking in shallow Wolfgang Lake and in deeper Aquatuk Lake (Fig. 4, part A). Although the trajectory of change for Wolfgang Lake closely follows PCA axis 1 (depth gradient), it is likely that this shift toward diverse benthic assemblages more commonly found in the shallow lakes is a reflection of the development of new aquatic benthic habitats rather than a decrease in lake depth per se. In the deeper Aquatuk Lake, the increase in species richness in the modern sediments also included the first occurrence of planktonic taxa (e.g., Asterionella formosa and Fragilaria tenera) and tychoplanktonic Aulacoseira subarctica and A. ambigua, and a variety of periphytic taxa (e.g., Amphora, small Navicula taxa, and Achnanthes lanceolata) that were not observed in the shallower study lakes. Therefore, the trajectory of change for Aquatuk Lake across PCA_spp space (Fig. 4, part A) did not follow the major axis of variation (i.e., lake depth, axis 1) but was nevertheless as ecologically striking as that for Wolfgang Lake.

Unlike High Arctic ponds where warming has been reported to trigger increases in diversity across multiple food web levels (Quinlan et al., 2005), in the HBL sites, diatoms were the only indicator group to register a significant (p = 0.01) increase in species diversity across sites in the modern sediments (Fig. 2), although not always of the same magnitude reported in High Arctic systems. Changes in species diversity with warming may not necessarily follow a predictable pattern in all lakes and in all regions (e.g., Rautio et al., 2000). The lack of significant change in species diversity observed for the chironomids and cladocerans (Fig. 2, b and c) may be due to a variety of reasons. Relative to the diatoms, these indicator groups are represented by a much lower number of taxa and therefore an expansion of littoral habitats with warming may not lead to a greater diversity of taxa to fill in these new niches, particularly if the pre-disturbance assemblages were dominated by taxa considered to be ecological generalists (e.g., Bosmina spp., Chydorus brevilabris, some Tanytarsus spp.; Kamenik et al., 2007; Keller and Pitblado, 1989; Walseng et al., 2003; Brooks et al, 2007). Given the recent onset of warming in the HBL region, the study lakes may be experiencing the initial stages of species responses to this new climate regime, and, to date, the diatoms have registered the largest and earliest response of the three indicator groups studied. This is likely because there are a great variety of specialist diatoms with different life strategies that are capable of occupying specific aquatic habitats (e.g., associated with vegetation, algae, sediment, rocks, etc.) that can develop with warming. Cladocerans tend to have more general habitat preferences (e.g., Walseng et al., 2008) relative to diatoms and chironomids, but diatoms have the greatest variety of habitat preferences of all three indicators. Therefore, an increase in the available aquatic habitat types may not result in an increase in the diversity of invertebrate communities to the same extent as among the diatoms (e.g., Douglas and Smol, 2010).

Cladocerans do not appear to be as sensitive in their response to warming in the HBL region as diatoms and chironomids. The minimal change and the lack of directional change that we observed for cladocerans in this study is similar to the muted cladoceran changes reported in a regional survey of 50 lakes in the central Canadian subarctic (Sweetman et al., 2008), relative to the striking directional response of diatoms to warming from these same lakes (Rühland et al., 2003). The results of the HBL study suggest that cladoceran responses to environmental change are not strongly linked to the recent changes in diatom assemblages, indicating that top-down controls did not play a role in these lakes. It would be expected that, as food sources for cladocerans increase with warming (i.e., increased primary production), that cladoceran abundance (e.g., Daphnia) may also increase as observed in a deep, subarctic Finnish lake (Rautio et al., 2000). In the HBL, it is likely a combination of the shallowness of the majority of the study lakes, the rather recent change in the regional climate regime (last few decades), and the naturally slower response of cladocerans relative to diatoms that best explain the observed muted response. Differences in the magnitude of change between modern and pre-disturbance biological assemblages in the HBL lakes suggest that the indicator groups were not strongly coupled (e.g., lakes with large diatom shifts had minimal cladoceran changes, an increase in periphytic diatoms was not matched with an increase in chironomids that scrape periphytic algae, etc.). The lack of a strong food-web effect in the HBL lakes is similar to trends reported in other multiple proxy studies where diatoms, cladocerans, and chironomids were found to respond independently across both spatial and temporal gradients (Heegaard et al., 2006; Rautio and Vincent, 2006; Sweetman et al., 2008). This suggests that currently top-down and bottom-up controls do not appear to play a role in organizing the biological communities of the HBL lakes and that each group of aquatic biota may be responding differently to climate-induced changes in lake properties.