Paleolimnology of the Lake of the Woods southern basin: continued water quality degradation despite lower nutrient influx

ABSTRACT

Reavie ED, Edlund MB, Andresen NA, Engstrom DR, Leavitt PR, Schottler S, Cai M. 2017. Paleolimnology of the Lake of the Woods southern basin: continued water quality degradation despite lower nutrient influx. Lake Reserv Manage. 33:369–385.

Despite decades of reduced nutrient inputs, Lake of the Woods, a large, shallow boreal lake on the US-Canadian border, shows little evidence of water quality improvements in the pelagic system. Here we analyzed sediments from 6 sites in the southern basin for diverse biogeochemical (loss-on-ignition, biogenic silica, pigments) and microfossil (diatoms, chrysophytes) remains to reconstruct the environmental history of the lake. Our objectives were to quantify the magnitude and direction of historical trophic change and evaluate reasons for an apparent lack of basin recovery following documented nutrient diversion. Evidence came from fossil indicator profiles and comparisons of these long-term trends with historical land use and monitoring data. Results indicate major changes in algal communities during and following peak nutrient loading in the mid-20^th century as well as more recent increases in colonial cyanobacteria and high-nutrient diatom taxa. Combined, fossil indicators reflect an anthropogenically enriched system that has undergone substantial ecological change, particularly since ∼1980, due to multiple drivers. Physical changes in lake thermal regime resulting from climate warming may be exacerbating internal phosphorus release from sediments, thereby lowering nitrogen:phosphorus ratios and enhancing cyanobacterial abundance. These drivers of lake condition in the lake may apply to other large shallow lakes that exhibit only limited biological recovery from reduced external nutrient loading.

KEYWORDS:

• Cyanobacteria
• diatoms
• eutrophication
• Lake of the Woods
• nutrient loading
• paleolimnology
• phosphorus
• pigments

Cultural lake eutrophication is a widespread problem that contributes to the increased occurrence of algal blooms and hypoxia, impacting human and aquatic life. As a result of policies to reduce nutrient flux from agriculture, sewage treatment plants, and fossil fuel combustion, there has been some remediation from the ecological and economic effects of eutrophication (Schindler et al. 2016). Unfortunately, in many cases the reduction of nutrient inputs has not been sufficient to return lakes to pre-impact conditions. Exacerbating factors such as warming temperatures and nutrients stored in sediments have frequently prevented full recovery (McCrackin et al. 2016). We consider such scenarios for the southern basin of Lake of the Woods.

Lake of the Woods (49°N, 95°W) is a large (385,000 ha), shallow (mean depth = 7.9 m) lake occupying parts of the Canadian provinces of Ontario and Manitoba and the US state of Minnesota (Fig. 1). Water movement is from south to north from the Rainy River to the outlet near Kenora, Ontario (Anderson et al., 2017), so human activities in the southern catchment are likely the dominant drivers of water quality in the southern basin. Historic and present land uses around the lake include agriculture, urbanization, recreation such as tourism, and industry including pulp and paper, mining, and forestry. In the southern regions, the effects of European settlement became notable after ca. 1905 in Lake of the Woods County, Minnesota, particularly in the southern watershed of Zippel Bay (Baldwin 1926). Forest clearance in northern Lake of the Woods County began around this time, and by 1915 a number of roads and drainage ditches had been constructed in the county. By 1925, farmland in the catchment totaled ∼65,000 ha. Human manipulation of water levels on the lake was initiated with the construction of the Norman Dam in 1895 at the outlet near Kenora. Not surprisingly, these anthropogenic activities in the catchment had profound impacts on water quality (Clark et al. 2014).

Figure 1. Core locations in Lake of the Woods, USA and Canada.

Records from the 1800s suggest Lake of the Woods is naturally productive (Anderson et al. 2017), a probable result of a southern catchment composed of nutrient-rich glaciolacustrine sediments. Examination of nutrient loads to Lake of the Woods indicates that during much of the 20^th century, loadings came from point sources, mainly municipalities and pulp mills along the Rainy River, which today contribute two-thirds of the annual phosphorus (P) load and 70% of the water supply to the lake (Hargan et al. 2011). Past nutrient additions likely increased primary production and degraded recreational and aesthetic values. Following the US Clean Water Act of 1972 (Litke 1999), inputs of nutrients from municipal sources and pulp mills to the lake were reduced (Hargan et al. 2011, Lake of the Woods Water Sustainability Foundation 2011), the P load from the Rainy River decreased from ∼1500 t/y (1960–1975) to ∼500 t/y (post-1990; Hargan et al. 2011), and limited monitoring suggests a decline in lake P concentrations (Fig. 2). Despite this reduction, symptoms of eutrophication are still observed with large populations of toxin-producing cyanobacteria in recent years (Chen et al. 2007, 2009), particularly in late summer and fall (DeSellas et al. 2009). Although anecdotal, this recent trend is similar to observed increases in cyanobacteria recorded from lake sedimentary records throughout north temperate-subarctic regions (Taranu et al. 2015).

Figure 2. Long-term data from Lake of the Woods and the surrounding region. Plots (a–c) detail measurements of N (total Kjeldahl, total nitrates, total organic), TP, and the N:P atomic ratios. P data include measurements from the southern basin and Fourmile Bay, just downgradient of the Rainy River outlet. Ratios in plot (c) include TN and dissolved inorganic N (DIN) relative to TP. The ranges of data indicating limitation of algal growth by either nutrient (Bergström 2010) are also provided. Plot (d) indicates available data on N compound deposition from atmospheric monitoring stations near the lake. Plot (e) summarizes lake catchment population and agricultural activity. Plot (f) provides minimum, maximum, and mean annual atmospheric temperatures in the region. Plot (g) provides trends in total annual precipitation and ice-off day. Plot (h) provides historical wind speeds from a nearby station (International Falls) plotted with a running average. Plots (i) and (j) summarize the first 2 axis scores from a detrended correspondence analysis of diatom assemblages from each core. Lowess lines (plots d-h) are based on a span of 5 years, and DCA scores are connected (not fitted) to form a line for each core. Details of data sources are provided in Supplement A.

A lower atomic ratio of total nitrogen to total phosphorus (TN:TP) is believed to result in a greater annual abundance of N-fixing cyanobacteria (in certain Canadian lakes; Orihel et al. 2012). There is some evidence that such changes have occurred in the lake with the cessation of point-source discharges in the 1970s (Fig. 2), particularly in the ratio of dissolved inorganic N to TP (DIN:TP), which is a better indicator of nutrient limitation than TN:TP (Bergström 2010). In particular, a shift from P limitation to N limitation or co-limitation is suggested. It may also be that the historical P loads in the lake sediments have become increasingly available for ongoing algal production through the process of internal loading (remobilization of P from the lake sediments; Søndergaard et al. 2003). While not yet fully understood in Lake of the Woods, elsewhere this process has resulted in water quality responses to nutrient remedial measures that are slower than expected, or inconsequential, because of the effect of recycled sediment P supporting algal growth (Sharpley et al. 2013, McCracken et al. 2016). With the recent discovery of periodic stratification and reduced hypolimnetic oxygen concentrations in the southern basin during extended periods of low wind activity (Heathcote et al. 2016), it is possibile that climate change may now favor the release of P from sediments. Such phenomena may be exacerbated by increased atmospheric temperatures that can strengthen summer stratification (Nürnberg et al. 2013) and elevate microbial degradation of sedimentary organic matter.

Long-term environmental data are critical to proper management of the lake. Even with an ongoing monitoring program (e.g., Pla et al. 2005, Clark et al. 2014), few data quantify how human activities have altered environmental quality or forecast how continued impacts may cause future environmental changes. Consequently, management strategies lack a long-term perspective of environmental degradation or its causes. Lake of the Woods has been investigated using various paleolimnological techniques (e.g., Pla et al. 2005, Yang and Teller 2005, Reavie and Baratono 2007, Rühland et al. 2010), with other studies currently in progress, but to date little has been revealed on the vast, shallow, southern basin. To address this, we conducted a paleolimnological analysis of the lake using 6 sediment cores from the southern basin to reconstruct long-term changes in water quality and rehabilitation (if evident). Specifically, we sought to determine (1) the extent of past ecological change, (2) the effects of rehabilitation efforts in remediating water quality, and (3) the likely driving factors resulting in changes in phytoplankton abundance and community composition. We hypothesized that undesirable algal blooms are a recent phenomenon, despite a reduction in nutrient flux from anthropogenic sources and that the timing of water quality change may correspond to observed atmospheric warming and its effects on internal nutrient loading.

Methods

Sediment sampling

We collected 6 widely spaced sediment cores from the southern basin of the lake to reconstruct the post-European settlement history (ca. 1890). Sediment profiles were collected from the main basin (Big Traverse Bay, 2 stations), Muskeg Bay, Sabaskong Bay, Little Traverse Bay, and Big Narrows (Fig. 1). Surface cores ∼1 m long, preserving an intact sediment–water interface, were collected from the ice surface from 28 February through 1 March 2012 using a piston corer. The cores were sectioned at 1 cm intervals and refrigerated until analyses. On completion of analytical work, the remaining sediment material was freeze-dried for long-term archiving at the St. Croix Watershed Research Station (Marine on St. Croix, MN), and prepared diatom slides were archived at the Natural Resources Research Institute (Duluth, MN).

Inorganic and organic content, sediment dating

Analysis of organic matter content by loss-on-ignition, following Dean (1974), was performed on all sediment intervals. Sediment water content was determined from weight lost following oven drying of sediment aliquots at 100 °C for 24 h. Weight loss after ignition in a muffle furnace at 550 °C for 2 h was used as an estimate of organic content. Weight loss after combustion for 2 h at 1000 °C provided an estimate of carbonate content (largely calcium carbonate; Boyle 2001).

Cores were analyzed for isotopes of lead (²¹⁰Pb) activity by alpha spectrometry and cesium (¹³⁷Cs) and beryllium (⁷Be) activity by gamma spectrometry to determine age, sediment mixing depth, and accumulation rates during the past 100–150 years. Details of the dating procedures and results are provided in Edlund et al. (2017).

Pigments

Pigments and biogenic silica were analyzed from approximately every second sediment interval. Fossil pigments including carotenoids, chlorophylls, and their derivatives were extracted in a mixture of acetone:methanol:water (80:15:5, by volume) for 24 h in darkness and under an inert N₂ atmosphere at 4 C from freeze-dried sediments following the standard procedures of Leavitt and Hodgson (2001). Briefly, carotenoid, chlorophyll, and pigment-derivative concentrations were quantified by reversed-phase high performance liquid chromatography (RP-HPLC) using an Aligent 1100 HPLC system equipped with a C-18 column (5 µm particle size; 10 cm length) and an Agilent model 1100 photodiode array spectrophotometer (435 nm detection wavelength). Concentrations are reported relative to total organic carbon as nmole pigment per gram of organic matter (Hall et al. 1999).

Biogenic silica

Biogenic silica (BSi), a proxy for historical diatom and chrysophyte algal abundance, was measured using weighed subsamples (30 mg) from each core, which were digested using 40 mL of 1% (w/v) Na₂CO₃ solution heated at 85 °C in a reciprocating water bath for 5 h (DeMaster 1979, Conley and Schelske 2001). A 0.5 g aliquot of supernatant was removed from each sample at 3, 4, and 5 h. After cooling and neutralization with 4.5 g of 0.021 N HCl solution, dissolved silica was measured colorimetrically on a Lachat QuikChem 8000 flow injection autoanalyzer as molybdate reactive silica (McKnight 1991).

Siliceous microfossils

Sample digestion, plating, and analysis

Diatom fossils were analyzed from approximately every second interval in each core. Diatom frustules were cleaned of organic material by strong acid oxidation to allow identification of diatom species. Coverslips were prepared using the Battarbee (1986) method and adhered to slides with Naphrax mountant. Diatoms were identified and enumerated using a light microscope at 1000–1250 × magnification with oil immersion. At least 400 diatom valves were counted per slide. Diatoms were identified to the species level or higher using standard floras and iconographs. Although diatoms were the primary siliceous fossils, chrysophyte stomatocysts and phytoliths were also enumerated to provide additional ecological information. For each core depth, count data were converted to percent (%) relative abundance (relative to total valves counted) and accumulation rates (valves cm²/yr).

Temporal shifts in diatom assemblages

To aid in temporally grouping historical trends in the diatom assemblages, cluster analysis was applied to the fossil assemblages using the chclust function in R (statistical package version 2.15.2; R Development Core Team 2010) using the rioja package (version 0.8-7; Juggins 2014). The CONISS algorithm (Grimm 1987) was used to perform clustering constrained to vertical stratigraphy. The embedded function bstick was used to perform a broken-stick analysis and determine the minimum number of significant clusters (Bennett 1996).

We condensed the diatom assemblage data using nonmetric multidimensional scaling (NMDS) to produce a simple metric of temporal changes in assemblage composition, timing of community changes, and similarity in assemblages among cores. NMDS, as implemented in the R package vegan (Oksanen et al. 2013), compared sample assemblages based on distance metrics (Rabinowitz 1975) and ordinated the phytoplankton data using Bray-Curtis dissimilarity on relative abundance diatom data. Analyses were performed on all cores in a single analysis, constrained to 2 axes. Similarly, a detrended correspondence analysis (DCA) was performed to derive sample scores for the first 2 major ecological gradients (axes 1 and 2), so that long-term changes in these scores could be related to environmental monitoring data, including atmospheric temperatures, wind, precipitation, N deposition, dates of ice melt, and regional human population and agricultural activity (details in Supplement A). Relationships between environmental parameters and diatom assemblages were measured by calculating the Spearman correlation coefficients among stratigraphic time series. The sample dates varied in each time series, so for correlation analysis, environmental measurements were standardized using a running average across time. For instance, for a 1980 diatom sample, a running-average temperature extrapolation for 1980 was used to estimate a corresponding temperature value.

Diatom-inferred environmental conditions

Past TP concentrations were inferred from the fossil diatom assemblages using an inference model based on Minnesota lakes (Ramstack et al. 2003, Reavie et al. 2005, Paterson et al. 2007, Hyatt et al. 2011). Weighted averaging calibration and regression were used to derive a diatom-based model using R and the package rioja. The full model is known to have good performance statistics (Reavie and Juggins 2011, Reavie and Edlund 2013), and details of diatom-based transfer function development and application are provided by Juggins and Birks (2012). In addition, steps were taken to ensure the diatom-based TP training set model was appropriate for application to the diatom assemblages in the Lake of the Woods cores (Reavie and Edlund 2013, Reavie et al. 2015), including analysis of analog compatability and application of canonical correspondence analysis constrained to TP to confirm fossil samples fit the modern TP gradient.

Results

Dating and sediment accumulation

Radioistope datingFigure 3., Figure 4. of the 6 cores indicated conformable, though relatively slow, sediment accumulation throughout the deeper areas of the southern basin (Edlund et al., 2017). Linear accumulation ranged from 17 to 32 cm since ca. 1900 and dry mass accumulation from 0.02 to 0.12 g cm⁻² yr⁻¹ over the same time period. Dry mass accumulation rates increased gradually by ∼2–3-fold at all core sites from settlement to present-day. The oldest ²¹⁰Pb dates in each core extended well before European settlement in the region.

Geochemical proxies

Water content was typically in excess of 68% and decreased downcore (Fig. 3). Organic content increased gradually in all cores following European settlement, most markedly so in the last few decades, although organic content was higher in the Muskeg core than other sites. No clear trend was observed for the proportion of carbonates in the sediments beyond increases in surface deposits and a generally lower carbonate content in Little Traverse Bay. Sediment accumulation rates increased after ca. 1950, especially for carbonates, possibly reflecting increased erosion of shoreline and catchment materials. Accumulation of inorganic materials plateaued since ca. 2000, whereas carbonate accumulation continued to increase at the Little Traverse, Big Narrows, Muskeg, and Sabaskong sites, all locations that are nearer to shore relative to the Big Traverse cores.

Figure 3. Proportions and accumulation rates of sedimentary components in cores collected from Lake of the Woods. Interval dates represent a constant-rate-of-supply age model applied to the ²¹⁰Pb data. Basin-specific line plots are identified by their corresponding color in the heading.

BSi comprised 2–6% dry weight of the lake sediment (Fig. 4). Elevated concentrations were noted in near surface sediments of Big Traverse 3 and 4, and Muskeg bays, whereas Sabaskong, Little Traverse, and Big Narrows had more gradual increases in BSi, mostly since the 1950s. When converted to accumulation rates, the flux of BSi increases toward the top of all cores with particularly marked increases in flux since 1950 in Big Narrows, Little Traverse Bay, and Big Traverse 3. By contrast, Big Traverse 4 and Sabaskong Bay exhibited continually increasing flux of BSi upcore, whereas Muskeg Bay had few long-term changes in BSi flux except in the most recent deposits.

Figure 4. Biogenic silica of 6 Lake of the Woods cores including sediment biogenic silica concentration (SiO₂ wt %; mg/100 mg sediment) and flux (mg/cm² yr). Also included are abundances of diatom valves.

Concentrations of diatom valves generally followed trends for BSi flux (i.e., increasing up-core in recent decades; Fig. 4). Valve concentrations are periodically high in Big Traverse 3 from ca. 1930 through ca. 1990, consistent with higher BSi content in those intervals. Note that assemblages deposited after 1980 are generally represented by smaller diatoms that require less silica per cell. Nonetheless, diatom densities were correlated significantly with sedimentary BSi (Table 1).

Table 1. Correlation coefficients among diatom-related variables measured from LoW sediment cores. No transformations were applied to variable data. ^*** = significant at P < 0.001; ^* = significant at 0.01 < P < 0.05.

	All samples	Since 1980
	(n = 120)	(n = 57)
Diatoxanthin (nmoles/g C) vs. % BSi	0.45^*^*^*	0.31^*
Diatom valves (/g dry wt.) vs. % BSi	0.62^*^*^*	0.34^*
Diatom valves (/g dry wt.) vs. Diatoxanthin (nmoles/g C)	0.34^*^*^*	0.33^*

Fossil algal pigments

Analysis of fossil pigments revealed 2 main stratigraphic patterns in the cores (Fig. 5; profiles of accumulation rates, which show nearly identical trends, are provided in Supplement D). First, many sites, including Big Traverse 3 and 4, Sabaskong, and Big Narrows, show a bimodal pattern of pigment concentrations, with transient peaks around the 1960s–1970s, decreases during the 1980s, followed by increased concentrations during last 15 years. These patterns are reflected in pigment profiles from most algal groups including diatoms (diatoxanthin), cyanobacteria (canthaxanthin, echinone), and general bloom indicators (lutein-zeaxanthin). The second and more recent change is characterized by a marked rise in myxoxanthophyll from colonial cyanophytes beginning around 1980. Timing of this later increase was similar to potential changes in the deepwater environment that favored chlorophyll preservation, as indicated by a slow rise in ratios of chlorophyll a to pheophytin a (Leavitt and Hodgson 2001). Overall, total algal abundance did not vary strongly in any core, as indicated by muted historical changes in the concentration of the chemically stable ubiquitous pigment β-carotene.

Figure 5. Concentration (nmol pigment/g organic C) of fossil pigments and pigment derivatives preserved in 6 cores from Lake of the Woods. Pigments and derivatives are grouped by algae type or indicator value.

Correlation analysis among the diatom pigment diatoxanthin, % BSi, and diatom valve concentrations indicated significant correlations among these variables (Table 1). Performing the same analyses constrained to post-1980 samples resulted in weaker (but still significant) correlations, however, suggesting a weaker modern relationship between BSi and diatom numbers that suggests lower per-cell silica requirements for more recent diatom species.

Siliceous microfossils

In general, the accumulation of diatoms and stomatocysts increased with time, particularly in the most recent decade (Fig. 6). Phytolith accumulation increased in the 4 sites representing smaller embayments (Little Traverse, Big Narrows, Muskeg, Sabaskong). Overall, the ratio of chrysophyte stomatocysts to diatoms declined during the 20^th century, indicating long-term nutrient enrichment (Smol 1985).

Figure 6. Downcore sediment trends for diatoms, chrysophyte stomatocysts, and phytoliths from Lake of the Woods. Data are presented as accumulation rates (/cm²/y) or relative abundance (%) to illustrate trends. The ratio of stomatocysts to diatoms is a unitless ratio. Basin-specific line plots are identified by their corresponding color in the heading.

Depth-constrained cluster analyses of diatom profiles showed remarkably consistent results. In all cores, clustering with broken-stick analysis identified 3 significant zones: the 1800s through the early 20th century, a mid-20^th-century period, and a recent cluster starting around 1980–2000. These most recent clusters represented the increasing relative abundance of taxa characteristic of high nutrient status (e.g., C. dubius, Stephanodiscus sp. 10; Reavie et al. 2014a) as well as species that may be sensitive to climate-induced changes in stratification (e.g., Discostella stelligera; Saros et al. 2012).

In all cores, pre-settlement assemblages were dominated by Aulacoseira islandica (Fig. 7). This taxon is a well-known spring diatom from large, north-temperate, shallow lakes (e.g., Reavie et al. 2014b). Pre-settlement subdominants included Aulacoseira subarctica (another winter–spring diatom) and Aulacoseira granulata, a summer diatom known to prefer higher nutrient conditions (Reavie et al. 2014a).

Figure 7. Stratigraphies of the relative abundance of dominant diatom taxa identified in sediment cores from Lake of the Woods. Interval dates represent a constant-rate-of-supply age model applied to the ²¹⁰Pb data. Taxa shown occurred at a relative abundance of at least 4% in at least one sample in any core. Note that x-axes vary in scale.

Before ca.1980, reorganization of the diatom communities occurred gradually at all core locations. In general, changes began in the early 20th century with a slow replacement of A. islandica by A. subarctica, A. granulata, and Aulacoseira ambigua (Fig. 7). A more substantial reorganization of the diatom assemblages occurred in the last 2–3 decades, with a shift to planktonic taxa known to prefer eutrophic conditions: Stephanodiscus medius, S. minutus, S. parvus, S. sp. 10, and Cyclostephanos dubius (Reavie et al. 2014a, Reavie and Kireta 2015).

NMDS analysis (Fig. 8) demonstrated that most diatom species assemblages exhibited similar trends in development over the last 150 years, from an assemblage dominated by Aulacoseira (upper, right section; e.g., forms of Aulacoseira islandica [AULISLTW, AULISLWA, AULISLWI], Aulacoseira subarctica [AULSUBA]) to one more represented by smaller centric diatoms (lower, left section; e.g., Stephanodiscus sp. 10 [STESP10], Cyclostephanos tholiformis [CYSTHOL]). Post-1980 trajectories of sample scores were less consistent, reflecting a rapid ontogeny of assemblages in recent decades.

Figure 8. Sample and taxa scores from a NMDS analysis of fossil diatom assemblages in the Lake of the Woods cores. Sample scores are traced over time from the oldest sample analyzed for diatoms to 2012 (the sampling year). The year 1980 is marked on each trajectory by a star. See text for abbreviations used for key species, and the full list of taxa and codes is provided in Supplement B.

Spatial variation in historical trends was also evident among cores, reflecting variation in the morphological forms of A. islandica. Discrete forms of A. islandica can be separated by size, ornamentation, and degree of silicification (Cleve-Euler 1951), all reflecting variation in silica supply, growth rates, and seasonal and life cycle responses in this species (Burge and Edlund 2015). A heavily silicified form with a relatively high diameter-to-length ratio (“thick wide,” AULISLWI) dominated in Sabaskong Bay, setting it apart from other cores. Also, the Sabaskong core contained no Aulacoseira distans v. nivalis (AULDINI) in contrast to the other cores. The northern cores (Big Narrows and Little Traverse) were similar to each other and dominated by a heavily silicified form of A. islandica (“thick wall,” AULISLTW) that did not exhibit broader diameters. The western cores (Big Traverse 3 and Muskeg) exhibited intermediate assemblages with a mix of forms.

Patterns in the diatom assemblages (DCA scores; Figs. 2i–j) varied most consistently as a function of minimum annual atmospheric temperature (Table 2; Fig. 2f), with higher modern minimum temperatures correlated with higher DCA axis scores and rapid changes in assemblages after 1980. The recent “rise” in the DCA scores may also be related to lower N:P ratios (Fig. 2c) and the recent decline in nitrate deposition (Table 2; Fig. 2d). The early 20^th century increase in population and agricultural activity (likely surrogates for point sources, for which data were sparse) within the catchment coincides with apparent water quality degradation from the mid-20^th century (Fig. 2b; Hargan et al. 2011). Precipitation generally increased in the region throughout the last century, but the timing of increase was decoupled from the more recent, rapid changes in the diatoms (Table 2; Figs. 2g, i, and j). The day of ice-off has fluctuated from day 98 to day 140 since the early 1980s, having no meaningful relationship with diatom assemblages (Figs. 2g, i, and j). Wind speeds from nearby International Falls declined from a running average of ∼4 m/s in the late 1940s to ∼3.3 m/s in the early 2000s (Fig. 2h), but these changes did not correlate strongly with the rise in diatom DCA scores (Figs. 2i–j).

Table 2. Spearman correlation coefficients of diatom assemblage DCA scores against environmental variables. Parameters included are those with sufficient long-term data for a reasonable comparison and those with potential for direct influence on diatom communities. The number of fossil diatom samples used in each comparison is provided in square brackets. Significance of each correlation (t-test) is shown at P levels of 0.05 (^*), 0.01 (^**) and 0.001 (^***). Units for each parameter are provided in Fig. 2.

Core location	DCA axis	Atm. NH4^+	Atm. NO3^−	Atm. inorganic N	Precipitation	Maximum temperature	Minimum temperature	Mean temperature	Ice-off day	Wind speed
Big Narrows	DCA1	0.62^* [11]	−0.77^*^* [11]	0.19 [11]	0.17 [16]	0.71^*^* [17]	0.64^*^* [17]	0.50^* [17]	−0.32 [10]	−0.37 [10]
	DCA2	0.69^* [11]	−0.74^*^* [11]	0.33 [11]	0.16 [16]	0.71^*^* [17]	0.68^*^* [17]	0.48 [17]	−0.33 [10]	−0.28 [10]
Big Traverse 3	DCA1	0.56 [11]	−0.70^* [11]	0.20 [11]	0.14 [19]	0.42 [19]	0.60^*^* [19]	0.15 [19]	−0.20 [10]	−0.21 [11]
	DCA2	0.56 [11]	−0.61^* [11]	0.24 [11]	0.40 [19]	0.22 [19]	0.58^*^* [19]	0.07 [19]	−0.23 [10]	−0.55 [11]
Big Traverse 4	DCA1	0.68^* [12]	0.01 [12]	0.45 [12]	0.37 [20]	−0.06 [20]	0.21 [20]	−0.35 [20]	0.19 [11]	−0.34 [12]
	DCA2	0.55 [12]	−0.44 [12]	0.25 [12]	0.34 [20]	0.20 [20]	0.74^*^*^* [20]	0.04 [20]	0.30 [11]	−0.41 [12]
Little Traverse	DCA1	0.53 [9]	−0.67^* [9]	0.23 [9]	0.40 [16]	0.47 [16]	0.70^*^* [16]	0.25 [16]	−0.33 [8]	−0.58 [9]
	DCA2	0.63 [9]	−0.68^* [9]	0.26 [9]	0.48 [16]	0.53^* [16]	0.77^*^*^* [16]	0.30 [16]	−0.38 [8]	−0.48 [9]
Muskeg	DCA1	0.85^*^* [8]	−0.65 [8]	0.40 [8]	0.45 [16]	0.14 [16]	0.61^* [16]	0.06 [16]	−0.55 [8]	−0.14 [8]
	DCA2	0.83^* [8]	−0.57 [8]	0.42 [8]	0.55^* [16]	0.12 [16]	0.52^* [16]	−0.01 [16]	−0.57 [8]	0.02 [8]
Sabaskong	DCA1	0.37 [10]	−0.77^*^* [10]	−0.05 [10]	0.09 [17]	0.37 [17]	0.42 [17]	−0.12 [17]	−0.56 [9]	−0.39 [10]
	DCA2	0.30 [10]	−0.77^*^* [10]	−0.08 [10]	0.43 [17]	0.43 [17]	0.80^*^*^* [17]	0.46 [17]	−0.42 [9]	−0.55 [10]
All cores	DCA1	0.49^*^*^* [61]	−0.27^* [61]	0.31^* [61]	0.15 [104]	0.26^*^* [105]	0.37^*^*^* [105]	0.04 [105]	−0.19 [56]	−0.03 [60]
	DCA2	0.49^*^*^* [61]	−0.48^*^*^* [61]	0.19 [61]	0.33^*^*^* [104]	0.34^*^*^* [105]	0.62^*^*^* [105]	0.20^* [105]	−0.32^* [56]	−0.40^*^* [60]

Diatom-inferred total phosphorus (DI-TP) concentration

In all cases,Table 2. modern (model) samples were sufficiently similar to each fossil sample to justify application of the TP model developed for Minnesota lakes. Further, all fossil sample scores were within the 95% distance from the TP axis, indicating good fit to TP for all fossil samples. Evaluation of the fossil assemblages relative to a redundancy analysis (RDA) of the model training set indicates that all fossil diatom profiles except Big Traverse 3 were significantly related to the TP gradient (Table 3), demonstrating that historical changes in diatom assemblages were related to variation in past water column P content. The apparent increase in productivity suggested by pigments that peaked in the 1960s–1970s is not similarly inferred by DI-TP in all of the cores.

Table 3. Correlations from a regression of diatom-inferred total phosphorus against axis 1 sample scores from a principal components analysis of sedimentary diatom assemblages. Asterisks indicate significant correlations (P < 0.05).

Core	R
Big Traverse 3	−0.15
Big Traverse 4	−0.7^*
Muskeg	−0.83^*
Sabaskong	−0.82^*
Little Traverse	−0.89^*
Big Narrows	−0.77^*
Buffalo	−0.91^*

Although the timing and magnitude of diatom-inferred changes in TP varied among cores, the general trend was one of increasing P concentrations starting around the 1930s or 1940s (Fig. 9). Big Traverse 3, Big Traverse 4, and Sabaskong cores indicate an increase in DI-TP, from ∼15 µg/L to as high as 25 µg/L from ca. 1945 until ca. 1970, after which DI-TP declined. The trend in these 3 cores was similar to that of peaks in biomarker pigments for diatoms (diatoxanthin), total cyanophytes (echinone), and colonial cyanobacteria (canthaxanthin). Inferences from Sabaskong and Big Traverse reflect a recovery from eutrophication from the 1970s (∼25 µg/L TP) through ∼2000 (∼14 µg/L TP), with higher P inferred for the most recent few years. Long-term trends of persistent eutrophication were recorded in Muskeg, Little Traverse, and Big Narrows sites, increasing from pre-impact concentrations of ∼10 µg/L to >30 µg/L in recent years. With the exception of the surface interval, DI-TP values appeared to underestimate measured TP concentrations. For example, despite peaks of ∼25 µg/L DI-TP in certain cores in the late 1960s, measurements exceeded 50 µg/L TP during those years. Discussion of possible issues with the Minnesota diatom model are provided in Supplement C.

Figure 9. Diatom-inferred TP concentrations based on sedimentary assemblages from Lake of the Woods. Basin-specific line plots are identified by their corresponding color in the heading. The inset plot expands the DI-TP axis and includes TP monitoring data (source MPCA Citizen Lake Monitoring Program).

Discussion

Based on integrated results from 6 cores in the southern basin of the lake, paleolimnological analyses demonstrated a mid-to-late 20^th century period of enrichment and a period of water quality degradation since ca. 1980, despite a reduction in external nutrient load notably from point sources. In general, phytoplankton abundance increased consistently during the early 20th century, with concomitant and coherent shifts in diatom community composition. Initial eutrophication likely reflects increasing post-settlement nutrient supplies while increases in sediment supplies came from shoreline erosion resulting from damming and land use changes such as logging. The increase in phytolith accumulation in 4 of the core sites nearer to shorelines (Little Traverse, Big Narrows, Muskeg, and possibly Sabaskong) was consistent with increasing erosion and deposition of terrestrial and riparian phytoliths (Piperno 2002). Multiple indicators suggest that the pelagic ecology of the lake has shifted dramatically since ca. 1980. Given the need for aquatic management recommendations, this recent change is a major focus of the following discussion.

Fossil diatoms provided strong evidence of nutrient fertilization of the lake. Before settlement, the dominant, pre-1900 presence of Aulacoseira islandica suggests a naturally mesotrophic condition (Reavie et al. 2014a) as well as a wind-driven physical control on limnology, which is needed to ensure the heavy valves of Aulacoseira remain suspended in the water column. Thereafter, diatom community composition shifted to greater abundance of species with higher TP optima including Cyclostephanos dubius, several small Stephanodiscus species, and Aulacoseira granulata (Reavie et al. 2014a). In concordance, the long-term decline in the ratio of stomatocysts to diatoms also suggested a shift to higher productivity because chrysophytes tend to be relatively more abundant under oligotrophic conditions (Smol 1985). The most recent sediments contain diatom assemblages not previously seen in the lake, representing continued high productivity and a condition reflecting drivers in addition to nutrients. Concerns that diatom preservation may be favoring these smaller, recent taxa relative to more silicified taxa (e.g., Aulacoseira) in upper sediment intervals prompted us to reexamine diatom slides for signs of dissolution or breakage with depth. We found no evidence of preservation issues; in other words, these smaller centric diatoms did not exhibit increased breakage or dissolution (“fading away”) with sediment depth. Hence, we are convinced that recent diatom trends reflect a true ecological shift.

Despite a decrease in nutrient flux to the lake (Hargan et al. 2011) and some evidence of recovery since the mid 20th century, the weight of indicator evidence suggests nutrient levels in the lake have not declined appreciably since ca. 1980. Pigments indicate conclusively a recent increase in cyanophytes, as observed elsewhere in subarctic lakes (Taranu et al. 2015). Although some evidence exists of increased algal abundance in the uppermost intervals, it is not occurring in all locations, and it does not seem that average year-round algal abundance is substantially increasing. Instead, the increase in cyanophytes was accompanied by a transition to a diatom assemblage of smaller, centric species, some of which are known to tolerate high epilimnetic nutrient concentrations (e.g., Stephanodiscus spp.; Reavie et al. 2014a, Supplementary Material), but not necessarily more tolerant than the Aulacoseira taxa that dominated in previous decades. The increase in diatom numbers seemed to include mostly smaller-celled forms and hence does not reflect a substantial increase in diatom biovolume. A decoupling between BSi and diatom numbers is evident from the weaker modern relationship between both diatom valve densities and the diatom pigment diatoxanthin and BSi, possibly a result of the shift to smaller, centric species that require less BSi per cell.

Different proxies revealed mechanisms affecting lake productivity. For example, pigment records of bloom-forming taxa (e.g., lutein-zeaxanthin) and diatoms (e.g., diatoxanthin) suggest 2 periods of elevated phytoplankton abundance since 1900 (1950–1970, post-1980). This pattern documented water quality degradation in the mid 20^th century resulting from excess nutrient flux from municipal and industrial discharges to the Rainy River (Hargan et al. 2011), followed by large reductions in point source nutrient loads beginning with remediation efforts in the 1970s. Since the 1980s, pigments indicated a second period of increased diatom productivity in some core locations, a pattern that generally matches fossil BSi content. Concentrations of pigments from colonial cyanophytes (myxoxanthophyll) became more common, consistent with observations of Microcystis aeruginosa and Aphanizomenon flos-aquae (per Chen et al. 2009), which have been observed in recent algal collections in high abundance (Watson and Kling, 2017).

What is causing the recent change?

Recent changes in Lake of the Woods are consistent with hypotheses that climate warming has altered stratification regimes and, possibly, internal nutrient loading. Warming atmospheric temperatures may be an important driver of recent changes; there was clear correlation between recent shifts in the lake diatom assemblages and atmospheric warming. Although a recent, slight increase in diatoms from the Cyclotella sensu lato complex (Cyclotella, Discostella) occurred in the Big Traverse cores, the increase we observed was far more subtle than that reconstructed in the deeper, northern basins of the lake (Rühland et al. 2010). The stratigraphic shift toward more Cyclotella in recent decades is generally correlated with lake warming (Rühland et al. 2008), but the recent rise in Cyclotella has generally occurred in deeper, dimictic lakes (e.g., Rühland et al. 2008, Saros et al. 2012). Therefore, whether we would expect a similar reorganization in the diatom assemblages in the southern basin of the lake requires better understanding of the effects of climate change on shallow lakes.

Increasing cyanophyte abundance in temperate lakes has been attributed to a decrease in the N:P ratio (Smith 1983), although this has been disputed (Downing et al. 2001). The stoichiometric N:P ratio may be relevant to these recent changes in Lake of the Woods; sampling data from recent decades in many bays (Clark et al. 2014) indicate unchanging or declining trends in N (nitrates and TN) in nearshore locations around the lake. Monitoring data to date (Fig. 2) suggest a shift in nutrient limitation (DIN:TP; Bergström 2010) occurred between the 1960s and 1990s. During the mid 20th century the lake was P limited, and in more recent years, regular conditions of N limitation or N/P co-limitation are likely playing a role in contemporary algal composition. Additional factors determining cyanophyte abundance, such as the important role of Fe (Molot et al. 2014) require further study.

Past activities around the lake resulted in large mass transfers of P (and N) to the lake. Although diffuse agricultural sources exist, the greatest source has been industrial and urban discharges, particularly in the Rainy River catchment (Hargan et al. 2011). From this study it is obvious that, despite the known reduction in nutrient flux to the southern basin from known point sources, nutrient enrichment and associated algal blooms continue to be a problem. Paleolimnological records of algae, sediment accumulation, and observations of extensive cyanobacterial blooms confirm these issues. Increased seasonal, episodic releases of P from lake sediments may be masking or buffering the positive effects of contemporary conservation measures (Sharpley et al. 2013).

Although the southern basin of the lake is a large, shallow system, comparable issues with remediation have occurred elsewhere. Loch Leven is a shallow lake that is relatively large by UK standards, and it has a long history of eutrophication problems. During the late 1970s to early 1990s, external P inputs were reduced by ∼60% by reducing industrial inputs, upgrading sewage treatment works, and implementing better farming practices (May et al. 2012). P concentrations in the lake dropped ∼75% in winter and 60% in spring, yet summer concentrations of P increased (Sharpley et al. 2013), suggesting a switch in dominance from external to internal P loading in the lake. Similarly, in Sweden's largest lake, Lake Vänern, phytoplankton biomass has increased since the 1980s while water P and N concentrations have decreased (Weyhenmeyer and Broberg 2014). It was determined that increasing lake temperatures were more important for phytoplankton growth than nutrient inputs from the catchment. Lake Vänern and Lake of the Woods are comparable in size and so may experience similar responses to atmospheric warming. Given the correlation between climatic and paleophycological data, Lake of the Woods may be experiencing analogous changes. Higher temperatures (Fig. 2) and the increasing ice-free season in the northern basins (Clark et al. 2014) are likely having some effect on ecological conditions. Although it is unknown whether lake stratification and low hypolimnetic oxygen (Heathcote et al. 2016) are new phenomena for the southern basin, it is conceivable that stratification has become more common with warmer summers, therefore increasing the frequency of low oxygen and episodic P resupply from the sediments. Increasing stratification is further suggested by the loss of heavily silicified diatoms that require isothermal mixing conditions (Aulacoseira; Lindenschmidt and Chorus 1998) in favor of lightly silicified centric species that would be better supported in an epilimnion.

Conclusions and recommendations

Even with reduced nutrient loads to a managed system like Lake of the Woods, our multi-proxy paleolimnological study indicates that confounding anthropogenic and environmental factors can contribute to continued water quality problems. Despite remediation attempts, the stoichiometry of nutrients in the pelagic regions of the lake are changing, likely caused by confounding climate effects and other factors facilitating the periodic resupply of nutrients from the sediments. To provide refined management recommendations and predict future water quality, a detailed historical nutrient budget characterizing P dynamics is needed (Edlund et al., 2017; James et al., 2017). Further, continued monitoring of water column profiles for temperature and oxygen and more detailed collections of nutrient data (in-lake concentrations and point and diffuse sources) will be valuable to refine estimates of lake condition.

Funding

Funding for this work was provided by grants from the MPCA under the Clean Water, Land and Legacy Amendment (CWLLA), and Clean Water Partnership (CWP) program.

Acknowledgments

Steve Heiskary (MPCA), Mark Tomasek (MPCA), and Joy Ramstack (SCWRS) provided long-term field support and helped compile the diatom model dataset. Nolan Baratono (MPCA), Randy Beebe (Wolfshead Logistics), and the Sportsman's Lodge provided field support. Jill Coleman Wasik, Erin Mortenson, Alaina Fedie, and Erin Mittag from SCWRS provided analytical support. We thank William Herb for the use of the wind data from International Falls, MN.