Lake Petit Saint-François (PSF; latitude 45.54, longitude −72.04) and Roxton Pond (latitude 45.47, longitude −72.66) are small (surface area < 200 ha), shallow (mean depth ≤ 3.1 m) lakes located in the Eastern Townships region of Quebec, Canada. The Eastern Townships are an important tourism and recreation area located approximately 150 km southeast of Montreal, south of the St. Lawrence River, and north of the U.S. states of Vermont, New Hampshire, and Maine. Although there has been considerable agricultural abandonment and forest regrowth in parts of this region during the twentieth century, agriculture has been a major land use over the last ~100 years in the catchments of our study lakes (MacDonald and Bennett, 2009). The study lakes cover a eutrophic to hypereutrophic range of contemporary mean summer TP and phytoplankton concentrations, measured as chlorophyll a (Chl a). PSF Lake is presently hypereutrophic (2009–2015 mean TP = 216 μg/L, and Chl a = 71 μg/L; http://www.mddelcc.gouv.qc.ca/eau/rsvl/index.asp) and contains almost no submerged macrophytes, whereas Roxton Pond is eutrophic (2004–2015 mean TP = 39 μg/L, and Chl a = 28 μg/L; http://www.mddelcc.gouv.qc.ca/eau/rsvl/index.asp) and still contained dense and widespread submerged macrophyte beds at the time of sampling. Cyanobacteria are an important part of the contemporary phytoplankton community in both of these lakes and elevated levels of microcystin toxin have been reported (Giani et al., 2005; Lévesque et al., 2014; Monchamp et al., 2014). Both study lakes have had their water levels artificially raised by small dams constructed in the early to mid-twentieth century. The height of the dams for lakes PSF and Roxton are 1.6 m (constructed in 1966) and 4.7 m (constructed in 1920), respectively (Quebec Hydrology Expertise Centre; http://www.cehq.gouv.qc.ca/index.asp).

Sediment cores for diatom, zooplankton, and pigment analyses were collected from near the deepest region of each lake in the summer of 2008 (1.9 m in PSF and 5.5 m in Roxton). A 78 cm long sediment core was obtained from PSF using a push corer that was then extruded at 1 cm intervals in the field to a depth (~15 cm) where the sediment was sufficiently consolidated to allow shipment back to the laboratory. In the laboratory, the PSF sediment core was split longitudinally and half the core was sub-sampled at 1 cm intervals for analysis, the other half was archived. A 39 cm long sediment core was obtained from Roxton Pond using a Glew core (Glew, 1989), which was extruded in the field at 0.5 cm intervals. The sub-sampled sediment was kept at 4°C in the dark until processed. The chronologies for the sediment cores were determined through radiometric ²¹⁰Pb and ¹³⁷Cs dating at McGill University using a Canberra gamma spectrometer with the constant rate of supply model (Appleby, 2001).

Diatom slide preparation and enumeration followed standard procedures (Battarbee et al., 2001). Briefly, samples were first treated with 10% HCl to remove any carbonates, rinsed several times with deionized water and then treated with boiling H₂O₂ to remove organics. Samples were subsequently rinsed with deionized water seven times, plated onto coverslips, and mounted with Zrax. For each sample, a minimum of 400 diatom valves were counted and identified following Cumming et al. (1995), Fallu et al. (2000), and Krammer and Lange-Bertalot (2004).

To characterize the littoral zooplankton assemblage, we quantified the proportion of chydorid subfossils as a function of all cladoceran carapace subfossils preserved in each sediment interval following the methods described in Korhola and Rautio (2001). The community composition of cladocera in the study region is driven largely by nutrient concentration and macrophyte abundance (Velghe et al., 2012). For PSF Lake, the chydorid data were previously published in study of a shorter sediment core (Richard Albert et al., 2010). We matched up the shorter sediment record from PSF Lake to data from the longer core using the ²¹⁰Pb-derived age models from each core.

We conducted pigment analyses to quantify the concentration of the cyanobacteria associated pigment echinenone following standard methods in paleolimnology (Leavitt and Hodgson, 2001). Briefly, to extract pigments from the sediments, we placed approximately 0.1 g of freshly freeze-dried sediment into 100% acetone, purged this sample with Argon, and stored it at −20°C for 24 h. The extract was subsequently filtered through a PTFE 0.2 μm filter and then injected into a Waters, high performance liquid chromatography (HPLC) system. Pigments were analyzed following the protocol of Zapata et al. (2000) and quantified using external standards purchased from DHI (Denmark).

As a proxy for local external stress on the lake ecosystems due to agricultural activities, we calculated phosphorus balances using county-level agricultural census data (counties are called “census division” in more recent censuses), which are available at 10-year intervals from 1871 to 2001 (e.g., decadal resolution of land-use change). Our phosphorus balances (detailed in MacDonald and Bennett, 2009) account for the net-phosphorus inputs to cropland soils via the application of livestock manure and mineral fertilizers after subtracting the phosphorus removed from cropland soils by harvested crops annually. These phosphorus balances provide an areal estimate of the magnitude of phosphorus accumulation in cropland soils (kg ha⁻¹ year⁻¹) and provide a long-term record of changing land use and management practices that could alter catchment soil phosphorus.

We calculated an approximate “phosphorus load index” (P load index) in g m⁻³ year⁻¹ by multiplying the areal phosphorus balances by the catchment area and scaling this measure by lake volume. Lake volume was estimated from bathymetric maps and we calculated the historical phosphorus load index using both estimates of modern lake volume and pre-dam estimates of lake volume as appropriate. The phosphorus load index is only indicative of changes in agricultural phosphorus management, and does not account for other anthropogenic sources of phosphorus (e.g., residential septic tanks systems); however, our analysis of aerial photographs from approximately 1945 to present indicates that agriculture has been the predominant land use type in the catchments of the study lakes. In addition, the counties are substantially larger than the catchments of the study lakes, but our earlier studies have shown strong relationships between the historical census land use records, catchment soil phosphorus concentrations, and lake eutrophication trends (e.g., MacDonald and Bennett, 2009; Richard Albert et al., 2010; Shaw Chraibi et al., 2011), suggesting that census data can provide a reliable historical estimate of land-use intensity and phosphorus loading in this region.

In addition to excessive nutrient loading, climate warming has also been shown to be an important predictor of cyanobacteria dynamics in lakes (Taranu et al., 2015). To examine if changes in mean annual temperature over the last century were related to changes in cyanobacteria concentrations in our study lakes we used the adjusted and homogenized climate data from the Brome, Quebec, weather station (latitude 45.18; longitude −72.57; altitude 206 m; sampling years 1891–2013; Vincent et al., 2012). The Brome weather station was the closest weather station to our two study sites (approximately 75 km from PSF and 40 km from Roxton Pond) where adjusted and homogenized climate data was available over a similar time period to our phosphorus load index estimates (~last 100 years).

To examine if PSF and Roxton have undergone significant shifts in their algal assemblage or concentrations we employed clustering analyses constrained by time. Stratigraphically constrained cluster analysis was applied using CONISS (Grimm, 1987) to examine clustering in the diatom assemblages using the rioja package (Juggins, 2016) in R (R Development Core Team, 2016). For the higher temporal resolution cyanobacterial pigment data, we tested whether PSF and Roxton Pond had statistically significant groupings in cyanobacteria concentrations over the sediment record using a regression tree approach constrained by time (De'ath, 2002) using the rpart package and MVPARTwrap wrapper function in the R statistical program (Ouellette and Legendre, 2013; Therneau et al., 2014; R Development Core Team, 2016). The cross-validation procedure available in MVPARTwrap helped determine the best number of statistically significant groupings (Legendre and Legendre, 2012). As recommended by De'ath (2002), the model containing the fewest number of groups within one standard error of the optimal model was selected. We then used an ANOVA and a post-hoc Tukey test to examine if the different groups identified in the temporally-constrained regression tree analyses had significantly different mean echinenone concentrations. If the regression tree groups did not have statistically different mean echinenone concentrations, we considered them part of the same cyanobacterial group of the lake whereas if the mean echinenone was significantly different we considered this period as a different cyanobacterial group for these shallow lakes.

Once historical cyanobacterial groups were identified based on the regression tree, we quantified the trends in echinenone concentration across the transition corresponding to the timing of the catchment P load index (from 1871 to 2001) and the observational temperature record by fitting linear, unimodal, threshold and sigmoidal models to both the environment (catchment P load index and mean annual temperature) and response (echinenone concentration) series using generalized nonlinear least squares regression (Seddon et al., 2014). The best model for each variable was then identified using the Akaike Information Criterion (AIC). To evaluate the relative importance of intrinsic vs. extrinsic dynamics over this transition period, we compared the best model for the response and environmental driver variables along a timeseries. If the two models corresponded to each other (e.g., step-change in response and step-change in driver) an extrinsic regime shift was inferred, whereas if the models did not correspond (e.g., step-change in response and linear increase in driver) an underlying intrinsic regime shift was presumed.

The ²¹⁰Pb inventories for each sediment core showed that our records spanned at least the last 150 years (Figure 1 and Figure S1). Diatom analysis indicated that PSF Lake was nutrient rich before the industrial era (i.e., pre-1850) with Aulacoseira ambigua and small benthic Fragilaria spp. (sensu lato) dominating the record (Figure 1A). In the mid-nineteenth century, however, there was a pronounced decrease in the proportion of planktonic diatoms and a shift toward epiphytic taxa that are associated with lower phosphorus conditions (e.g., Achnanthidium minutissimum; Hall and Smol, 1992; Vermaire et al., 2011). This record suggests that PSF transitioned from a historically eutrophic turbid-water state to a clear-water, macrophyte-dominated state. During this clear-water, macrophyte-dominated phase, we noted a corresponding decrease in cyanobacteria pigments, suggesting that this was an ecosystem-wide event. Around 1910, however, the diatom record transitioned back to taxa associated with greater nutrient concentrations, which was also reflected in the loss of macrophyte-associated chydorid taxa and an increase in cyanobacteria pigments.

In contrast, Roxton Pond was in a lower-nutrient state before the industrial era (pre-1850). The dominant planktonic (e.g., Discostella stelligera, Tabellaria flocculosa) and benthic diatom taxa (e.g., A. minutissimum) observed in that period are known to be associated with low phosphorus conditions (Hall and Smol, 1992; Dixit et al., 1999; Figure 1B). During this time, macrophyte-associated chydorid taxa were relatively abundant in the bottom half of the core and cyanobacteria concentrations were low, suggesting a clear-water, macrophyte-dominated state. Following an increase in historical catchment phosphorus loading around 1950, there was a rise in the relative abundance of planktonic diatoms known to be indicative of higher phosphorus conditions (e.g., Fragilaria crotonensis, Aulacoseira ambigua; Hall and Smol, 1992; Dixit et al., 1999) and decline in the relative abundance of chydorid taxa.

Regression tree analysis of the echinenone concentration constrained by time suggests four groups for the PSF record and seven for the Roxton Pond record (Figure 2). The trends in echinenone concentration suggest that both PSF and Roxton Pond have undergone abrupt shifts between higher and lower cyanobacteria concentrations (Figure 3). Although there is agreement between the general trends in diatom community assemblage and echinenone concentration, the record of echinenone concentration suggests that more shifts have occurred possibly due to the echinenone being analyzed at greater resolution than the diatom samples. Although the regression tree analyses suggest four and seven echinenone groupings for PSF and Roxton Pond, respectively, the ANOVA analyses with subsequent Tukey tests indicate that the lakes have alternated around three (in the case of PSF) or four (in the case of Roxton Pond) cyanobacteria states over their record (Figure 3).

In PSF, four shifts were detected in cyanobacteria concentration. Between ca. 1600 and 1775 PSF was in a very high cyanobacteria regime before transitioning to a moderate cyanobacteria regime between 1775 and 1835 (Figure 3). PSF then underwent a rapid decline in cyanobacteria between 1835 and 1940 before transitioning back to the high cyanobacteria regime between 1940 and 2008. The first three shifts pre-dated our record of environmental time series but the final shift corresponds to a rapid increase in phosphorus loading to the catchment due to agricultural activities.

Similar to PSF, Roxton Pond showed multiple transitions between cyanobacteria groups but followed a more expected pattern transitioning from low cyanobacteria concentrations in the past to higher cyanobacteria concentrations in the latter part of the record. Overall, four distinct cyanobacteria groups were identified in Roxton Pond over the last two centuries (Figure 3). Groups with cyanobacteria minima were observed between the early 1800s to 1875 and from 1935 to 1955 whereas groups with moderate cyanobacteria occurred between 1875 to 1935 and from 1955 to 1995. In ca. 1995 Roxton Pond transitioned to groupings with greater concentrations of cyanobacteria reaching a cyanobacteria maxima group between 1995 and 2005, before declining slightly to more moderate cyanobacteria concentrations.

During the transitions in both PSF (post-1900) and Roxton (post-1870) we observed a nonlinear increase in both cyanobacteria concentrations (echinenone pigment) and the P load index. In PSF Lake the cyanobacteria trend was best described by a step function model (Figure 4), whereas changes in the P load index through time were best described by the segmented model (Figure 5). In Roxton Pond, a segmented linear model best described changes in cyanobacteria concentrations through time with an abrupt increase in the latter part of the record (Figure 4). On the other hand, the historical P load index for the Roxton Pond catchment was best described by a step function model (Figure 5). Trends in the mean annual temperature data from the Brome weather station were best described by a segmented model and suggests more rapid warming in the region after ~1980 (Figure 6).

The abrupt, nonlinear increase in cyanobacteria in PSF and Roxton Pond occurred either during a time of linear increase in the environmental drivers in the case of PSF or approximately 10–20 years after abrupt transitions in the environmental drivers in the case of Roxton Pond (Figure 7). In fact, the rapid nonlinear increase in echinenone pigment occurred at Roxton Pond during a time of stabilizing P load index and linear climate warming trend (Figure 7). These trends suggest that the threshold response of cyanobacteria concentrations in both lakes was driven by intrinsic factors within the ecosystem.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.