The Athabasca Oilsands Region (AOSR) has the largest bitumen reserves in the world (Lynam et al., 2015), with development starting in 1967, intensifying during the 1980s, and achieving near exponential growth beginning in the 1990s (Schindler, 2010, 2013). As the industrial footprint in the AOSR expanded, emissions of carbon dioxide, sulphur oxides, nitrogen oxides and other pollutants increased (Aherne and Shaw, 2010; Kelly et al., 2010; Kurek et al., 2013), leading to concerns about the impacts of the AOSR industrial activities on the surrounding terrestrial and aquatic ecosystems (Schindler, 2013). Initial environmental focus was on the impacts of sulphur and nitrogen compounds as acidifying agents leading to lake acidification, as seen in eastern North America (Sullivan et al., 1990; Cumming et al., 1994).
To date, the analysis of fossil diatom assemblages to infer past lake-water pH trends (Smol et al., 1986) from boreal lakes in Alberta close to the AOSR, as well as further downwind in Saskatchewan, showed limited evidence of regional surface-water acidification (Hazewinkel et al., 2008; Curtis et al., 2010; Laird et al., 2013). Although, upland forest soils in the AOSR should be sensitive to acid deposition, due to their extremely low base cation weathering rates, previous studies indicate that the risk of acidification is largely mitigated by the high deposition of base cations, even within 20 km of the industrial center where acid deposition is very high (Watmough et al., 2014; Fenn et al., 2015). Instead, environmental concerns now focus on the potential eutrophication impacts resulting from deposition of reactive nitrogen, Nr (Fenn et al., 2015). Evidence from the sedimentary studies of lakes in close proximity to the AOSR indicated elevated abundance of diatom species typical of higher nutrient conditions in a number of the lakes (Hazewinkel et al., 2008; Curtis et al., 2010), as well as increased rates of diatom production from cores in basins downwind in Saskatchewan (Laird et al., 2013). Together these studies speculate that enhanced production may result from a complex interaction between climate change, internal nutrient loading, and deposition of nitrogen (N) compounds derived from the AOSR.
Anthropogenic deposition of N from agriculture and fossil fuel combustion has been linked to nutrient enrichment of lakes, particularly oligotrophic lakes in high-altitude and high-latitude regions (Elser et al., 2009a; Baron et al., 2011; Lepori and Keck, 2012). Other studies have also implicated a potential synergy between climate warming and nutrient stimulation of Nr deposition (Thompson et al., 2008; Bergström et al., 2013; Spaulding et al., 2015). The importance of the deposition of inorganic N (nitrate and ammonium) on aquatic ecosystems is dependent on many factors, including soil cover and vegetation composition which influence soil N content, microbial N2 fixation and export to lakes (Baron et al., 2011; Lepori and Keck, 2012). In particular, lakes with limited catchment vegetation may receive a higher proportion of deposited N and may respond more quickly to enrichment than water-bodies in forested catchments (Fenn et al., 2003; Baron et al., 2011). Low levels of N deposition have led to pronounced ecological changes in some oligotrophic alpine lakes in the U.S. Colorado Front Range of the Rocky Mountains, as revealed through both limnological (Elser et al., 2009b) and paleolimnological studies (Wolfe et al., 2001; Saros et al., 2005b). Despite potentially important effects, the role of AOSR-derived N as a nutrient stimulus on lakes in the boreal region of northwest Saskatchewan is not known.
Fertilization effects of Nr should be greatest in lakes where growth of phytoplankton is limited by the supply of N. Although there is a continuing debate on the relative role of N and P as controls of lake production and eutrophication (Lewis and Wurtsbaugh, 2008; Schindler, 2012; Schindler et al., 2016), many studies show that N can limit phytoplankton growth alone or in conjunction with supply of P (Elser et al., 1990; Downing and McCauley, 1992; Guildford and Hecky, 2000; Nydick et al., 2004; Bergström and Jansson, 2006; Lewis and Wurtsbaugh, 2008). A long-term study of Mirror Lake, New Hampshire, U.S. indicates it has remained nutrient poor for over 30 years and co-limited by N and P (Winter and Likens, 2009). A large survey of lakes in the U.S. estimated that ~18% were nitrogen limited and another ~19% were colimited by N and P (Baron et al., 2011), although nutrient limitation status may also vary among seasons (Morris and Lewis, 1988; Kolzau et al., 2014).
To address the potential importance of N deposition on the production of lakes downwind from the AOSR, we selected lakes in a priori defined groups of N-limited and P-limited lakes in regions known to receive enhanced Ndeposition from the AOSR (Environment Canada’s IMP, 2011a). Historical changes in fossil diatom assemblages in these sites were compared to those from reference Nlimited and P-limited lakes located outside of the potential AOSR-impacted zone to isolate the effects of Nr deposition from those of regional climate change. This study analyzed fossil diatom assemblages in sediment cores to evaluate whether historical changes in these primary producers could be related to nutrient enrichment or climate influences since ~1890 CE.
Lake selection and study sites
Study basins were selected from a dataset of 262 lakes sampled during September-October in 2007-2009 and 2011 by the Saskatchewan Ministry of Environment within ~300 km of Fort McMurray, and represent an expansion of the survey provided by Scott et al. (2010). Individual lakes were sampled in 1 to 4 years. Initial site selection was based on choosing sites with small to moderate lake size, ~9-250 ha (~76% of lakes); depth, ~4-29 m (~51% of lakes); and catchment area: lake surface area of ~3-50 (~83% of lakes). Secondary selection removed lakes with outlier chemistry (high dissolved organic carbon, DOC; total nitrogen, TN or total phosphorus, TP; alkalinity and specific conductivity), while lakes with <6 mg L–1 alkalinity were removed to reduce potential sensitivity to acidification (Laird et al., 2013). Study sites were located either within a zone modelled by Environment Canada (2011a, Fig. 3.4, p. 50) to receive AOSR emissions, < ~150 km from AOSR (termed potential impact lakes) or a reference region outside the modelled deposition zone, > ~200 km from AOSR (Fig. 1). Of the 43 lakes remaining after initial screening, basins were then identified as potentially N- or P-limited based on molar dissolved inorganic nitrogen to total phosphorus (DIN:TP) ratios: <3.5 DIN:TP (N-limited or N + P colimited; Bergström, 2010); and >7 (majority >10) DIN:TP (P-limited) in both impacted and reference regions. Morris and Lewis (1988) indicated from bioassay experiments that the most accurate index of limitation (out of 9 indices, including TN:TP) was based on DIN:TP, as bioavailable nitrogen is better approximated by DIN (ammonium and nitrate). Bergström (2010) and Kolzau et al. (2014) also found DIN:TP ratios to be the best predictor of nutrient limitation. The resulting suite of lakes were then divided into a 2 x 2 factorial design, with at least 4 lakes in each group: four reference N- or N+P-limited lakes (hereafter Ref N-limited), four Ref P-limited lakes, four impact Nor N+P-limited lakes (henceforth Imp N-limited) and four Imp P-limited lakes. Lakes that were only sampled once, or had high variance in the multiple estimates of DIN:TP were not included in the final suite of lakes, except in the case of the impact N-limited lakes, as preliminary screening greatly reduced the number of available lakes in this category.
A non-metric multidimensional scaling (nMDS) ordination of the physical (area, depth) and chemical characteristics (chlorophyll a, Chl a; ammonium, NH4-N, TN, TP, DOC, dissolved inorganic carbon, DIC; chloride, Cl; sulfate, SO4, silica, Si; alkalinity and pH) was used to identify outlier lakes (those furthest from ordination centroid), resulting in the selection of a total of 22 candidate lakes with similar characteristics in each of four groups (Tab. 1 and Supplementary Tab. 1). Study and alternate lakes (selected in case of poor chronological control) included: four Ref N-limited lakes, four Ref P-limited lakes, five Imp N-limited lakes (+ two alternatives) and five 5 Imp P-limited lakes (+ two alternatives). Due to poor weather conditions and travel restrictions, one of our northern Ref P-limited sites and the two alternate Imp Plimited sites were not sampled.
Regional climate of the study region is characterized as Subarctic, with mean annual temperature at ~2.3°C, mean annual precipitation from ~450 to 530 mm, with 300-320 mm of rain occurring between May and September (Scott et al., 2010). Vegetation in the region consists primarily of jack pine (Pinus banksiana)/ lichen forest with balsam fir (Abies balsamea) and trembling aspen (Populus tremuloides) in areas with well-drained mineral soils; in the more poorly drained regions black spruce (Picea mariana) dominates. Further details of the ecoregions, geology and vegetation can be found in Scott et al. (2010).
Sediment coring and sampling
Collection of cores took place in February (11 lakes) and April (8 lakes) of 2014 (Fig. 1). Sites were accessed via a helicopter and cored from the ice using a Glew gravity corer with 7.62 cm internal diameter tube (Glew et al., 2001). Core lengths varied from ~27-52 cm (Tab. 1) and were sectioned into 0.5-cm intervals upon return to base camp. Core samples were shipped to Queen’s University, Kingston, Ontario where they were stored in a 4°C cold room until further analysis. The weight of all sediment samples was determined to calculate the total weight of sediment prior to any sub-sampling.
Radiometric dating and sample resolution
Twenty-one intervals were prepared for 210Pb analysis and counted using gamma spectroscopy at the Paleoecological Environmental Assessment and Research Laboratory (PEARL) at Queen’s University, Kingston, Ontario with samples closely spaced in the upper 10-cm of sediment (every 1 cm) and more coarsely spaced in deeper deposits (every 2 cm) down to 30 cm. Samples were freeze dried and the dry weight and percent moisture in each sample was determined. For gamma spectroscopy, dried sediment was weighed into a plastic tube, sealed with epoxy, and equilibrated for a minimum of 2 weeks to allow 214Bi to equalize with 226Ra, the parent isotope of 210Pb. Activities of 210Pb, 137Cs, and supported 210Pb (via 214Pb and 214Bi) were determined for each sample using an ORTEC well detector, as described in Schelske et al. (1994). Unsupported 210Pb activities were used to estimate the chronology of the cores using the constant rate of supply (CRS) model (Appleby and Oldfield, 1978) using the MATLAB programme ScienTissME (M. Scheer, unpublished data) with calibration to IAEA standards 312 & 385. Unsupported 210Pb was calculated by subtracting supported 210Pb. Age-depth models were calculated in R version 3.1.2 (R Core Team, 2010) using a cubic spline with a monotonicity constraint.
Sediment samples from each core were dated prior to sampling for diatoms. The target sampling strategy focused on sediments deposited after ca. 1890 CE and at a temporal resolution of ~5 years. Cores from Imp N lakes 10Y and 6E had to have additional samples prepared and counted based on unsupported 210Pb not being reached by 30 cm.
For each analyzed sediment interval, ~0.2-0.3 g of wet sediment was subsampled into a 20-mL glass vial to which was added a 1:1 molar mixture of concentrated nitric (HNO3) and sulphuric (H2SO4) acid to remove organic matter. Samples settled 24 h before the acid above was aspirated and the residue rinsed to a constant pH with distilled water. Four aliquots of each sample were pipetted onto coverslips, and coverslips were air-dried overnight before heating on a warming plate and mounting with Naphrax® onto glass microscope slides. Diatoms were identified and counted along transects on slides using a Leica (DMRB model) microscope fitted with a 100x fluotar objective (numerical aperture of objective = 1.3) and using differential interference contrast optics at 1000 x magnification. Approximately 400 diatom valves were enumerated per slide. Diatoms were identified to the species level or lower, using the following taxonomic references: Krammer and Lange-Bertalot (1986, 1988, 1991a, 1991b), Lange-Bertalot and Melzeltin (1996), Camburn and Charles (2000) and Fallu et al. (2000).
Sedimentary concentrations of diatoms in each sample were determined following Battarbee and Kneen (1982). Briefly, an aliquot of a known concentration of microspheres was added to each diatom sample, prior to subsampling for mounting on coverslips. The microspheres were enumerated along with the diatoms and used to calculate estimates of number of diatoms per gram dry weight of sediment. The flux rates of total planktonic diatoms and specific planktonic taxa groups were calculated by multiplying the concentration (# valves per g dry weight) by the sedimentation rate (g cm2 y–1) from the CRS chronology model.
The degree of species turnover in the diatom assemblages within each core since ca. 1900 CE was assessed using two methods. The first method used the gradient length of a detrended correspondence analysis (DCA) (Smol et al., 2005) using Canoco 5 (Šmilauer and Lep, 2014). The second method compared the community similarity of each sample assemblage to that from sediments deposited ca. 1900 CE using the Bray-Curtis similarity coefficient on diatom taxa present at an abundance of both 2 and 5% using Primer version 6.1.11 (Clark and Warwick, 2001). The ca. 1900 CE samples in the cores ranged from an age of 1892 to 1910 (1922 for Lake 2L), with a mean age of 1902 CE.
All reference lakes (both N- and P-limited) are located within the Boreal Shield Ecozone (Tab. 1, Fig. 1). The majority of the impact lakes (N- and P-limited) are also located within the Boreal Shield Ecozone, with three located within the Boreal Plain Ecozone. The majority of the lakes have surface areas of between ~20 to 90 ha, with 1 lake having surface area of ~9 ha and four that are larger at ~130 to 170 ha (Tab. 1). Most lakes are relatively shallow, with maximum depths of ~5 to 11 m and four slightly deeper at ~15 to 29 m. Lakes are circumneutral (pH 6.8-7.3), with low DOC (range 1.8-7.1 mg L–1), low alkalinity (4.2-18.1 1 mg L–1) and low specific conductance (13.8-35.5 μS cm–1).
Phosphorus levels of the lakes were within the oligotrophic to mesotrophic range (Wetzel, 2001), with estimated mean TP of ~3-18 μg L–1 (Tab. 1). The Ref and Imp P-limited lakes were within the oligotrophic range with mean TP of ~3-8 μg L–1 and the Ref and Imp N-limited lakes were within the mesotrophic range with mean TP of ~12-18 μg L–1. Mean Chl a data from the P-limited lakes (both Ref and Imp) were within the oligotrophic range (Wetzel, 2001), with overall mean Chl a of 2.1 mg m–3 (range 1.3-3.1) for the Imp P lakes and overall mean Chl a of 2.4 mg m–3 (range 2.0-3.2) for the Ref P lakes (Tab. 1). Whereas the N-limited lakes (both Ref and Imp), fall with the mesotrophic range, with overall mean Chl a of 5.2 mg m–3 (range 4.4-6.2) for the Ref N lakes and mean of 6.0 mg m–3 (range 3.9-7.3) for the Imp N lakes.
Chronology and resolution
The majority of cores in the study lakes showed an exponential decline in total 210Pb activity with sediment depth (Supplementary Fig. 1) enabling the development of strong depth-time chronologies. The coefficient of determination (r2) ranged from 0.69-0.99, with only 3 lakes with an r2 <0.8 (lakes 2L, 10E, 10Y). Background (or supported) levels of 210Pb (total 210Pb activity equals total 214Bi activity) were reached by 18-40 cm (median = 22 cm) in the N-limited lakes (Ref and Imp), and by 5-20 cm (median = 12 cm) in the P-limited lakes (Ref and Imp). Sediment accumulation rates were generally much higher in the N-limited lakes. The relationship between log unsupported 210Pb and cumulative dry mass was generally strong (r2 = 0.77-0.96). Although two lakes (Ref P- 2L and Imp N-10E) had a weaker relationship between unsupported 210Pb and cumulative dry mass (r2 = 0.53-0.56), the dating models are sufficient to access whether changes in species composition were concomitant with development of the Athabasca Oil Sands.
The number of samples analyzed since ca. 1890 CE varied between the cores due to the varying sediment accumulation rates. The median number of post-1890 CE samples in the N-limited lakes was 25 (range 21-36) and 20 (range 9-25) in the P-limited lakes. The average resolution between samples for the majority of the lakes (12 of 15) ranged from ~3 to 6 years (median = 5 yrs), with slightly lower resolution in Imp P-13L, at an average resolution of 7 years between samples. Only two lakes, Ref P-17F and Ref P-17P, were not analyzed at a sub-decadal resolution due to low sediment accumulation rates, with average resolution of ~10 and 14 years, respectively.
Diatom species assemblages and similarity analyses
An average of ~160 diatom taxa (range 94-269) were identified within each of the study lakes; however, the majority of these were rare, with only 13-55 (average 26) reaching >2%, and 7-22 taxa reaching abundances greater than 5%. The planktonic component of the assemblages varied across the lakes from ~ 3% in Ref N-2F to 87% in Imp P-11G (Fig. 2). There were few (if any) notable trends in the total % of planktonic taxa (Fig. 2), total planktonic concentrations (Supplementary Fig. 2), total planktonic flux (Supplementary Fig. 3), or in the dominant benthic taxa (Supplementary Figs. 4 and 5), either within or between the four lake groups defined based on DIN:TP ratios and geographic position either inside or outside the depositional area. Total planktonic concentrations for each lake are highly related to the total planktonic flux (r = 0.6-0.99), with the exceptions of lakes 2Q, 2L and 10Y (r = 0.2-0.44). A few lakes did record small variations in the taxonomic relative abundance of the dominant planktonic species in more recently deposited sediments (Fig. 3, Supplementary Figs. 4 and 5).
The predominant planktonic diatom taxa present in the cores varied between the lakes (Supplementary Figs. 4 and 5). A majority of the lakes (12 of 15) contained small centric Discostella (mostly D. stelligera (Cleve & Grunow) Houk & Klee and some D. pseudostelligera (Hustedt) Houk & Klee), ranging from ~2 to 54% of the assemblage with only 5 lakes > ~20% (Fig. 3A). Other Cyclotella (michiganiana Skvortzow, bodanica lemanica (Müller ex Schröter) Bachmann, tripartia Håkansson, ocellata Pantocsek) were typically present at lower percentages, except in Imp P-11G where C. tripartia comprised ~20-40% of the assemblage and in Imp P-13L where C. ocellata was ~10-20% (Supplementary Fig. 5d). Other planktonic taxa included representatives of the genus Aulacoseira (ambigua Simonsen, subarctica (Müller) Haworth, distans (Ehrenberg) Simonsen group, crenulata (Ehrenberg) Thwaites, tenella (Nygaard) Simonsen, perglabra (Øestrup) Haworth), although these species rarely reached >20% of the assemblage, except in Ref N- 17J in which A. ambigua and A. subarctica when combined comprised ~30-40% of assemblage (Supplementary Fig. 5a) and in Imp N-9C in which A. ambigua comprised ~20-30% (Supplementary Fig. 5c). A number of lakes also contained pennate planktonic taxa (Asterionella formosa Hassall, Fragilaria tenera (Smith) Lange-Bertalot, Tabellaria flocculosa (Roth) Kützing strain 3p sensu Koppen), which typically comprised less than ~10% of the assemblage, except the recent sediments in Imp N-10E which reached >15% (Fig. 3B, Supplementary Fig. 5c).
Fluctuations in the relative abundance of planktonic species were not consistent within or between the four lake groups in either direction or timing, but examples are noted, as comparable changes have often been related to recent climate change. The most notable Discostella increases of ~10% were recorded after ca. 2000 CE in Ref N-17V and Imp P-13N (Fig. 3A, Supplementary Fig. 5). Similar percent increases occurred in Imp N lakes 9C and 10Y prior to ca. 2000 CE (ca. 1980-2000 CE and ca.1930-1950 CE, respectively), while declines in Discostella were recorded in the Ref P lakes of variable timing (Fig. 3A, Supplementary Fig. 5). The most apparent increases in pennate planktonics occurred in Ref N-17J (A. formosa and T. flocculosa str. 3p) and Imp N-10E (Fragilaria crotonensis Kitton and T. flocculosa str. 3p) lakes after ca. 2000 CE (Fig. 3B, Supplementary Fig. 5). Smaller relative increases in A. formosa were recorded in Imp N-6E (after ca. 2000 CE) and Imp P-10A (after ca. 1950 CE), with very small increases in F. tenera in Imp P-11G after ca. 1970 CE. While declines in pennate planktonics occurred in Ref N-2Q following ca. 2000 and in Imp N-10Y after ca. 1980 CE. Planktonic Fragilaria (tenera, nanana Lange-Bertalot, capucina gracilis (Øestrup) Hustedt) comprised part of the Ref N-2Q assemblage and T. flocculosa str. 3p and F. tenera in Imp N-10Y (Supplementary Fig. 5). Concentrations of the planktonic species groups generally followed the trends of the percent abundances (Supplementary Figs. 6 and 7) indicating that the changes in percent composition are generally representative of the fluctuations recorded in the core assemblages. Rates of flux of the planktonic species groups generally reflect the same patterns seen in the individual lake concentration data (r = 0.68-0.99; Supplementary Figs. 8 and 9), with exception of lakes 2L, 2Q, 10Y and 17F (r = 0.26-0.55).
Small Fragilariaceae Greville were the predominant benthic taxa in many of the lakes (Supplementary Figs. 3 and 4), reaching combined average abundances in each core from ~5 to 79%. Staurosirella pinnata (Ehrenberg) Williams & Round was often the most predominant small Fragilariaceae, with an average in each core ranging from ~4 to 74%. Five of the lakes reached abundances >50% of S. pinnata; all were N-limited lakes, two reference (2F and 17V) and three impact (6E, 10E, 10Y). There was no significant relationship between max depth of the lake and the percent abundance of the small Fragilariaceae (r = -0.44, P=0.10), although there was a significant negative relationship of total benthic taxa and maximum depth (r = -0.62, P=0.013). Benthic taxa comprised > ~70% of the diatom assemblage in lakes <10 m maximum depth, except for Ref N-17J which only achieved ~40% benthic taxa.
Variability in the diatom species composition of the sediment assemblages was large across all lake groups, even with commonality of many taxa (Supplementary Figs. 4 and 5). Despite this pattern, there were few marked changes in diatom assemblages in individual lakes since ca. 1890 CE. Analysis of species turnover based on the DCA estimate was rarely greater than 1.0 (Supplementary Tab. 2) and, for most lakes, 20th and 21st century assemblage composition was highly similar (>80%) to that recorded in ca. 1900 CE based on the analysis of species with either 2% (not shown) or 5% relative abundance (Fig. 4). Only three lakes had any samples below 70% similarity, Imp N-10Y and Imp P lakes 10A and 13N. Lake Imp N-10Y had a few samples below 70% similarity in the 1940s as a result of a decrease in S. pinnata and an increase in Discostella (Supplementary Fig. 4). Lake Imp P-13N dropped to ~55% similarity in the 1940s as the result of an increase in S. pinnata and other small Fragilariaceae, and a decrease in D. stelligera (Supplementary Fig. 4). Lake Imp P-10A had the lowest similarities of all lakes due to greater abundance of S. pinnata and Aulacoseira distans humilis in the early 1900s, later fluctuations in S. pinnata and increases in A. formosa after ca. 1950 CE (Fig. 4).
Analysis of our study lakes in a factorial design (AOSR-impacted, reference; N- and P-limited) revealed few effects of either industrial development or climatic warming on diatom assemblages in northwest Saskatchewan, Canada. Analysis of metrics of diatom community change (% planktonic taxa, similarity analysis, species turnover) showed only minimal variation in assemblage composition in most study lakes (Figs. 2 and 4). The relative proportions and concentrations of many taxa were relatively stable over the past ~125 years (Fig. 3; Supplementary Figs. 4-7), other than generally small but nonsynchronous changes in centric taxa (Discostella and Cyclotella species) and pennate taxa, such as Asterionella formosa and Fragilaria crotonensis (Fig. 3), possibly reflecting individual lake response to climate. The absence of a consistent pattern of diatom changes associated with receipt of reactive nitrogen or intrinsic nutrient-limitation status of the lake suggests that AOSR emissions had no demonstrable effect on diatom composition.
Lack of large changes in diatom assemblages
Limited response of diatoms to the presumptive AOSR development may reflect the relatively low input of Nr to the boreal region of northwestern Saskatchewan. Distinct ecological response of diatoms to atmospheric deposition of industrially-derived N has been recorded in high-elevation alpine lakes in watersheds where limited terrestrial vegetation may favour N export directly to lakes (Wolfe et al., 2001; Saros et al., 2005b; Elser et al., 2009b; Baron et al., 2011). Estimated critical loads for diatom response in these Rocky Mountain basins range from ~1 to 3 kg N ha–1 yr–1 (Baron, 2006; Baron et al., 2011; Saros et al., 2011; Sheibley et al., 2014), whereas strong N effects in well vegetated eastern North American catchments have been estimated to require critical loads of ~3.5-6 kg N ha–1 yr–1 (Baron et al., 2011). These higher critical loads to elicit a change in aquatic ecosystems within a well-forested basin may reflect local uptake and storage of industrial N by N-limited catchment vegetation, as soil cover and vegetation composition has been found to influence N saturation and leaching (Baron et al., 2011; Lepori and Keck, 2012).
The forested boreal region of northwestern Saskatchewan exhibits higher land cover than alpine catchments, consequently atmospherically deposited N is more likely to be retained within the watershed (Environment Canada IMP, 2011b). In general, atmospheric nitrate and ammonium deposition from AOSR activities is greatest within ~25-30 km of the industrial center (Laxton et al., 2012; Proemse et al., 2013; Watmough et al., 2014; Fenn et al., 2015). To date, N loads (together with sulfur, S) have not caused widespread lakewater acidification due largely to mitigation by concomitant deposition of base cations (Watmough et al., 2014). Nitrogen cycling studies in Pinus banksiana and Populus tremuloides stands, suggests that soil N leaching was negligible for both types of forest plots, even within close proximity to the industrial center (Laxton et al., 2012). Vegetation throughfall data of DIN deposition was estimated at 2.1 kg ha–1 yr–1 ~25 km from the industrial center, which decreased to ~1.1 kg ha–1 yr–1 at sites located 113 and 129 km away from the main source (Fenn et al., 2015). These distant loads are similar to critical loads found to enrich high-elevation lakes within watersheds that have limited vegetation, but may be insufficient to saturate terrestrial plant demand and favour N export to Saskatchewan lakes downwind from the AOSR.
Limited diatom response to AOSR emissions could be a reflection of benthic diatoms being an important component of many of the fossil assemblages in our study basins (Fig. 2, Supplementary Figs. 4 and 5). The majority of our study lakes were dominated by benthic diatom taxa, with lakes generally having a planktonic component of <40%, with only two lakes > ~60% planktonic taxa (Imp N-9C and Imp P-11G). While littoral zones tend to be areas of high production and biodiversity, benthic algal production is often limited by light rather than nutrients, and may not respond as readily to nutrient fertilization (Vadeboncoeur et al., 2001; Cantonati and Lowe, 2014; Vadeboncoeur et al., 2014). The high complexity of benthic habitats and the close proximity of benthic algae to dissolved and particulate organic matter in the sediments has been suggested to be the main reason for this pattern (Saros et al., 2005a, 2005b; Bennion et al., 2001; Davidson and Jeppesen, 2013). Although epilithic algal communities have been noted to be influenced by water-column nutrient concentrations (King et al., 2000; DeNicola et al., 2004), estimates of total phosphorus from diatom assemblages dominated by benthic taxa have been found to have weaker responses to nutrient levels (Cumming et al., 2015).
Limited diatom response to Nr emissions may also be a reflection of inter-annual or seasonal changes in limiting nutrients (Morris and Lewis, 1988; Axler et al., 1994; Bergström et al., 2008; Kolzau et al., 2014; Dolman et al., 2016). For example, unproductive Swedish lakes in regions of higher deposition of Nr were pushed into Nlimitation only as the DIN pool declined over the summer season (Bergström et al., 2008). The degree and duration of lake stratification, as well as lake depth, can also influence the seasonal progression of limiting nutrients. Shallow lakes and those tending towards polymictic regimes can be limited by N (or light) in the summer months as a result of denitrification, and release of P from the sediments in summer (Kolzau et al., 2014; Dolman et al., 2016). Whereas, in deeper dimictic or monomictic lakes, a seasonal shift from P to N limitation may be limited by the isolating effect of stratification and any release of P from the sediments being trapped in the hypolimnion (Kolzau et al., 2014). While our deeper lakes (max depths ~15-29 m) were generally P-limited, the majority of our study lakes have max depths ranging from ~5-11 m, and were either N- or P-limited (Tab. 1).
Muted response of diatoms to atmospherically deposited Nr has been recorded in previous studies. For example, Spaulding et al. (2015) recorded synchronous, but small increases (~1-3%) of either A. formosa or F. crotonensis in relatively shallow (max depth ~7-8.5 m) alpine lakes, which were concurrent to estimates of Nr deposition from delta δ 15N in the sediments. In contrast, the diatom assemblages in their shallowest lake (max depth ~3.7 m) did not show similar changes in diatom assemblages (Spaulding et al., 2015).
The planktonic taxa, A. formosa and F. crotonensis have been identified previously as indicators of anthropogenic Nr deposition in a number of alpine lakes (Wolfe et al., 2003; Saros et al., 2003, 2005b). Increases in A. formosa and other pennate planktonics have also been suggested to be related to changes in the thermal stability of lakes (Rühland et al., 2015), as well as being linked to early increases in anthropogenic lake eutrophication (Anderson et al., 1995).The generally small increases of these planktonic taxa in our study lakes, with exception of the larger rise in Ref N-17J and Imp N-10E, were asynchronous and occurred across the different lake groups, suggesting that N deposition did not drive these changes. Interestingly, increases generally occurred in the more mesotrophic lakes with smaller surface areas (~9-41 ha vs ~60-170 ha), thus nutrient status and physical aspects of the lakes cannot be discounted as factors influencing these changes.
Is there a climate signal in the small phytoplankton changes?
Atmospheric mean annual temperature has increased during the 20th century by ~2°C in both northern Alberta (Laird et al., 2013) and northern Saskatchewan, particularly after the mid-1970s (Adjusted and Homogenized Canadian Climate Data archive, http://ec.gc.ca/dcchaahccd). These warming trends are consistent with those observed in other regions such as the Canadian Prairies and the boreal transition zones to the south of the present study (Schindler and Donahue, 2006). Analysis of seasonal temperature records from boreal regions of northwest Ontario indicated that the largest increases occurred during the winter months (Parker et al., 2009), with the result that more winter precipitation is now falling as rain (Schindler and Donahue, 2006) leading to changes in the timing and magnitude of snow melt and stream discharge. In principle, such changes in temperature and precipitation could alter lake production even without an increase in the deposition of N.
Longer growing seasons and higher water-column temperatures may manifest as increases in phytoplankton production in boreal lakes (Karlsson et al., 2005; Bergström et al., 2013), whereas cold temperatures can constrain growth even when excess nutrients are available (Thompson et al., 2008; Bergström et al., 2013). In this study, we found little evidence for coherent changes in the timing, magnitude or direction of change in small centric taxa (Discostella, Cyclotella), where previous studies have linked recent increases to climate warming (Rühland et al., 2008, 2015). Several of our study lakes recorded increases in both percent abundance and concentrations of Discostella taxa (Ref N-17V, Imp N-9C and 10Y, Imp P-13L and 13N), but timing of these changes were variable, and other lakes indicated declines in the relative abundance of this taxon (Ref P-17P and 2L). Rates of flux of Discostella were also variable within and between the lake groups. Similarly, there was no apparent pattern in which lakes indicated increases in Discostella in terms of morphological characteristics of lake surface area and maximum depth. While climate warming has affected surface- water temperatures and duration of stratification across the globe, these affects can be modulated by differences in lake morphometry, ice cover, geomorphic factors and variability in climate factors (O’Reilly et al., 2015; Kraemer et al., 2015). A global assessment of warming of lake-surface waters indicated a high degree of spatial heterogeneity with lake characteristics strongly mediating the response to climate effects, whereas previous regional-scale studies suggested a higher coherency (O’Reilly et al., 2015). While there is a stronger coherence of ice-out dates among lakes within a region (Magnuson et al., 2004), differences in lake morphometry and variation in timing of planktonic blooms can result in highly variable responses of phytoplankton communities (Boeff et al., 2016). Thus, while many lakes around the world have experienced increases in small Cyclotella taxa (including taxa now in Discostella), thought to be the related to climate warming manifested through increased strength or duration of stratification (Rühland et al., 2008, 2015), other studies have shown a lack of coherency of change within a region (Perren et al., 2009; Hobbs et al., 2010; Saros et al., 2011, 2012; Boeff et al., 2016). While such high variability in response makes it difficult to identify a common causal mechanism (Vogt et al., 2011), incorporation of phytoplankton autoecology with individual study lake processes can improve on the ability to decipher potential mechanisms.
Variability in the seasonal timing of diatom blooms and relation to stratification varies across lakes (Boeff et al., 2016; Wiltse et al., 2016). For example, the timing of higher populations of D. stelligera occurred in two Experimental Lake Area boreal lakes from early spring into early summer (Wiltse et al., 2016) and occurred during both un-stratified spring and summer stratified periods in three Maine lakes (Boeff et al., 2016). Boeff et al. (2016) found there was very little relationship between ice-out records and the relative abundance of D. stelligera, suggesting that direct links between the timing of ice out and potential duration of summer stratification may not always be a primary driver. Furthermore, Wiltse et al. (2016) suggest that increases in the abundance of Discostella may be the result of increased duration of spring mixing with earlier ice out. Whereas in other lakes, the degree of thermal stratification may be the main driver of Discostella abundances (Saros et al., 2016). Nonetheless, changes in light and nutrient availability may also alter a direct relationship between D. stelligera and thermal stratification (Saros and Anderson, 2015). Our deepest lakes (>18 m depth) indicated very little change in the percentage of Discostella/Cyclotella and increases were variable in magnitude and timing in the other shallower lakes (Ref N-17V, Imp N-6E, Imp P-13N). Variability in the timing of blooms and changes through time of this taxon suggests caution in indiscrimately identifying common direct or indirect linkages to climate.
Recent increases in other planktonic taxa may be the result of interactions of climate influences with other lake chemical and biological processes. For example, increases of A. formosa in Lake Windermere were suggested to be related to nutrient enrichment, warming and increased overwintering populations (Thackeray et al., 2008). Whereas, direct links to changes in thermal stability are often difficult, as A. formosa has been found to be in high concentrations during both un-stratified and stratified conditions (Wolfe et al., 2003; Boeffe et al., 2016).
Analysis of diatom assemblages in a robust factorial design demonstrated that potential nutrient enrichment as a result of nitrogen deposition has not greatly altered the abundance or species composition of diatoms in boreal lakes downwind from the AOSR. Changes in diatom assemblages were minimal across all groups of the study design. Similarly, although temperature has increased in this region and may have had small effects on the diatom assemblages through influences of water temperatures on phytoplankton growth, water column stability, or other direct or indirect effects or interactions, there was a lack of regional coherence in the direction and timing of most species change, including Discostella and small Cyclotella, which have been previously related to climate changes. Similarly, whereas five small, shallow (< ~10 m depth and less than 40 ha area), and generally mesotrophic lakes exhibited recent increases in A. formosa or F. crotonensis, the specific factors which influenced these changes and the Discostella/Cyclotella changes are uncertain and require further study.