Saskatchewan Condition Assessment of Lotic Ecosystems (SCALE): A multivariate tool for assessing the integrity of Northern Great Plains wadeable rivers and streams

We selected sites in a stratified random design to cover the range of abiotic conditions in southern Saskatchewan (Fig. 1). In this design, ecoregions were used as strata to ensure inclusion of the range of landscape (e.g., surficial geology) and climatic conditions resulting in 50 potential sites in each ecoregion. We conducted this study in the Prairie and Boreal Plain ecozones of south and central Saskatchewan, within the Cypress Upland, Mixed Grassland, Moist Mixed Grassland, Mid-Boreal Lowlands, Mid-Boreal Uplands, Aspen Parkland, and Boreal Transition ecoregions (see biodiversity.sk.ca/eco.htm).

The second strata we selected was stream order (Strahler 1964), with sites ranging from first-order streams to seventh-order streams, and we spread our sampling effort across 10 sites of each stream order in each ecoregion. Analysis of data for eighth order and higher streams of southern Saskatchewan is in Phillips et al. (2015) and are not included here. The cumulative length of stream type categorized by order in the Prairie and Boreal Plain ecozones of southern Saskatchewan is summarized below (Saskatchewan Ministry of Environment 2006):

As such, the current study estimates stream health of an estimated (83,086 + 75,297 + 15,140)/(83,086 + 75,297 + 15,140 + 1,214) = 99.3% of the stream length in southern Saskatchewan.

The two steps used in site selection were consultation and mapping with GIS. A series of consultations with provincial experts identified anthropogenic activities that are known or suspected to potentially impact stream ecosystems across the ecosystems in the region. Potential reference streams identified for sampling were then outlined on 1:250,000 maps and stream orders identified using the Strahler (1964) method based on 1:50,000 National Topographic Survey of Canada series maps (Natural Resources Canada 2018). Road crossings possessing open and free-flowing culverts or bridges over streams were identified as access points for sampling. Two or three potential sites more than 100 m upstream of each road crossing were in each potential reference stream within a subcatchment. A potential site was considered to be a stream reach with a longitudinal distance approximately six times its width (∼6−100 m; Newbury 1984). Only run reaches were sampled as they were the dominant habitat type and riffles were highly uncommon (Phillips unpubl. data). Where riffles did occur, they were typically cattle-crossings that received focused impact from livestock and would not likely reflect the broader reach-scale instream conditions. The next step was field verification involving site visit and site assessment if no discernable point source human activity was identified at the site. A total of 486 sites were visited; however, if sites were dry at time of visit, they were not sampled or included further in the study.

Reference and test sites were visited only once, in autumn (15 August–15 October, 2006–2009) as it provided fewer populations of rapidly growing and reproducing multivoltine benthic macroinvertebrate taxa, and a greater proportion of univoltine taxa representing longer-term conditions in the waterbody (Beatty et al. 2006). Further, autumn provided a period of low-water conditions across Saskatchewan allowing for greater access.

The RCA requires approximately 250 sites (Reynoldson and Wright 2000) to characterize the variability among streams and build predictive models. We assembled a total of 280 potential reference sites and sampled them over a four-year study period from 2006 to 2009 (Fig. 1), together with 8 test sites known to have a measurable human activity that could affect water chemistry, hydrology, or habitat (Table 1). Specifically, we selected a waste water stressor site that is 100 m downstream of a major urban center’s waste water treatment plant (City of Regina) on the Wascana Creek, four sites ∼100 m downstream of reservoirs on the Whitesand River, Frenchman River, Avonlea Creek, and Souris River, a site on the Wascana Creek ∼100 m downstream of a major urban area it passes through (City of Regina; but upstream of the waste water treatment plant), a site with 100% cropped agriculture in the watershed upstream on a Pipestone River tributary, and a site with high levels of polycyclic aromatic hydrocarbon (PAH) contamination on the Moose Jaw River. Although each of these sites will have co-occurring stressors in addition to the stressors identified, they represent extremes of specific human activities we could identify in this region of the NGPs.

Traveling kick-and-sweep samples (500-μm mesh) were conducted at each site along a transect perpendicular to stream banks following the methods of Phillips et al. (2016). These consisted of composite sweeps (10 seconds per sweep over an area ∼1 m) against the left bank, ¼, ½, and ¾ of the stream width to the right bank, and against the right bank relative to a perspective looking upstream. All five position sweeps were integrated into a single sample per transect. Samples were concentrated using a 500-μm mesh sieve and immediately preserved with 95% ethanol. Organisms were sorted using a stereoscope at 7 × magnification. For samples estimated to have a total number of organisms exceeding 1000 individuals, the sample was subsampled by evenly spreading the samples on a 250-μm mesh sieve and then removing half the sample to sort. Subsample counts were then multiplied by the fraction removed to estimate abundance in the original sample. Use of a Marchant-box-style subsampler was not possible due to high amounts of macrophytes and filamentous algae in many samples. Specimens were identified to the lowest possible taxon designation (usually genus and species, but family for Chironomidae, Ceratopogonidae, Stratiomyidae, Tabanidae, and Tipulidae, subcohort for Hydrachnidia, class for Oligochaeta, and phylum for Nematoda) using keys for North America (Clifford 1991; Webb 2002; Webb et al. 2004; Merritt et al. 2008). Voucher series were deposited in both the Water Security Agency Invertebrate Voucher Collection (Saskatoon, Saskatchewan) and the Royal Saskatchewan Museum (Regina, Saskatchewan).

The environmental data collected included variables characterizing time of sampling, landscape geology, hydrology, water chemistry, habitat characteristics, and land cover (Table 2). These variables were divided into two categories distinguishing those least-likely affected by human activities and those that reflect human development in the watershed. From this, the environmental attributes least affected by human activities were used to distinguish among biological community types (Table 3). Those most affected by anthropogenic activity were used to evaluate the biological community descriptors that were most indicative of specific stressor gradients (i.e., land cover variables; Table 3).

Study site locations did not align with previously delineated watershed boundaries for Saskatchewan; therefore, contributing areas upstream of study sites were delineated individually. The flow network upstream of each site was created using the ArcGIS Network Analyst extension and the Saskatchewan Stream Network (Water.stream_network_v2, Sask. Env. Stream Network version 2, circa May 2006, based on the National Topographic Series 1:50,000 water data). Canadian Digital Elevation Data (CDED, Natural Resources Canada 2016, resolution of 0.75 arc seconds) was then overlain on the flow network for each site and the ArcGIS ArcHydro Extension applied to delineate the estimated contributing watershed polygon for each of the 299 point locations. Visual inspection determined that the watershed polygons for 30% of the sites were obviously incorrect, likely due to the relative lack of accuracy of the CDED and stream network datasets. Watersheds with obvious errors were manually corrected in ArcGIS to capture the accurate watershed polygons.

Portions of many watersheds in southern Saskatchewan have both dynamic contributing areas that increase in size with the magnitude of individual runoff events (PFRA Hydrology Division 1983), and areas that do not contribute flows downstream. To account for dynamic and non-contributing areas, the watershed polygons delineated for each study site were further reduced by overlaying the effective drainage area (Agriculture and Agri-Food Canada 2013) which may be expected to entirely contribute runoff to the main river channel during a runoff event with a return period of two years (Godwin and Martin 1975). A complete discussion of drainage boundary delineation methods can be found in Hydrology Report #104 (PFRA Hydrology Division 1983). The study site watershed polygon line was snapped to the effective drainage area boundary wherever it delineated an area larger than the effective drainage area.

The evaluation of land use, land cover, and physical characteristics of watersheds upstream of each site were evaluated using the “Tabulate Area” tool in ArcMap version 10.1. Calculations of land cover in each polygon were derived from the AFSC_56m_2006 dataset produced by the Agricultural Financial Services Corporation (AFSC, Agricultural Financial Services Corporation 2006) (Table 3). Land cover variables were summarized as the % composition of the total land cover in the effective drainage area of a polygon.

Soil composition within each site watershed was quantified using Agriculture and Agri-Food Canada’s (2005) digital Saskatchewan Soil Resource Data, Saskatchewan Soil Information System Ver. 1 (SKSISv1). Ecoregion membership and composition of each site watershed were determined using the Ecoregions ArcGIS file (National Ecological Framework for Canada 2006). Surficial geology of each polygon was quantified using 1:250,000 scale GIS map (250K_surficial) from Agriculture and Agri-Food Canada Quaternary geology map series, showing surficial terrain deposits classified by depositional environment and geomorphology (Agriculture and Agri-Food Canada 2008). Soils, ecoregion, and surficial geology variables were summarized as % composition of the total effective drainage area in each polygon (Table 3).

Hydrology variables for each site were determined using Water Survey of Canada Hydrometric Stations near each site (wateroffice.ec.gc.ca/). A 1:2 year median value was calculated for each of the following variables: Minimum Mean Daily Discharge (m³ · s⁻¹), Peak Mean Daily Discharge (m³ · s⁻¹), and Annual Mean Volume (dam³). Finally, Effective Drainage Area (EDA, km²) was estimated using the two-year return period contributing area described above (PFRA Hydrology Division 1983; Table 3).

Turbidity (NTUs), Dissolved Oxygen (mg · L⁻¹), pH, and Specific Conductivity (μS · cm⁻¹) were all measured at the time of benthic macroinvertebrate collection for each site. However, Dissolved Oxygen was removed from further analysis as it is susceptible to variation in time of day or season, while Specific Conductivity was removed due to concern it could be influenced by human activity. Further, pH was removed for further analyses as it was found to have too small a range to meaningfully explain variation in benthic macroinvertebrate communities (mean ± 1 SD, 8.33 ± 0.58; Table 3).

Site habitat condition is a composite of Environmental Protection Agency (EPA) Rapid Bioassessment variables for streams and wadeable rivers (Barbour et al. 1999). Specifically, we included Embeddedness, Channel Flow Status, Sediment Deposition, Flow Alteration, and Bank Stability. Each of these variables is scored from 0 to 20 from low condition to optimal condition (see Barbour et al. 1999). All variables were then summed to a maximum score of 100 at each site that was collectively referred to here as Site Habitat Condition.

We designed and calculated a metric we termed the Cross-Sectional Habitat Availability (CSHA) at each site as a variation on hydraulic radius, in an effort to estimate how much space was available to benthic macroinvertebrates at each site (Fig. 2). To estimate this distance, we first measured the Channel Width (CW; m) based on the CABIN protocol for bankfull width (Environment Canada 2012) along the study site transect, then the depth of the channel at 1/3 Channel Width (z₁; m), ½ Channel Width (z₂; m), and 2/3 Channel Width (z₃; m) at time of site benthic macroinvertebrate collection. We then applied these dimensions in Eqs. 1–4 below to estimate benthic-distances (hypotenuses) for four sections of benthic habitat which were subsequently summed (Eq. 5) to provide total distance we termed the CSHA (see Fig. 2) as follows:where d₁ = benthic distance from left bank to 1/3 CW (m), w₁ = 1/3 CW (m), and z₁ = depth of channel at w₁ (m).where d₂ = benthic distance from 1/3 CW to ½ CW (m), w₂ = ½ CW - 1/3 CW (m), z₁ = depth of channel at 1/3 CW (m), and z₂ = depth of channel at ½ CW.where d₃ = benthic distance from ½ CW to 2/3 CW (m), w₃ = ½ CW - 1/3 CW (m), z₂ = depth of channel at ½ CW (m), and z₃ = depth of channel at 2/3 CW (m).where d₄ = benthic distance from 2/3 CW to right bank of channel (m), w₄ = 1/3 CW (m), and z₃ = depth of channel at 2/3 CW (m).where d_n = benthic distance (m) at each cross sectional division (Fig. 2).

We used the land use layer described above to quantify the amount of human activity upstream of a site in the contributing areas of each watershed and construct a collection of reference condition sites with the least amount of human activity possible. Specifically, we considered the % urban land use, % cropland, % pasture land use, and % total land under human activity (combination of urban, cropland, and pasture land), as well as the number of landfills, oil wells, and bridges or road crossings (Government of Saskatchewan 2019) as the measure of human activity in a watershed. We calculated land use % composition from the AFSC_56m_2006 dataset (Agriculture and Financial Services Corporation 2006), the number of landfills and lagoons from Saskatchewan Ministry of Environment landfill and lagoon site data (SaskH₂O 2019), and bridges or road crossings from the National Road Network dataset (National Resources Canada 2019). We established an objective of ∼90 reference sites from our 280 potential reference sites and set our best-possible reference condition based on the sites with the least amount of these seven human activity variables above. Although Reynoldson and Wright (2000) recommended that as few as five sites may be necessary per biological group, Bowman and Somers (2005) demonstrated that five reference sites provide inaccurate and imprecise measures of community and metric variation. Further, this study of the NGP covers a broad range of abiotic characteristics, and despite the benefits of having as little human disturbance as possible in selecting a small number of reference sites, doing so would reduce representation of abiotic variables and applicability of the overall model to all wadeable rivers and streams of the NGP. As such, we chose to increase the accuracy and precision of this RCA following the guidelines of Bowman and Somers (2005) to set a minimum cutoff of 20 sites for a biological group. Ultimately, selection was based on relaxing the amount of urban land use, cropland land use, pasture land use, and total land under human influence in the upstream watershed, as well as total number of landfills, oil wells, bridges, or road crossings in the contributing watershed until a minimum of 20 sites per biological group were retained for reference condition of all the available sites in the study.

The next step in an RCA study design is to identify reference site biological community groups at minimally disturbed sites. To begin with, taxa occurring in fewer than five sites or those that had a total abundance of less than 10 individuals were removed from analyses of community composition as these community summaries (e.g., ordination and cluster analysis) are sensitive to rare taxa. Community abundance data were log_(n+1)-transformed prior to analyses, and three sites were removed from further analyses as the only invertebrates present were a few (<5) Chironomidae individuals. Including these sites could have led to poor characterization of community type likely reflecting a period in which the stream dried up prior to sampling.

To achieve these biological groups, we clustered reference sites using the quantitative symmetric dissimilarity metric known as the Kulczynski distance linking metric (McCune and Grace 2002) to relativize data points by sample unit totals, since there was a large range in total abundance among samples (range = 18,464). This distance linking metric is a variation on the Bray–Curtis coefficient, but with a “built-in standardization” (Faith et al. 1987). Clustering of the reference sites into biological groups was completed using an agglomerative hierarchical fusion method with complete linkage (Clarke and Warwick 2001). Separate clusters of reference sites were then characterized as biological groups for further analysis and discussion, and an analysis of similarities (ANOSIM) test (Clarke and Warwick 2001) was used to quantify the assemblage difference or overlap between grouping outputs (performed in PRIMER 6, Version 6.1.13, PRIMER-E Ltd.).

Secondarily, it is necessary to identify abiotic characteristics that can be used to distinguish reference biological communities and then subsequently use those abiotic characteristics to assign test sites to the most appropriate reference group. Of the 43 environmental abiotic variables collected in this study, we chose 14 of these variables most likely not to be influenced by human activity to search for underlying relationships separating the reference site biological groups (date of sampling, surficial geology, soils composition, ecoregion, effective watershed area, stream order, median annual volume, median peak flow, median minimum flow, CSHA, site velocity, specific conductivity, turbidity, and pH). Basic Spearman’s rank correlations were calculated to determine redundancy of abiotic variables prior to further community abiotic variable analysis to avoid false inflation of predictive capacity and multicollinearity of predictor variables using the hmisc package in R (Harrell 2016) and an α <0.05. We used a cutoff of r > 0.6 and r < −0.6 to establish whether variables were correlated and then retained the variable with most biological relevance, established through best professional judgment focused on which variable of a highly correlated couplet is most easily accessible for applied use in the future (e.g., annual volume vs. stream order).

Discriminant functions analysis (DFA) was used to evaluate which of the remaining environmental abiotic variables (Table 3) best explain the separation between cluster analysis-based biological groups in the MASS package in R (Ripley et al. 2013). Due to unequal covariances in environmental variables, quadratic discriminant functions were used (Clarke et al. 1996).

We compiled commonly used benthic macroinvertebrate summary metrics which incorporate ecological knowledge about invertebrate assemblage responses to human stress (Resh and Jackson 1993). We conducted metric selection specific to each biological group by considering 33 metrics using lowest-taxonomic designation-level data (typically genus). We began with number of Diptera, number of Diptera families, number of EPT, % EPT, number of EPT genera, % collector gatherers, % detritivores, % filterers, % herbivores, % omnivores, % predators, % scrapers, % shredders, number of taxa, total abundance, Simpson’s diversity, Shannon’s diversity, Jaccard’s evenness, number of Coleoptera, % Coleoptera, number of amphipods, % Gammarus lacustris Sars, 1863, % Hyalella azteca (Saussure 1858), Family Biotic Index (from Barbour et al. 1999), and NMDS community Axis 1.

Next, we selected metrics for further inclusion in the model by conducting a preliminary ordination and ranking of the metrics from their correlation with the overall ordination structure and then correlating metrics with all other metrics to reduce the number of redundant metrics (Reynoldson et al. 2001). We retained a metric if the probability of its correlation with the other metrics was low (p > 0.05) for each biological group independently and retained the NMDS Axis 1 as a consistent metric in each biological group to capture changes in community that may go undetected in summary metrics.

The metrics we used in the construction of a reference database for each biological group are summarized in Table 4. For Biological Group 1, the metrics retained were number of species, total abundance, and % detritivores (Table 4). The analysis for Biological Group 2 resulted in the selection of the metrics % Hyalella azteca, total abundance, and number of Coleoptera (Table 4). Finally, % EPT, % shredders, and % detritivores were selected for representation in Biological Group 3 (Table 4).

All correlations for metric selection were performed using Spearman’s rank correlations in the hmisc package in R (Harrell 2016) and an α < 0.05.

The next step in the RCA is to compare biological metrics at test sites to those at reference sites with similar environmental characteristics deemed important in structuring biological communities in the previous modeling step.

We predicted membership of test sites to each of the biological groups, using backward discriminant functions analysis (bDFA) and compared test site metric variables to those of their appropriate reference groups using TSA (Bowman and Somers 2006). TSA formally evaluates the magnitude of difference between test sites and reference sites using a noncentral hypothesis test, defining the normal range as the probability region enclosing 95% of the reference sites (Kilgour et al. 1998). In this evaluation, a small probability (p ≤ 0.05) suggests the test site is impaired, while a large probability (p > 0.05) suggests the site is not impaired (Kilgour et al. 1998). However, further to this, we followed Bowman and Somers (2006) three-tiered condition applying the following conditions: impaired (p ≤ 0.05), possibly impaired (0.05 < p < 0.95), and within reference condition (p > 0.95).

Eight test sites known to have human activities (waste water effluent [n = 1], reservoirs [n = 4], high urban development [n = 1], high crop agriculture [n = 1], and PAH contamination [n = 1]) were selected to evaluate the performance of the four metrics in discerning impact (Table 4). There was no a priori way to determine if these test sites were actually impaired, but they were selected on the likelihood of impairment given known human development. However, comparing these sites with observed stressors against reference site groups will provide insight into whether the metrics selected are sensitive to disturbance.

We used bDFA to assign test sites to their respective biological groups of reference sites based on stream order and ecoregion parameters described above. Test sites were then evaluated against respective biological groups using TSA.

Specifically, the Wascana Creek downstream of the Regina Waste Water Treatment Plant (WWTP), was compared to the reference community of Biological Group 2 (Table 4). As a contribution to the understanding of potential impacts of reservoir operation, Avonlea Creek at the Avonlea water control structure and the Souris River downstream of the Rafferty Dam were compared to Biological Group 1, the Frenchman River at the Eastend Reservoir was compared to Biological Group 2, while the Whitesand River downstream of Theodore Dam was compared to Biological Group 3 (Table 4). A single urban-dominated site at Regina in the Wascana Creek downstream of the city but upstream of the waste water treatment was evaluated against reference Biological Group 2, and a crop-dominated site on a Pipestone River tributary was evaluated against reference Biological Group 1 (Table 4). Finally, we selected a single site contaminated with PAHs in the Moose Jaw River near an asphalt refinery as a unique site, wherein sediment concentrations of PAHs exceed the Canadian Council of Ministers for the Environment (CCME) interim sediment quality guidelines (ISQG, CCME 1999, site 11, Table 4 and 5). This provided a situation of known anthropogenic activity and a mechanism at an extreme end of perturbation to evaluate against reference condition.

We developed a second line of evidence to verify metric responses to specific stressors, by interpreting regressions between the metrics we identified above and gradients of abiotic variables within each biological group. Specifically, all sites collected in each biological group which had complete data for macroinvertebrate communities and abiotic variables were included in this analysis (Groups 1, 2, and 3 consisted of 71, 68, and 33 sites, respectively). We used both metrics that we had validated to be relevant in specific biological groups, as well as selected alternative metrics that were highly correlated with the metrics selected for the TSA. We included all metrics (% detritivores, number of species, total abundance, NMDS Axis 1, % Hyalella azteca, Log_(n+1) number Coleoptera, % shredders, and % EPT) calculated for each biological group.

Next, we included human-influenced abiotic variables for this gradient analysis as follows: number of road crossings, number of lagoons, number of landfills, and eight land cover metrics (% cropland, % urban, % pastureland, % native prairie, % hayland, % forest, % water, and % wetlands) in the upstream contributing watershed, then reduced the total number of variables using a preliminary correlation to combine correlated stressor metrics. The final list of seven abiotic variables included a single landfill/lagoon/road crossings variable, % cropland, % urban, % uncultivated (summary of native pasture, hayland, and pasture), % forest, % water, and % wetland variables to correlate against the biological metrics.

We used linear regression analysis from several R packages (R Development Core Team 2013) to illustrate the relationships between metrics and human stressor variables. Each continuous variable was log_(n+1)-transformed, while each proportional metric (e.g., % Hyalella azteca) and environmental variable (e.g., % land use) were arc sin (x/100)^0.5-transformed (Podani 2000) prior to analysis.

To further investigate the benthic macroinvertebrate taxa that characterize the relationship between NMDS Axis 1 and land use, we conducted supplementary Spearman’s rank correlations in the hmisc package in R (Harrell 2016) and an α <0.05.

Using this minimum requirement of 20 sites per biological group described above, we generated a final collection of 83 reference sites with the best possible reference condition characterized as follows: <5% urban land use, <40% cropland land cover, <50% pasture land use, < 70% total land under human influence in the upstream contributing watershed, as well as no greater than two landfills, oil wells, bridges, or road crossings. Ultimately, we were able to assign the reference sites to three biological groups (Fig. 3).

Classification based on cluster analysis of the community data at reference sites resulted in three possible biological groups of 76 taxa found in the reference sites (Fig. 4). Groups 1, 2, and 3 consisted of 25, 33, and 25 sites, respectively. A subsequent one-way ANOSIM between these three groupings provided a global R = 0.381 (p < 0.1%), indicating difference between assemblages but with some overlap (Clarke and Warwick 2001). Group 2 had the highest abundance and lowest Shannon’s diversity, and Group 3 had the lowest abundance and the highest Shannon’s diversity (Table 6). The most commonly occurring taxa in all groups were Chironomidae, but the mayfly genus Caenis spp. occurred in high proportions in each group as well (Table 6).

Biological Group 1 was characterized primarily by Aspen Parkland, Cypress Hills Upland, and Mixed Grassland ecoregions, and all sites occurred in low first- and second-order streams (Table 7). In contrast, Biological Group 2 primarily occurred in Boreal Transition, Mixed Grassland, and Aspen Parkland ecoregions, with most of the sites in third- and fourth-order streams (Table 7). Finally, Biological Group 3 was nearly entirely characterized by Boreal Transition and Aspen Parkland ecoregions, as well as larger fourth- and fifth-order rivers (Table 7).

In applying TSA to the 8 sites with known human activity, we found that two sites with upstream reservoirs in Biological Group 1 demonstrated moderately lower (possibly impaired) total abundance and number of species relative to reference (Table 8); specifically, the site on the Avonlea Creek downstream of the Roleau Dam had a total abundance p = 0.17, and number of species p = 0.81, while the site downstream of the Rafferty Dam on the Souris River had a total abundance p = 0.03 and number of species p = 0.99. Further, the community of benthic macroinvertebrates (based on NMDS Axis 1) was possibly impaired relative to reference downstream of the Rafferty Dam on the Souris River (p = 0.13), but not downstream of the smaller reservoir on Avonlea Creek (p = 1.00, Table 8). The third test site compared using TSA in Biological Group 1 had the distinct characteristic of its contributing watershed being completely dominated by cultivated cropland (Table 8), and although the total abundance was possibly impaired (p = 0.59), the other metrics are within reference (Table 8). In this unique situation, all other stressors we measured are absent and we measured no change in the community (p = 1.00, Table 8).

Furthermore, there was no relationship between % cropland cover and any of the metrics applied in our stressor gradient analysis for Biological Group 1 (Table 9). The metric % detritivores had a weak negative relationship with the number of landfill/lagoon/road crossings in the contributing watershed (R² = 0.18, p < 0.01; Table 9). Although the metric % Hyalella azteca was not included in the four metrics selected for this biological group (as it was correlated with total abundance in Biological Group 1 reference sites), it was also positively correlated with landfill/lagoon/road crossings (R² = 0.24, p < 0.001; Table 9). Percent EPT was also weakly correlated with forest cover in Group 1 (R² = 0.10, p < 0.01; Table 9).

In Biological Group 2, TSA provided that the Wascana River downstream of the Regina Waste Water Treatment Plant (WWTP) had a significantly lower Log_(n+1) number of Coleoptera and significantly different NMDS Axis 1 relative to reference sites (p = 0.01 and p < 0.001, respectively; Table 10). The Moose Jaw River site possessing high levels of PAH contamination demonstrated significantly lower Log_(n+1) number of Coleoptera (p = 0.01) and potentially stressed NMDS Axis 1 metric (p = 0.51; Table 10). The Frenchman River near Eastend, downstream of a reservoir, and the Wascana River in Regina downstream of high urban cover but upstream of the WWTP both had ecosystem health metrics not significantly different from reference sites (Table 10).

The gradient analysis for Biological Group 2 indicated that % Hyalella azteca had a weak relationship with the number of landfill/lagoons/road crossings (R² = 0.18, p < 0.01) and % cropland upstream of sites (R² = 0.11, p < 0.01; Table 9). Total abundance had a stronger positive relationship with the amount of uncultivated land upstream of sites (Fig. 5a, Table 9). NMDS Axis 1 however, demonstrated negative relationship with % uncultivated (R² = 0.26, p < 0.001; Fig. 5b) and % wetlands (R² = 0.08, p < 0.05), and positive relationships with % forest (R² = 0.12, p < 0.05; Table 9). The alternative metrics for Biological Group 2 showed that % shredders had a negative relationship with the % uncultivated land in the upstream watershed (R² = 0.15, p < 0.001), while % detritivores had a negative relationship with % forest cover in the upstream watershed (R² = 0.06, p < 0.05; Table 9). Finally, % EPT demonstrated a weak negative relationship with the number of landfill/lagoons/road crossings (R² = 0.10, p < 0.01), % uncultivated (R² = 0.15, p < 0.001) and % wetlands (R² = 0.14, p < 0.01; Fig. 5c), and a positive relationship with the % forest in the upstream watershed (R² = 0.17, p < 0.001; Table 9).

In summary, linear regression showed that for Group 2 the % uncultivated land in the contributing watershed had a significant relationship with total abundance (R² = 0.23, p < 0.001; Fig. 5a) and community structure (NMDS Axis 1; Fig. 5b; R² = 0.26, p < 0.001), while there was a significant negative relationship between % EPT and the number of wetlands upstream of a site (R² = 0.14, p < 0.01; Fig. 5c).

The benthic macroinvertebrate metrics of the Whitesand River site downstream of Theodore Reservoir in Biological Group 3 were not significantly different from reference sites (Table 11). However, the gradient analysis in Biological Group 3 showed that % EPT is weakly positively related to the % forest upstream of sites (R² = 0.14, p < 0.05; Table 9). The % shredders and % detritivores had no significant relationship with any of the abiotic variables analyzed here (Table 9). Total abundance was found to be positively related to the number of landfill/lagoons/road crossings (R² = 0.12, p < 0.05), % uncultivated (R² = 0.25, p < 0.01; Fig. 5d), and % forest (R² = 0.16, p < 0.01; Table 9); however, Log_(n+1) number of Coleoptera demonstrated the opposite response to % EPT, having a significant positive relationship with % uncultivated (R² = 0.25, p < 0.01) and negative relationship with % forest in the upstream watershed (R² = 0.19, p < 0.01; Table 9). NMDS Axis 1 had a negative relationship with % uncultivated (R² = 0.29, p < 0.001; Fig. 5e) and positive relationship with % forest in the upstream watershed (R² = 0.27, p = 0.002; Fig. 5f; Table 9). In the alternative metrics, % Hyalella azteca had a positive relationship with % uncultivated cover (R² = 0.20, p < 0.01), but a negative relationship with % forest in the upstream watershed (R² = 0.18, p < 0.0; Table 9).

There were strong significant correlations with individual taxa along the uncultivated and forest gradients in Group 3. We found the crawling water beetle Haliplus apicalis Thomson, 1868 (Coleoptera: Haliplidae) to have a strong negative correlation with the % forest cover (r = −0.674, p < 0.001), as well as Hyalella azteca (r = −0.427, p < 0.05), and the elmid riffle beetles (Coleoptera: Elmidae, r = −0.584, p < 0.01). Further, we found the uncultivated gradient to be characterized by positive correlations with the small square-gilled mayflies Caenis latipennis Banks, 1907 (Ephemeroptera, Caenidae, r = 0.711, p < 0.05), Haliplus apicalis (r = 0.624, p < 0.001), and Hyalella azteca (r = 0.516, <0.05).

Setting reference criteria for the NGPs, as for any biomonitoring program, is a compromise between historical conditions and what is realistically achievable for current circumstances. Truly pristine waterbodies untouched by human activities no longer exist in the strictest sense (Boivin et al. 2016), and therefore setting management goals for pristine conditions is untenable. As such, suggestions that current water quality and ecological health could be compared to that which would occur in the absence of human activity is not helpful and requires a more realistic view of what a biomonitoring program can achieve (Ode et al. 2016). Although this change in perspective is often characterized as a “shifting baseline,” and a lamentable disservice to the environment (Papworth et al. 2009), it is realistically all that may be possible today (Stoddard et al. 2006). There is even scientific argument for this human alteration of the environment being so great that the earth has passed into a new epoch because of it (Waters et al. 2016), and the environmental consequences captured in this “Anthropocene” are broad-reaching to global declines in biodiversity and climate change.

The limitation of minimally disturbed sites in biomonitoring is a common issue in regions of high human activity and result in a situation where exceptional ecosystem health is challenging to quantify (e.g., Bouchard et al. 2016). To accommodate this reality in management, indices of biotic integrity in Maine, Ohio, Vermont, and Minnesota have been investigating the concept of Tiered Aquatic Life Use (TALU) frameworks wherein altered river systems can be held to an attainable ecosystem health goal that recognizes the effects of river alteration but establishes goals that are realistic with further improvements to water quality management practices (Bouchard et al. 2016) along a biological condition gradient (Davies and Jackson 2006). This may be a valuable narrative to standardize condition found in Saskatchewan to that of other regions in the future and bridge interpretation of ecosystem health in the CABIN program with that of the Saskatchewan Condition Assessment of Lotic Ecosystems (SCALE) program presented here despite their respective differences in underlying community composition, habitats sampled, and ecoregions.

However, establishing site-specific objectives for community composition and diversity based on reference condition has advantages even if this condition is not pristine. Specifically, mitigation of human impact can improve metrics of ecosystem health such that diversity is greater than what is observed in comparable waterbodies. Improved condition in one site enables management of water quality, habitat, or abiotic conditions that could increase ecosystem health of other similar sites. This is particularly realistic to expect when working with reference sites that already have quite tolerant assemblages of benthic macroinvertebrates.

Though this model works in its present form, many questions and hypotheses remain. For instance, Chironomidae dominated the benthic macroinvertebrate community and are a very diverse group with ∼ 190 species known from the region from the South Saskatchewan River to Montana (Mason and Parker 1994; Parker and Glozier 2005; Phillips et al. 2013). A valuable future direction for the expansion and refinement of the model developed here would be to identify chironomids to genera or species to help refine community classification and detect impacts. As a family, Chironomidae are represented in all functional feeding groups (Merritt et al. 2008) and cover the range of tolerances provided in Barbour et al. (1999).

Consistent with our initial observations that benthic macroinvertebrate taxa in the NGPs are more tolerant than other regions where ecosystem health models have been developed, reference sites in the NGPs were typified by taxa such as Chironomidae, Oligochaetes, and Caenis sp. mayflies—all of which are more tolerant to anthropogenic stress (Barbour et al. 1999). Despite significant degradation in even the least impacted sites, we were able to develop a biological classification of reference sites that can be used to evaluate human stress at test sites and set site-specific ecosystem health objectives. Though the communities of benthic macroinvertebrates in the wadeable streams of the NGPs are not without human influence, we have constructed models for site-specific objectives that can be used to begin management and allow hypothesis testing for mitigation of anthropogenic perturbation.