Study area
Our study lakes (n = 27) in the Mixedwood Plains ecozone are distributed across northeastern Ontario and southern Quebec (Figure 1). Lake sampling was conducted as part of the LakePulse program during the summer of 2017, following protocols described in Huot et al. (2019) and detailed in NSERC Lake Pulse Network (2021). For each lake, detailed physiographic information (geographic location and lake morphology) as well as water chemistry data were collected and are presented in Table 1.
Table 1
Geographic and water-chemistry data of the 27 Mixedwood Plains study lakes
Lake Name
|
Lake ID
|
Latitude
|
Longitude
|
Area
(ha)
|
Size
(class)
|
Depth (m)
|
Stratification (class)
|
Human Impact (class)
|
Secchi
(m)
|
DOC
(μg L–1)
|
Chl a
(μg L–1)
|
Specific conductance (µS cm–1)
|
TN
(μg L–1)
|
TP
(μg L–1)
|
Lac Saint-Augustin
|
08-097
|
46.749625
|
-71.392292
|
69.6
|
Medium
|
5.5
|
Mixed
|
High
|
0.9
|
5.7
|
82.5
|
631.5
|
620.1
|
194.6
|
Lac des Chicots
|
08-120
|
46.796216
|
-72.523542
|
71.1
|
Medium
|
21.0
|
Stratified
|
High
|
1.0
|
5.7
|
35.1
|
97.5
|
406.8
|
33.2
|
Upper Rock Lake
|
08-148
|
44.494194
|
-76.413815
|
84.0
|
Medium
|
43.0
|
Stratified
|
Low
|
4.7
|
8.3
|
3.6
|
154.2
|
70.9
|
11.0
|
Devil Lake
|
08-151
|
44.577578
|
-76.440764
|
1100.5
|
Large
|
41.0
|
Stratified
|
Low
|
5.5
|
5.1
|
0.9
|
169.5
|
60.2
|
9.4
|
Christie Lake
|
08-152
|
44.603507
|
-76.467922
|
32.6
|
Small
|
6.0
|
Stratified
|
Low
|
4.1
|
7.7
|
4.2
|
62.3
|
66.0
|
11.8
|
Lac Green
|
08-153
|
45.651003
|
-76.495042
|
55.4
|
Medium
|
11.5
|
Stratified
|
High
|
4.4
|
3.5
|
0.9
|
271.5
|
152.9
|
10.9
|
Lacs Healy
|
08-154
|
45.844902
|
-76.653917
|
12.3
|
Small
|
21.0
|
Stratified
|
Low
|
3.5
|
6.3
|
3.7
|
62.5
|
188.7
|
13.5
|
Stoco Lake
|
08-160
|
44.473281
|
-77.289620
|
575.1
|
Large
|
13.5
|
Stratified
|
Low
|
1.8
|
11.2
|
2.7
|
198.5
|
193.1
|
24.2
|
Seymour Lake
|
08-163
|
44.380086
|
-77.814912
|
751.5
|
Large
|
7.0
|
Stratified
|
High
|
1.6
|
6.7
|
7.0
|
250.7
|
222.6
|
31.6
|
Rice Lake
|
08-166
|
44.177900
|
-78.179116
|
9966.2
|
Large
|
6.0
|
Mixed
|
High
|
1.3
|
5.9
|
11.6
|
236.6
|
259.1
|
6.6
|
Head Lake
|
08-170
|
44.740317
|
-78.912505
|
979.1
|
Large
|
5.1
|
Mixed
|
Low
|
2.4
|
4.6
|
4.9
|
140.2
|
177.4
|
11.1
|
Wilcox Lake
|
08-179
|
43.949049
|
-79.436034
|
58.4
|
Medium
|
17.0
|
Stratified
|
High
|
1.7
|
7.2
|
5.3
|
759.3
|
339.8
|
18.4
|
Bass Lake
|
08-181
|
44.603455
|
-79.506646
|
571.4
|
Large
|
9.3
|
Mixed
|
High
|
2.5
|
5.1
|
4.6
|
283.3
|
219.9
|
13.8
|
Little Lake
|
08-182
|
44.426424
|
-79.671388
|
235.4
|
Medium
|
1.7
|
Mixed
|
High
|
0.8
|
9.6
|
16.6
|
477.7
|
455.3
|
29.0
|
Orr Lake
|
08-185
|
44.608247
|
-79.802979
|
320.3
|
Medium
|
2.2
|
Mixed
|
High
|
1.8
|
5.5
|
4.0
|
271.5
|
273.1
|
16.4
|
Schmidt Lake
|
08-189
|
44.167604
|
-81.311053
|
25.1
|
Small
|
2.0
|
Mixed
|
High
|
0.4
|
12.4
|
11.3
|
271.5
|
409.9
|
30.1
|
Pinehurst Lake
|
08-190
|
43.269235
|
-80.390052
|
10.9
|
Small
|
7.5
|
Mixed
|
High
|
2.3
|
8.4
|
6.9
|
294.1
|
287.8
|
18.0
|
Alder Lake
|
08-191
|
43.353353
|
-80.537259
|
16.1
|
Small
|
1.5
|
Mixed
|
High
|
0.5
|
4.4
|
51.3
|
271.5
|
1724.4
|
56.4
|
Dankert Lake
|
08-193
|
44.205289
|
-81.052711
|
23.3
|
Small
|
8.0
|
Stratified
|
High
|
2.4
|
8.9
|
4.0
|
203.9
|
328.8
|
22.0
|
Copps Lake
|
08-197
|
44.258274
|
-80.965301
|
11.4
|
Small
|
10.1
|
Stratified
|
High
|
2.4
|
4.3
|
0.9
|
376.8
|
138.5
|
8.6
|
Gillies Lake
|
08-205
|
45.204954
|
-81.326861
|
224.8
|
Medium
|
31.0
|
Stratified
|
Low
|
2.4
|
6.2
|
1.0
|
283.3
|
138.0
|
9.9
|
Clam Lake
|
08-208
|
44.072133
|
-81.413948
|
40.7
|
Small
|
14.0
|
Stratified
|
High
|
2.2
|
12.4
|
8.3
|
459.0
|
343.1
|
11.9
|
Cameron Lake
|
08-209
|
45.216979
|
-81.557960
|
164.4
|
Medium
|
13.0
|
Stratified
|
Low
|
3.2
|
7.4
|
2.3
|
304.9
|
162.3
|
11.2
|
Sucker Lake
|
08-211
|
45.723555
|
-81.878484
|
228.2
|
Medium
|
3.2
|
Mixed
|
Low
|
2.7
|
6.1
|
5.3
|
184.6
|
237.3
|
12.3
|
Bass Lake
|
08-213
|
45.889627
|
-81.937381
|
273.9
|
Medium
|
11.0
|
Stratified
|
Low
|
4.0
|
6.5
|
4.2
|
301.8
|
196.3
|
25.5
|
NA
|
08-214
|
44.007233
|
-78.029579
|
11.5
|
Small
|
3.0
|
Mixed
|
High
|
1.5
|
10.3
|
5.0
|
277.4
|
393.4
|
15.7
|
NA
|
08-215
|
42.874465
|
-82.183136
|
12.9
|
Small
|
1.3
|
Mixed
|
High
|
0.1
|
8.8
|
81.0
|
271.5
|
684.0
|
148.3
|
Human impact (HI) scores were assigned to each study lake based on weighted differences in the land-use composition of the watershed, ranging from 0 (no disturbance) to 1 (high degree of disturbance; Huot et al. 2019). The lakes were split into low (HI index < 0.15) and high human impact (HI index ≥ 0.15; max score is 1) classes; the threshold in HI scores was calculated on the full Lake Pulse dataset using a univariate regression tree, with the water chemistry data as response variables and HI as the predictor (Griffiths et al. 2022). The lakes were also split into three size classes: small (< 0.5 ha), medium (≥ 0.5 and < 5 ha), and large (≥ 5 ha). Thermal stratification class was assigned based on the stratification status of the lake at the time of sampling. If the lake thermal profile was not available (for 5 of 27 lakes), the stratification status was assigned based upon a univariate conditional inference tree of stratification class as a function of maximum lake depth. This was run on the entire set of LakePulse lakes with thermal profiles. Lakes without thermal profiles but that were > 5.7 m deep were considered as likely stratified (Griffiths et al. 2021).
Besides, historical (1841 CE – 2017 CE) monthly records of air temperature were obtained for each of the 27 study sites from Environment and Climate Change Canada using the R rclimateca package (Dunnington, 2018). Only climate stations with at least a complete year of monthly data were selected to create the mean annual temperature estimates. Rather than assigning air temperature and precipitation estimates based on the nearest station (which can be hundreds of km away), estimates from stations within 75 km of the study site were spatially interpolated, forming rings of raster values (estimated temperature) around the input station. This step ensured that the estimated temperatures for any site were reflective of station input data.
Case study: Lac des Chicots
Lac des Chicots was chosen as our case study lake to conduct detailed diatom and geochemical analyses due to regional conservation concerns. The lake has been subjected to a series of environmental issues over the past few decades, including cyanobacterial blooms and a biological invasion of Myriophyllum spicatum (“Municipalité de Sainte-Thècle” n.d.). Additionally, the lake likely experienced natural events that pre-date the European expansion into this region. Several parts of the lake were littered with submerged tree trunks (resulting in the translated name “Stubs Lake”), likely due to a potential subsidence or flooding of the surrounding terrain prior to European settlement (Veillette 1973).
Lac des Chicots (46.8017° N, 72.5175° W), located at 143 m a.s.l. in the municipality of Sainte-Thècle (Québec, Canada), has a surface area of 71.1 ha and a maximum measured depth of 21 m. The lake measures 2 km in length and 0.8 km in width, with an elongated form running along a north-east/south-west axis. It can be subdivided into three sections, separated by narrows, with several embayments. The lake is supplied by four inlets, including the discharges of Lac Auguste-Leblanc, Lac à la Peinture, and Lac Rose, as well as Vandal Stream, and it has one outlet to the Rivière des Envies.
Lac des Chicots was classified as mesotrophic to eutrophic on the four occasions of summer sampling from 2008 to 2017 (Table S2; Tremblay et al. 2014; Huot et al. 2019; “Gouvernement du Québec” n.d.). It has a shallow Secchi depth (1.5 ± 0.4 m) and high chlorophyll a (chl a) content (14.3 ± 6.6 μg L–1) in comparison to other Mixedwood Plains lakes (Table 1), suggesting a high primary productivity. The lake was hypoxic or anoxic below 4 m depth in the summer of 2008 (G Cabana pers. commun.). Land-use analysis (Figure 2) of the watershed of Lac des Chicots using ArcGIS indicates that around 41% of its watershed is covered by forests, with another 41% used for agricultural production. Urban landscapes account for ~10% of the watershed area, with 7.9% as urban land use (structures) and 3.5% as roads.
History of local settlement in Sainte-Thècle
At the beginning of European settlement and development in the Sainte-Thècle region, Lac des Chicots was of considerable importance for local transportation, fishing, drinking water, and agriculture. In the 1860s, the first permanent settlements were founded around the lake (“Municipalité de Sainte-Thècle” n.d.). In the 1870s, land clearing and logging took place, and public buildings (a parish and a school) were built to serve the growing population. Meanwhile, infrastructure was installed, with roads constructed in the 1870s and railways in the 1880s. By the end of the 19th century, more than 1,000 inhabitants lived in this region. Developments and population growth further accelerated in the 20th century in the town of Sainte-Thècle. The local population reached almost three thousand residents by 1986 (Statistics Canada 1986), while dropped to ~2,500 starting from the 21st century (Statistics Canada 2006).
Coring and dating
Sediment cores were extracted from approximately the deepest point of each of the 27 lakes in summer of 2017. To retrieve cores that reached pre-industrial times, we aimed to collect cores that were at least 32 cm in length when possible, a length targeted based on a survey published sedimentation rates across Canada (Griffiths et al. 2021). To depict the general trend of ecological change that had occurred in the Mixedwood Plains ecozone, we used the “top-bottom” approach (Smol 2008). Specifically, we collected the surface sediment (0-1 cm interval) and a sediment interval from the bottom of the core (X-4 to X-3 cm interval, where X is the total length of the core) from each lake. Indicators preserved in the surface samples represent modern conditions (the “tops”), and those deposited in the bottom samples reveal historical conditions (ideally before the industrialization of North America; the “bottoms”). The bottom 2 cm of each sediment core was discarded due to potential smearing by the core plug. All cores for the top-bottom analyses were extruded in the field using a vertical extruder.
The full sediment core collected at Lac des Chicots was stabilized onsite with Zorbitrol (Topkins et al. 2008) for further analyses at research facilities in Québec and Ontario. This core was scanned using Siemens SOMATOM Definition AS+ 128 computed-tomography (CT) core scanner at the Institut National de la Recherche Scientifique (INRS) in Québec City. The full sediment core was then sent to the University of Laval and split longitudinally in two. The working half was scanned on the non-destructive Cox Analytics micro-X-ray Fluorescence (µ-XRF) ITRAX core scanner located at INRS and subsampled at 1-cm intervals. Subsamples were subsequently freeze-dried for dating and diatom analyses.
210Pb analyses using gamma spectrometry was completed at the Paleoecological Environmental Assessment and Research Laboratory (PEARL) at Queen’s University. 210Pb activities were measured for both the top (0-1 cm) and bottom (X-3 to X-2 cm interval, a slightly lower interval than used for diatoms due to the sharing of sample across multiple researchers and projects) samples for the 27 study lakes, while the Lac des Chicots full core was dated using a Constant Rate of Supply Model (Appleby 2001). One of the regional samples (lake 08-163) was excluded from the top-bottom analyses due to high unsupported 210Pb levels in the bottom sample (greater 210Pb activities than 214Bi), suggesting that the bottom sample may not have reached the pre-industrial age (all the other study lakes met this requirement; Griffiths et al. 2022). The surface sample of lake 08-163 was used in ordination and other analyses requiring just the top samples. The dating model of our case study lake is provided as Figure S1, and its decay curve is within our expectations (i.e., unsupported 210Pb generally decreased exponentially with depth).
Diatom processing
The top and bottom sediment intervals of all study lakes were processed for diatom analyses. For our case study lake, Lac des Chicots, 23 intervals spanning the full length of the core were examined. Diatoms were prepared for microscopic analysis following standard techniques described in Battarbee et al. (2002b). In brief, dried sediment subsamples of 0.01- 0.02 g were taken from each interval. A 10% hydrochloric acid solution was used to dissolve carbonates, and the process of dilution, settling, and aspiration was repeated at least five times until the samples reached a neutral pH. Then, a 30% hydrogen peroxide solution was added to digest organic matter. Samples were put into a water bath and heated at around 80°C for 8 hours for three consecutive days, and the solution was diluted again until they reached a neutral pH. The cleaned diatom suspensions were pipetted onto coverslips and were mounted in Zrax onto glass microscope slides. For each sample, a minimum of 500 diatom valves was enumerated with a Leica DM 2300 microscope (1000x magnification). Diatoms were counted to species, and variety level when possible, following Kramer and Lange-Bertalot (1986; 1988; 1991a; 1991b) and Lavoie et al. (2008).
Numerical analysis
All numerical analyses were conducted in the R platform (R Core Team 2021). First, to visualize the relationships between environmental factors and diatom assemblages, we conducted a Principal Component Analysis (PCA) on the modern diatom species data using the vegan package (Oksanen et al. 2020), and passively plotted environmental factors on the ordination plot. We divided the environmental factors into two groups: (1) physiographic and water quality variables, and (2) human-impact factors. Physiographic variables include geographic location and morphological information, while water quality variables describe the physical (e.g., Secchi) and chemical variables (e.g., nutrients) that may affect lake conditions for diatom communities. Human-impact factors include human impact (HI) index, population in the watershed, and land-use variables (e.g., percent agricultural land use in the catchment). To meet an assumption of the PCA (normal distributions of environmental variables), we tested the normality of the distribution for all environmental variables and transformed those that were not normally distributed, using Box-Cox, logarithmic, or square-root transformations (the type of transformation and transformation coefficient for each parameter are provided in Table S1). When necessary, variables with zero values were offset to positive values by adding 0.00001.
We investigated the positioning of the major diatom functional groups along the first two PCA axes by plotting them passively on the ordination. The major functional groups were divided based on their TP optima into oligotrophic (TP optima < 12 μg L-1; Tremblay et al., 2014) or mesotrophic/eutrophic taxa (TP optima ≥ 12 μg L-1; Tremblay et al., 2014), as well as their habitats into benthic, benthic-planktic (primarily benthic species that sometimes can be entrained in the plankton such as benthic fragilarioids), planktic, or tychoplanktic taxa (with the habitat assignations based mainly on descriptions in Round et al. 1990). Additionally, the dominant species, defined as those with a strong influence on the ordination (having a loading ≥ 0.45 on either PC axis 1 or axis 2), were also passively plotted on the ordination to examine the major drivers of the principal components of variation. We applied Hellinger transformation to all diatom data using the vegan package (Oksanen et al. 2020).
Then, to investigate how the diatom assemblages of Mixedwood Plains lakes have changed over time, we plotted the site scores of the bottom samples (pre-industrial conditions) passively on the ordination using a weighted-averaging model. For each top-bottom pair of diatom assemblages we calculated the Bray-Curtis (BC) dissimilarity index to quantify the magnitude of temporal beta-diversity. To assess how these sites have shifted between their pre-industrial and modern conditions with respect to the two major axes of variation in the diatom assemblages across the landscape (PC axis 1 and 2), we subtracted the site score of each bottom sample on the PCA ordination from its respective top and plotted their differences onto a biplot of the changes across PC axis 1 and PC axis 2.
Lastly, we explored the distributions of a series of biotic and abiotic variables, including changes in diatom functional groups, Bray-Curtis dissimilarity index values, and water chemistry parameters (e.g., nutrients, ions), across three categorical variables: human impact class (high vs. low), stratification class (mixed vs. stratified), and lake surface area size class (large, medium, or small). Differences of these biotic and abiotic variables between categorical groups were tested for statistical significances (P-values) using Wilcoxon signed-rank test or Kruskal-Wallis one-way analysis of variance.
Case study: Lac des Chicots
To track the limnological history of Lac des Chicots and to interpret its changes in the context of the Mixedwood Plains ecozone, we passively plotted its trajectory on the PCA ordination. Then, to draw a more complete inference of its past conditions, we developed a diatom stratigraphic diagram, calculated the planktic to benthic diatom ratio and Hill’s N2 diversity value, as well as analyzed other biological indicators (chrysophyte scales and cysts, sedimentary chlorophyll a). We also investigated several geochemical indicators, including sediment density, calcium (Ca), titanium (Ti), total organic carbon (TOC), and total nitrogen (TN).
The diatom stratigraphic diagram was produced using R package analogue (Simpson and Oksanen 2020), and only species that occurred at relative abundances of > 5% in at least one sample were plotted. Stratigraphic zones were identified by applying a stratigraphically constrained incremental sum-of-squares (CONISS) analysis (Grimm 1987), and the significance of biostratigraphic zones were determined by a broken stick model (Bennett 1996). When calculating the planktic to benthic diatom ratio, we divided the planktic diatom counts by the combined counts of benthic-planktic and benthic taxa (together as benthic group), while excluding those tychoplanktic taxa. Hill’s N2 diversity was also plotted for the examined intervals in our stratigraphic diagram (Hill 1973).
Chrysophyte cysts and scales were counted whenever they were encountered in the diatom samples. The ratio of cysts and scales to diatom valves can be useful as a measure of lake production, as chrysophytes are typically found in higher relative abundances in oligotrophic lakes (Smol 1985). The relative abundances of scaled chrysophytes can also be used to indicate the water mixing regime since most of them are planktic species. Whole-lake primary production was inferred by spectrally-inferred sedimentary chl a, which tracks trends in whole-lake primary production in lakes (Michelutti et al. 2010).