The Prairie region of Canada is a key agricultural region of the nation, but one with important stressors impacting water quality and security. Within the prairies, rates of internal phosphorus (P) loading and trophic status of lakes are known to be high, and cyanobacterial blooms are common, owing in part to nutrient-rich soils, human activities in the watersheds, and often shallower lakes. However, there is important variation in the region in trophic status, internal loading, and sediment geochemistry. We characterized the variation in lake sediment composition, water chemistry, lake and watershed attributes expected to affect sediment P cycling within the prairie region (the southern portion of the provinces of Manitoba, Saskatchewan, and Alberta; Boreal Plains and Prairies Ecozones). We aim to help understand what drives variation in sediment P speciation, and porewater SRP. Multiple attributes expected to affect P cycling varied substantively between the two ecozones (e.g., lake trophic status, pH, salinity, and dissolved sulfate were higher in Prairie ecozone lakes than in the Boreal Plains ecozone). While the studied ecozones differ in terms of soils, land use, lake water chemistry, lake trophic status, and sediment P fractions, the ecozones themselves were not associated with significant differences in porewater SRP.

Field sampling was conducted during the summer of 2017, 2018, and 2019 (Huot et al., 2019), using a single sample approach intended as a snapshot of summer variability. In sum, we present data from 60 lakes for sediment P sequential fractionation, 58 lakes for select dissolved ions, 32 lakes for sediment porewater soluble reactive P (SRP), and 29 lakes for select sediment trace metals (Fig. 1). Porewater and trace metal analyses were restricted based on sampling time and cost (due to lab access limitations during the pandemic). These samples were selected using a stratified random sampling approach applied by LakePulse (Huot et al. 2019), with select samples missing due to shipping or lab issues.

Water column sampling

Water samples were collected at an ‘index site’ (usually in the deepest part of the lake). We collected the surface water samples within the euphotic zone (estimated with a Secchi disk) with a 2-metre-long integrated tube sampler (Huot et al., 2019; Canadian Lake Pulse Network, 2021). We used a Van Dorn bottle to sample the hypolimnion of stratified lakes (1m from the bottom). The collected lake water was stored in an acid-washed carboy inside a cooler with ice packs until further analysis at the mobile lab on the shore. The water samples for SRP were filtered through a 0.45 μm syringe filter and stored in the freezer at -20℃, until analyzed.

Sediment coring

Sediment cores (6.6 cm internal diameter) were collected using a National Lakes Assessment (NLA) gravity corer (USEPA, 2017) at the index site. Upon return to the shore, the cores were carefully extruded and sliced (to avoid disturbing the sediments) using a vertical extruder (Aquatic Research Instruments). The sliced sediment layers were frozen at -20℃ until analysis. For this study we used only the topmost slice (0-2 cm depth).

Porewater extraction

Sediment samples were weighed while still frozen and left to thaw for about 4 hours in a glovebox (Kraal et al., 2012). After thawing, we vacuum-filtered the sediments in a Thermo ScientificTM NalgeneTM reusable filter holder (0.45 μm) for 1-2 minutes to extract the porewater (Gruzalski et al., 2016). The remaining sediments after porewater extraction were weighed, transferred to labeled zip lock bags, stored in the freezer, and subsequently lyophilized (i.e., freeze-dried) for further analyses. The P mass extracted as SRP from porewaters was added back to the loosely bound P fraction and sediment bulk P (i.e., total phosphorus or TP) to avoid biasing the estimates of sequentially extracted P fractions, as described in (de Toledo & Baulch, 2023). Samples that did not have their porewater extracted were also lyophilized. This freezing step does have the potential to influence the sediment structure and chemistry, affecting P release (Bechmann et al., 2005; Özgül et al., 2012). For instance, Liao et al. (2019) demonstrated that freeze-thaw cycles increased phosphorus release from Yarlung Zangbo River sediments by 12%.

Soluble reactive phosphorus (SRP)

Porewater samples were transferred to acid-washed Falcon tubes, stored in the fridge, and analyzed within 24 hours of extraction using standard methods. We added a solution of ammonium molybdate and antimony potassium tartrate to the sample and used ascorbic acid to form a blue-coloured complex using a WESTCO SmartChem 170 spectrophotometer (method PHO-001-A, EPA 365.1). The SRP concentration in the lake water column was quantified at the GRIL-Université du Québec à Montreal (UQAM) analytical laboratory using a standard ascorbic acid protocol (Wetzel & Likens, 2000) on an Ultrospec 2100 pro spectrophotometer.

Sequential fractionation of sediment phosphorus

Phosphorus sequential fractionation is a technique that extracts groups of different P species (i.e., operationally defined P fractions) based on their similar reactivity towards a chemical extractant (Psenner & Pucsko, 1988). Freeze-dried samples were analyzed for total sediment P (sed TP) and P fractions at the St. Croix Watershed Research Station laboratory (MN-USA) using the method proposed by (Psenner & Pucsko, 1988) with a modification suggested by (Hupfer et al., 2009). This is one of many similar approaches to distinguish P fractions, each which have their own limitations, hence fractions are often considered “operationally defined”.  The method selected here, was used due to suitability for use across a wide range of sediment types (de Toledo & Baulch, 2023). The concentration of bulk sediment TP was quantified independently on replicate samples following sediment digestion (one hour at 85℃) with 30% H₂O₂ and a 0.5M HCl solution (Engstrom, 2005). The expected P fractions were loosely-bound P (exchangeable P susceptible to pH shifts; loose-P), redox-sensitive iron bound P (Fe-P), aluminum bound P (Al-P), labile organic P (lab-org-P), calcium bound P (Ca-P), and residual P (residual-P) (Table A1) (Psenner & Pucsko, 1988; Hupfer et al., 2009). For further description of the method and QA/QC, please see (de Toledo & Baulch, 2023). 

Bulk sediment elements 

We assessed bulk sediment composition by quantifying select elements in lake sediments (Al, Ca, Mg, Mn, Na, S, and Fe) at the laboratory of CRC in Ecotoxicology and Global Change at Université de Montréal (de Toledo & Baulch, 2023). An aliquot of 15-20 mg of freeze-dried sediments was weighed in acid-washed Teflon vials. Then, 600 𝜇L of ultra-trace metal concentration HNO₃ were added to the samples and kept at room temperature overnight. The sediments were digested in a benchtop sterilizer for 3 hours at 121℃ at 15 psi. Once the samples cooled off, 250 𝜇L of H₂O₂ trace metal grade were combined with 150 𝜇L of HCl ultra-trace metal grade, and the samples were left overnight. The digest was then transferred to a trace metal cleaned vial and analyzed by a triple quadrupole inductively coupled plasma mass spectrometer (ICP−QQQ Agilent 8900).

Lake water chemistry

Dissolved oxygen (DO) and pH were measured across the water column using a RBR multi-parameter water quality sonde (RBR multi-channel logger Maestro). For this study we used the mean pH measurements in the surface water (0-2m depth, a.k.a. epi-pH), and the mean DO concentration (mg ⋅ L^-1) in the bottom meter. The concentration of select dissolved ions (SO₄, Ca, and Na) was measured in epilimnetic water samples at the Biogeochemical Analytical Service Laboratory at the University of Alberta by Ion chromatography (US EPA Anion Method 300.1) using a Dionex DX-600 Ion chromatograph.

Geospatial data

Each lake watershed was delineated using a Canadian digital elevation model (DEM) and a flow direction calculation through the Arc Hydro extension in ArcGIS software (Huot et al., 2019). The land use raster was created from various open data sources (e.g., Canvec manmade features, Canvec resources management features, Agriculture and Agri-food Canada Annual Crop Inventory and land use 2010, and USDA Crop Data Layer) (Huot et al., 2019). The delineation reflects the entire watershed of each study lake. The land use/cover classes were grouped into simplified categories (e.g., different types of cultivated crops were grouped into the agriculture category) (Huot et al., 2019). We grouped the soil types obtained from the slc_v2r2_canada.shp shapefile (Soil Landscapes of Canada, version 2.2 at https://sis.agr.gc.ca/cansis/nsdb/slc/v2.2/index.html) into higher soil orders (Soil Classification Working Group, 1998). The surface geology data were obtained from the website of the Geological Survey of Canada (GEO_POLYS.shp shapefile, Canadian Geoscience Map 195,2014, 1 sheet, https://doi.org/10.4095/295462). We used QGIS to plot the maps (QGIS Association, 2022).

Data analysis Summary statistics and comparisons

We used JASP (Love et al., 2019) to produce summary statistics for the concentration of porewater SRP, the concentration of bulk sediment P (sediment TP) and the extracted P-fractions (Table A2), the % area of land use, soil orders, and surface sediments in the watershed (Table A3), the concentration of dissolved ions and bulk sediment elements (Table A4), the concentration of TP in surface waters (Table A5), lake pH (Table A5), DO in the bottom meter (Table A5), salinity (Table A5), and lake area and depth (Table A5). Furthermore, we compared the Log₁₀-transformed concentration of porewater SRP between the prairie region (Boreal Plains and Prairies ecozones) and the remainder of other Canadian ecozones, between the Boreal Plains and Prairies ecozones, and among different classes of lake depth (measured at the index site), lake area, and trophic state using the non-parametric Mann-Whitney and Kruskal-Wallis tests, as the data violated assumptions for parametric tests (e.g., normal distribution and variance equality).  We ran Mann-Whitney and Kruskal-Wallis tests to assess the variability of individual P-fractions and other variables between the two ecozones and across different classes of trophic state, followed by pairwise post-hoc comparisons (Dunn’s, corrected for multiple comparisons; Holm method). In addition to that, we applied the Spearman’s rank correlation to assess the relationship between porewater SRP and all the other variables (excluding bulk sediment P and the extracted P fractions; but see PCA for associations). These statistical tests and raincloud plots were all performed in JASP (Love et al., 2019). All analyses met assumptions of the respective tests, once data transformations were completed. 

Analysis of Similarities (ANOSIM) to assess variability of sediment P composition

We treated sediment P composition as an assemblage of the sequentially extracted P fractions representing each lake in a 6-dimensional space. Then, we used Analysis of Similarities – ANOSIM (Clarke, 1993) on PAST software (Hammer et al., 2001) to compare the distance among lakes from the same groups and among different groups (e.g., ecozones, trophic state, depth, and surface area). Despite its name, ANOSIM is a non-parametric multivariate technique that tests the dissimilarity between two or more groups (Clarke, 1993; Zuur et al., 2007). It is a nonparametric test, requiring continuous variables. We selected ANOSIM because of its high sensitivity in detecting differences between the group’s centroids, dispersion, degree of skewness, and correlation structures (Anderson & Walsh, 2013). We selected Bray-Curtis as a distance measure because it is a proportion coefficient that emphasizes the similarity in species composition (i.e., extracted P-fractions) among samples (McCune & Grace, 2002; Anderson & Walsh, 2013). The mean within-group ranked distance was compared with the mean between-group ranked distance (Clarke, 1993; Anderson & Walsh, 2013) by calculating an ANOSIM R statistic, which is constrained to the range -1 to 1. Positive values suggest that dissimilarity is higher between-groups than within-groups, whereas negative values indicate that dissimilarity is higher within-groups than between-groups (Clarke, 1993). The test significance was assessed by 9999 permutations of group membership, and a post-hoc test performed pairwise comparisons among all groups, generating Bonferroni-corrected p−values (Hammer et al., 2001). The null hypothesis tested was there is no difference between groups (Clarke, 1993; Anderson & Walsh, 2013).

In summary, we used ANOSIM to compare the six extracted P-fractions as an aggregate indicator of lake sediment P composition between the Boreal Plains (n=32 lakes) and Prairie (n=28 lakes) ecozones and across different groups of lakes based on morphometry (e.g., lake depth and surface area),  and lake trophic state (e.g., based on TP concentration measured on the water column). We used three lake depth classes, measured at the index site (shallow: < 4.35 m depth n=36, intermediate: 4.4 – 10 m depth n=17, and deep: 14 – 20.5 m depth n=7) and three lake size classes (small 0.1 – 0.5 km² n=17; medium 0.5 – 5 km² n=23; large >5 km² n=20). Additionally, we compared sediment P composition among four different classes of lake trophic state, based on the TP concentration measured on the water column (mesotrophic 10−20 μg ⋅ L^-1 n=4, meso-eutrophic 20−35 μg ⋅ L^-1 n=6, eutrophic 35−100 μg ⋅ L^-1 n=21, and hyper-eutrophic >100 μg ⋅ L^-1 n=27; (Canadian Council of Ministers of the Environment, 2004). ANOSIM can be sensitive to unequal group sizes (Anderson & Walsh 2013); hence results should be interpreted with caution.

Detrended Correspondence Analysis (DCA) and Principal Components Analysis (PCA)

We ran Detrended Correspondence Analysis – DCA (Hill & Gauch, 1980) on PAST (Hammer et al., 2001) to estimate the length of the environmental gradient before selecting an ordination technique (Hill & Gauch, 1980; McCune & Grace, 2002; ter Braak & Šmilauer, 2015). As the length of our environmental gradient was short (1 unit of standard deviation), we selected Principal Components Analysis (PCA) because it is based on a linear response model, which is the most appropriate response model for analyzing short gradients (ter Braak & Prentice, 1988; Lepš & Šmilauer, 1999; Legendre & Gallagher, 2001; Kenkel, 2006; Zuur et al., 2007). Principal Components Analysis – PCA (Pearson, 1901), is one of the most widely used ordination techniques aimed at multidimensional data reduction and extraction of the strongest data patterns (McCune & Grace, 2002). As PCA is also referred to as an indirect gradient analysis (Zuur et al., 2007), it outputs the variable loading for each axis, which supports the interpretation of the axes as gradients, and the biplots allow the identification of relationships among variables and lakes.

We applied PCA on a set of response variables comprised of the Log₁₀-transformed concentration of bulk sediment P (sediment TP) and the extracted sediment P-fractions. We selected the correlation matrix (normalized variance-covariance), which means that all variables were standardized by division by their standard deviation (Hammer et al., 2001; Bialik et al., 2021). Subsequently, we added supplementary variables, and alternated among three different sets (Table A3, Table A4, and Table A5). As these supplementary variables were not used in the ordination, they did not influence the results. But their correlation with the environmental gradients represented by the PCA axes was calculated and plotted as vectors on the ordination diagram (Hammer et al., 2001). We plotted three different sets of supplementary variables onto the ordination: one comprised the percentage land use, soil orders and surface geology (i.e., representing watershed variables), the second included bulk sediment elements and dissolved ions, and the third one contained lake water chemistry indicators (pH, dissolved oxygen – DO in the bottom meter, and salinity) and lake morphometry (lake area and depth, measured at the index site) (Table A3, Table A4, and Table A5). We also included the concentration of porewater SRP in each one of the sets to assess its relationship with the other supplementary variables, besides the extracted P fractions. Additionally, all the supplementary variables, with the exception of pH, were centred and standardized by column mean and standard deviation to improve the assessment of their relationship with each ordination axis (McCune & Grace, 2002; Kenkel, 2006). Some of the supplementary variables used in PCA had missing data. Missing data were handled by iterative imputation in PAST software. So, PCA was first run with missing values replaced by their column means. Then, the program calculated regression values for the missing data, using an iterative process until reaching convergence (Ilin & Raiko, 2010).