A multi-proxy paleolimnological survey was performed on 13 lakes in the Hudson Bay Lowlands (HBL) of northern Ontario in order to provide a regional analysis of recent environmental changes in this poorly studied sub-Arctic region. In contrast to the amplified warming experienced by most of the circumpolar Arctic since the mid-19th century, the climate of the Hudson Bay (HB) region has remained relatively cool and stable for hundreds of years. However, since approximately the 1990s, the HBL has experienced rapid and large increases in air temperature and declines in sea ice. Diatom, cladoceran, and chironomid remains preserved in the recent (surface) and pre-1850 sediments of 13 lakes were used to examine whether this new climate regime has resulted in species assemblage changes across multiple trophic levels. Our results indicate clear limnological responses to warming among the freshwater biota of HBL lakes; however, the magnitude of this change varied among both biological indicators and sites. As expected, diatoms exhibited the greatest degree of change, closely followed by chironomids, with relatively little change observed among cladoceran assemblages. Planktonic diatoms were more common in modern assemblages, often including plankters that were previously not recorded in the bottom sediments, and in fact all indicator groups recorded a change in benthic/littoral taxa in the recent sediments indicative of warming-induced increases in habitat availability due to decreased lake ice cover.
During the past century, high-latitude regions have experienced amplified anthropogenic warming trends relative to the rest of the world (ACIA, 2004; Kaufman et al., 2009). Rising Arctic temperatures have changed fundamental properties of lentic ecosystems in ways that often precipitate striking responses among freshwater biota (e.g., Smol et al., 2005). However, in contrast to other northern regions, the subpolar Hudson Bay Lowlands (HBL) have registered only minimal warming (or even cooling) since at least the mid-20th century (e.g., Chapman and Walsh, 1993). This notable lack of warming has been attributed to a negative feedback response of the moderating role of Hudson Bay (HB) sea ice on the regional climate (Rouse, 1991). Historically, western HB (bordering the HBL) has been the last part of the bay to become completely ice-free in the summer, thus maintaining a cooler climate relative to surrounding regions, including those at higher latitudes (Rouse, 1991; Gough et al., 2004). Paleolimnological studies east of HB have yielded algal and invertebrate profiles with little biological change over the past few centuries (e.g., Laing et al., 2002; Paterson et al., 2003), as well as over longer time scales (Ponader et al., 2002; Pienitz et al., 2004; Fallu et al., 2005)—trends consistent with the observed lack of recent warming. However, recent studies have concluded that the HB region crossed a major climate threshold during the mid-1990s, precipitating rapid and substantial increases in regional air temperature, reductions in HB sea ice, and changes in sea ice phenology, resulting in a switch to a positive feedback response (Gagnon and Gough, 2005; Hochheim and Barber, 2010; Bhiry et al., 2011).
Due to its southern latitude, the HBL's status as part of the Canadian Arctic is often overlooked, and it is vastly understudied relative to other polar regions (Lemelin et al., 2010). Nonetheless, the sensitivity of the HBL to potential future warming is of concern, given the vast amounts of carbon stored in the region's extensive peatlands (e.g., Tarnocai, 2006). The abrupt onset of warming over the past few decades is already having notable effects on the structure and function of many aspects of the HB marine ecosystem. Notable examples include changes to the productivity of the marine food web (Hoover, 2010), with cascading effects on high-trophic-level sea birds (Mallory et al., 2010), as well as the replacement of polar bears as the top predator by relative newcomers to the region—killer whales (Ferguson et al., 2010; Peacock et al., 2010).
In fresh-water systems, the recent shift to warmer conditions in the HBL has already had notable ecological effects on fish populations over approximately the past decade. For example, in late summer of 2001, abnormally high air temperatures were observed (maximum 30 °C and exceeding 25 °C for at least 17 days) that resulted in uncharacteristically warm waters in the Sutton River and unusual thermal stratification conditions in upstream Hawley Lake (Gunn and Snucins, 2010). These unusually warm conditions led to thermal stress and a massive fish die-off event of anadromous brook charr (Salvelinus fontinalis) as they were traveling from the cold HB up the Sutton River to spawn. Although studies such as this are indicative that fresh-water systems of the HBL are already experiencing the effects of the recent rise in temperature, detailed temporal data are lacking. The paucity of data from the HBL is especially pronounced with respect to limnological studies. For example, with the exception of a few recent studies in Wapusk National Park, Manitoba (e.g., Bos and Pellatt, 2012; Symons et al., 2012; Bouchard et al., 2013), little information is available on the fresh-water ecosystems of the HBL. Available paleoenvironmental assessments are likewise sparse from the region. These include Holocene-scale, multiple-proxy paleoecological studies from a lake west of the HBL in Manitoba (Camill et al., 2012) and from a peat bog in the Attawapiskat River watershed in the southern HBL (Bunbury et al., 2012). Examinations of the effect of the recent and rapid onset of high-magnitude warming, and diminishing sea ice upon adjacent fresh-water aquatic ecosystems in the HBL have recently been reported from four high-resolution diatom records from the Sutton River region (Rühland et al. 2013); however, a regional survey from a variety of lakes (deeper and shallow) across multiple trophic levels has not yet been examined.
Here, we present an exploratory spatial and temporal survey of biological assemblages preserved in the sediments of 13 lakes located in the Sutton River region of the HBL, northern Ontario. We use a “top-bottom” paleolimnological approach to examine whether the present-day (top sediment) algal (diatoms) and invertebrate (chironomid and cladoceran) assemblages differ from pre-warming (bottom sediment) assemblages (Smol, 2008). The objectives of our study were to assess how lake ecosystems near the Sutton River in the HBL have responded to recent regional warming trends across taxonomic groups from different trophic levels. Given the relatively low number of study lakes (13), this study is exploratory in nature. Our questions include the following: (1) Can we detect temporal changes in the assemblage composition and species diversity of algal and/or invertebrate assemblages? (2) Do algal and invertebrate assemblages show similar temporal and spatial trends? (3) Is the magnitude of response similar across lakes and across trophic levels? If not, what are the possible reasons? (4) How do the patterns of biological change relate to measured environmental variables (e.g., distance from coastline, lake depth), and which lake types (e.g., deep, shallow, inland, coastal, etc.) register the greatest (and least) responses among the three indicators? (5) How do the trends we observe in the HBL compare to other lake systems where similar types of analyses have been completed?
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).
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.
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.
To ensure that the bottom sediment samples were deposited prior to the onset of warming in the HBL, 210Pb 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). 210Pb activities from these four lakes together with 210Pb 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 210Pb activities in the top and bottom intervals to 214Bi (a proxy for background or supported 210Pb activity levels) allowed for a determination of whether the bottom samples represent background conditions.
SEDIMENT SAMPLE ANALYSES: LOSS ON IGNITION, DIATOMS, CHIRONOMIDS, CLADOCERANS
Loss on Ignition
Loss on ignition (LOI) was used to estimate the percent organic matter content (Heiri et al., 2001) of both top and bottom sedimentary intervals for North Raft, Aquatuk, Hawley, Spruce, Wolfgang, Sam, and Opinnagau East Lakes. For each sample, approximately 0.3 grams of freeze-dried sediment was combusted at 550 °C for 4 hours (following Heiri et al., 2001). There was insufficient sediment to undertake LOI analyses for Billbear, Julison, Cassie, Kinushseo, Stuart, and Raft Lakes.
Sediment samples (∼0.3 g wet sediment if available, otherwise ∼0.03 g freeze-dried sediment) were prepared for diatom analysis using a strong acid solution (50:50 molar ratio of concentrated nitric acid [HNO3] to concentrated sulfuric acid [H2SO4]) following standard procedures (outlined in Battarbee et al., 2001). For many of the 13 lakes, the prepared diatom samples contained an excess of siliciclastic material and very low concentrations of diatom valves in the bottom (i.e., pre-warming) samples. To attain sufficient diatom counts for analysis, excess clastic material was removed and diatom valves were concentrated on the microscope slide by treating with a density gradient separation technique using the heavy liquid sodium polytungstate (SPT) (Tapia and Harwood, 2002) at a density of 2.3 g/cm3. Wherever possible, a minimum of 350 (usually more than 500) diatom valves were counted for each sedimentary interval using a Leica DMRB light microscope fitted with differential interference contrast optics at l000x magnification. Diatoms were identified to the lowest taxonomic level (often to variety) possible using a selection of taxonomic references including Krammer and Lange-Bertalot (1986–1991), Camburn and Charles (2000), and Antoniades et al. (2008). For each sedimentary interval, diatom counts were expressed as a percent relative abundance of the total number of valves counted. To eliminate very rare taxa, only diatoms that occurred in at least 1% relative abundance in at least two lakes were retained for analyses.
Chironomid head capsules were recovered and identified following the detailed procedures outlined in Walker (2001). For all top and bottom intervals, a subsample of wet sediment (typically 2 g but up to 25 g for a few samples with particularly high water content) was deflocculated in 5% KOH at 70 °C for 20 minutes and rinsed with deionized water through a 100-µm sieve. The sieved sample was rinsed with deionized water into a 100-mL beaker, and subfossil chironomid head capsules were picked from a Bogorov chamber using a dissecting microscope and forceps. Specimens were permanently mounted onto microscope slides using Entellan®. Our goal for each sample was to identify a minimum of 50 head capsules, the number deemed satisfactory for chironomid analyses (Larocque, 2001; Quinlan and Smol, 2001). Specimen identification to the lowest taxonomic level (often to the genus level, but whenever possible to the species level), primarily following Wiederholm (1983) and Brooks et al. (2007), was performed using a Leica DMRB light microscope with bright field optics at 400x magnification. For each sedimentary interval, chironomid head capsules were expressed as a percent relative abundance of the total number of head capsules enumerated. A cut-off criterion to eliminate rare taxa was not applied to the assemblage data.
Preparation of microscope slides for the cladoceran analysis followed standard procedures (Korhola and Rautio, 2001) using ∼1.0 g of wet sediment per sample if available, otherwise ∼0.01 g of freeze-dried sediment. Cladoceran remains were identified using a Leica DMRB light microscope with bright field optics (l00× to 400× magnification), and the primary taxonomic resources were Korosi and Smol (2012a, 2012b). The remains were tabulated separately by type (i.e., carapaces, headshields, postabdominal claws, etc.), and the total number of individuals represented was calculated using the most abundant remain for each taxon (Frey, 1986). To obtain a representative estimation of the cladoceran assemblage, the target count sum was 70 individuals (Kurek et al., 2010), and slides were counted in their entirety to remove any potential bias from a nonrandom distribution of remains. For each interval, cladoceran counts were expressed as a percent relative abundance of the total number of individuals counted. To eliminate the influence of very rare taxa, only those present at >2% relative abundance in at least two samples were included in the statistical 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 PCAspp . 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 PCAspp 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 PCAenv . 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 PCAspp and PCAenvp 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).
Radioactive decay of 210 Pb was at background levels (i.e., approached the 214Bi 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. PCAenv 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 PCAenv axis 2 represented a mixed gradient of conductivity (and related variables), dissolved organic carbon (DOC), and color. These environmental gradients in the PCAenv are consistent with the distribution of species in the PCAenv for all indicator groups, with PCAspp 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 PCAspp 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 PCAenv (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 PCAspp 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 PCAspp plot).
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.
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, PCAspp 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 PCAspp 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 PCAspp 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 PCAspp 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 PCAspp 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 PCAspp 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 PCAspp 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.
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 γ13C 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 PCAspp 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 PCAspp 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.
Multi-proxy paleolimnological analyses provide evidence of early responses among freshwater biota to the recent shift in the HBL climate regime. Although the magnitude of change was variable across sites, with no discernible difference between the deep and shallow lakes, diatoms demonstrated the most pronounced change, closely followed by chironomids, whereas only modest changes were observed among cladoceran assemblages. The most notable trend common to all indicator groups was a change in benthic/httoral taxa in the modern samples that was indicative of a warming-induced expansion of available aquatic habitats with decreased ice cover. In the deeper lakes, warmer winters and a longer summer growing season due to a shorter period of ice cover have likely increased the frequency and duration of thermal stratification events (already evident in increased abundances of planktonic diatom taxa). However, if rising temperatures in the HBL continue along the current trajectory, earlier, longer, and stronger periods of thermal stratification may result in more pronounced response across trophic levels, with increases in the relative abundances of pelagic cladoceran taxa and continued increases in planktonic diatom taxa. The lack of strong linkages among our three indicator groups suggests that top-down controls currently do not play a strong role in the HBL lakes and that these lakes are experiencing the very early stages of aquatic ecological change. However, with continued warming, it would be expected that higher magnitude biological changes will quickly follow. For some of the shallower lakes, warming may ultimately lead to complete desiccation, as already observed in some regions of Canada's High Arctic.
We thank Albert and Gilbert Chookomolin for their excellent guidance in our fieldwork and for their generosity in sharing valuable traditional knowledge and insights into the Hudson Bay Lowlands region; Joan Bunbury, Ron Ingram, Chris Jones, Shannon MacPhee, and Maara Packalen for fieldwork assistance; Hearst Air, particularly pilots Georges and Mike Veilleux, for their first-rate air transport service in this remote region; the staff at the Ontario Ministry of the Environment's Dorset Environmental Science Centre for chemical analyses and database support. This work was supported by the Ontario Ministry of the Environment through the Climate Change and Multiple Stressor Research Program at Laurentian University, the W. Garfield Weston Foundation Fellowship from the Wildlife Conservation Society Canada (to K. Hargan), and by the Natural Sciences and Engineering Research Council of Canada (to J. P. Smol).
The number of biological indicators counted in the 13 HBL lakes and the % organic matter content through loss on ignition (LOI).
Taxon numbers and corresponding taxonomic names for diatoms, chironomids, and cladocerans used in the principal components analysis figures (Fig. 4, parts A-C).