Open access

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

Publication: FACETS
30 March 2023

Abstract

Human access to surface water resources in the Northern Great Plains (NGP) is challenged by availability and quality, and ecosystem health objectives for these characteristics have not been well developed. Here, we present a predictive multivariate model using the reference condition approach to inform goals for ecosystem health assessment. Benthic communities and abiotic variables were collected at 280 potential reference sites and 8 test sites, and of these, reference sites with least amount of human activity (n = 83) were classified into three community groups and summary metrics. Discriminant function analysis and cross-validation determined that stream order and ecoregion predicted 68.7% of the sites correctly, thus enabling comparison of sites with unknown condition to reference site groups. We then evaluated metrics through Test Site Analysis and stressor gradient analysis in each biological group. Beetle and amphipod fauna were found to be important for condition assessment in addition to traditional metrics of species richness, abundance, detritivory, Ephemeroptera/Plecoptera/Trichoptera dominance, and assemblage composition. These results provide least disturbed reference condition and ecological insights into land use impacts in the NGP. Ultimately, this model is an effective tool for evaluating biotic condition, enables prioritizing river management strategies, and can quantify the efficacy of mitigation measures.

Introduction

Water security is threatened globally by stressors brought about by resource extraction, waste water, widespread land cover change, urbanization, and engineered works such as reservoirs (Vörösmarty et al. 2010). The development of methods to monitor the impact of these anthropogenic perturbations efficiently and effectively are crucial to conserve ecosystem health for future use (Vörösmarty et al. 2010). Communities of aquatic organisms and the traits they possess are often used to evaluate the biological condition of freshwater habitats as they are responsive to changes in environmental characteristics and express the ultimate ecological consequences of anthropogenic perturbation (Resh and Jackson 1993). Benthic macroinvertebrates in particular have been used in biomonitoring studies since the early 1900s when summary metrics of their sensitivity to organic pollution were incorporated into a European-focused Saprobien system of river condition (Kolkowitz and Marsson 1908). Since this initial application, many multimetric (e.g., Karr et al. 1986; McCormick et al. 2001; Klemm et al. 2003), Nearest-Neighbour Analysis (e.g., Linke et al. 2005), Bayesian (e.g., Qian et al. 2003), and multivariate (e.g., Wright 1984; Hawkins et al. 2000; Reynoldson et al. 2001; Bowman and Somers 2006) approaches to biomonitoring have been developed and applied globally with the intention of quantifying ecosystem health in measures that can be communicated among environmental managers and the public. However, regardless of method used, the definition of biological condition and the ability to assess impairment from healthy condition relies on an understanding of what community of organisms and their traits are expected in the absence of human activity (Ode et al. 2005). These baselines are categorized commonly as reference condition in biological monitoring, and the method of evaluating condition against reference condition has been labeled as the reference condition approach (RCA; Wright 1995).
A fundamental challenge in the RCA, or any control baseline in environmental management for that matter, is to define a reference benchmark such as historical condition, least disturbed condition, minimally disturbed condition, or best attainable condition (Stoddard et al. 2006). Historical condition is a reference prior to human activity, least disturbed accepts some human activity and seeks the best available condition, while minimally disturbed condition requires only minor levels of human activity, and best attainable condition is the most optimistic condition given mitigation intervention at least disturbed sites (Stoddard et al. 2006). Although historical, minimal and best attainable conditions can be appealing for their relevance to what condition could be in the absence of human activity, on a practical level, many regions such as the Northern Great Plains (NGPs) do not yet have adequate historical condition assessment, representation of minimally impacted sites, nor examples of best attainable conditions given mitigation solutions (Stoddard et al. 2006). Therefore, the least disturbed condition is likely the most realistic characterization of reference ecosystem health possible today.
Aquatic macroinvertebrate community (species composition, abundance of individuals, and species diversity) change depends on natural characteristics of hydrology and fluctuation of flow condition (Poff et al. 1997), as well as the type and amount of water chemistry and pollution present in the system (Hilsenhoff 1988). Measurable, predictable change enables development of ecological health metrics. Other jurisdictions (e.g., United States, Britain, European Union, and Australia) have developed macroinvertebrate aquatic health measures that together with water chemistry, are used to establish and monitor surface water quality. However, macroinvertebrate measures for the NGPs have been more difficult to develop, because the natural extremes (e.g., winter, drought, erosion, high nutrient and productivity; Covich et al. 1997) of this region result in a community characterized by taxa already more tolerant to abiotic extremes than those found elsewhere in North America (Williams 1985). Therefore, macroinvertebrates in the NGP may be insensitive to pollution relative to other regions. For example, benthic metrics such as % EPT developed to respond negatively to increasing turbidity elsewhere in North America (Roy et al. 2003; Freeman and Schorr 2004; Evans-White et al. 2009) have been found to respond positively to increasing turbidity in the NGPs due to differing composition of taxa in the metric (e.g., burrowing mayflies, Phillips et al. 2016). As such, feasibility of an aquatic health measure using NGP macroinvertebrates requires investigation before existing metrics from other jurisdictions are applied, and these metrics must be validated for particular stressors to account for variability across ecoregions.
Rivers serve as the principal source of renewable freshwater supply for humans and freshwater ecosystems globally (World Water Assessment Programme 2009) and are of great importance in the NGPs (Dodds et al. 2004). Despite this, the NGPs in western Canada possess a landscape largely devoid of unaltered historical conditions and the tools to assess impacts of potential anthropogenic perturbations of lotic aquatic ecosystems have not been developed. Intensive agriculture, water management, and development have impacted most rivers and streams since European colonization, and the stressors these activities impart do not have ecosystem health measurement tools to reflect the impact they may be having. Here, we developed a benthic macroinvertebrate-based, multivariate biomonitoring tool using least disturbed reference conditions and Test Site Analysis (TSA) condition assessment for wadeable rivers and streams in the NGPs. The steps in its development are as follows:
1.
Characterize scope of stream ecosystems and human stressors to include in the study.
2.
Define least disturbed condition based on stressor gradients.
3.
Classify biological communities at reference sites and determine what underlying abiotic characteristics best discriminate between classifications. These relationships are used to match test sites of interest to appropriate reference groups.
4.
Select summary metrics of the benthic macroinvertebrate community that capture significant community structure and eliminate redundant descriptors.
5.
Designate test sites to appropriate reference groups using predictive modeling and evaluate biological condition based on metrics identified above using TSA.
Our goal was to apply a RCA based on benthic macroinvertebrate communities to develop an ecosystem health monitoring tool for NGPs rivers that can be applied to assess the effects of anthropogenic perturbation. This study develops a predictive model to evaluate the amount of variation in the benthic macroinvertebrate communities that can be explained by abiotic characteristics of rivers and streams at reference sites with the least amount of human activity. Further, we apply these abiotic characteristics to assign test sites to appropriate reference benchmarks and evaluate biological condition based on community metrics to predict what ecosystem health metrics should resemble under the least amount of human activity.

Methods

Study region

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).
Fig. 1.
Fig. 1. Map of the study area in Saskatchewan, Canada. Black triangles indicate location of study sites within ecoregions across southern Saskatchewan. Figure 1 was built at the Water Security Agency, Saskatoon, Canada, using a NAD83 Universal Transverse Mercator Projection. The main layer (Ecoregions) was downloaded from the National Ecological Framework for Canada (2006): A National Ecological Framework for Canada – GIS data – Agriculture and Agri-Food Canada (AAFC)). The supporting layers are downloaded from the Statistics Canada 2016 site (2016 Census Boundary files (statcan.gc.ca). All basemap files are provided under the Open Government Licence – Canada | Open Government, Government of Canada and the Statistics Canada Open Licence (statcan.gc.ca).
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):
1.
Ephemeral streams: Classified as stream order 0–1, comprising 83,086 km.
2.
Temporary streams: Classified as stream order 2–4, comprising 75,297 km.
3.
Perennial rivers: Classified as stream order 5–7, comprising 15,140 km.
4.
Large rivers: Classified as stream order 8, comprising 1,214 km.
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.
Table 1.
Table 1. List of test sites with their respective stressors and biological grouping membership.
Test SiteStressorRiverLocationBiological Group
1Waste waterWascana RiverRegina2
2ReservoirWhitesand RiverTheodore Dam3
3 Frenchman RiverEastend Reservoir2
4 Avonlea CreekAvonlea1
5 Souris RiverRafferty Dam2
6UrbanWascana RiverGolf Course2
7100 % CropPipestone River Trib.Moosomin1
8PAH* contaminationMoose Jaw CreekMoose Jaw2
*
PAH – polycyclic aromatic hydrocarbons.

Biological data

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).

Environmental data

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).
Table 2.
Table 2. Land cover variables evaluated in watersheds upstream of benthic macroinvertebrate collection sites throughout southern Saskatchewan for evaluation of reference condition criteria and metric evaluation.
Land cover typeDescription
Annual CroplandLand used for annual cropping
Native PastureInvcludes native and seeded grazing land but not riparian areas
Improved PastureIncludes native and seeded grazing land but not riparian areas
HayLand used for cut forage (alfalfa, clover, grass, mix)
ForestTreed land including cutovers and forest burns
WetlandSaturated landscapes with wetland vegetation species
WaterPermanent waterbodies
BarrenNon-vegetated areas including badlands, salt/mud flats, and industrial facilities
Built upUrban and populated areas
Table 3.
Table 3. Physical abiotic variables of benthic macroinvertebrate assemblage reference sites in rivers and streams of southern Saskatchewan, and retention for further analysis in reference site classification after evaluation for inclusion in reference site biological group abiotic discrimination based on collinearity.
VariableUnitsMeasurement ScaleVariable CharactiristicsRetained?
Time 
 Date of samplingn/aDay of year0–365N
Landscape    
 Surficial geologyPCA Axis 1Contributing Watershed0–100Y
 Soils compositionPCA Axis 1Contributing Watershed0–100N
 EcoregionnoneContributing Watershed7 categoriesY
Hydrology    
 Effective watershed aream2WatershedContinuousN
 Stream ordernoneSite6 categoriesY
 Median annual volumedam3SiteContinuousN
 Median peak flowm3 · s−1SiteContinuousN
 Median minimum flowm3 · s−1SiteContinuousN
 Cross-Sectional Habitat Availabilitym2SiteContinuousN
 Site velocitym · s−1SiteContinuousN
Water chemistry    
 Specific conductivityμS · cm−2SiteContinuousN
 TurbidityNTUSiteContinuousY
 pH SiteContinuousN
 Dissolved Oxygenmg · L−1SiteContinuous*
 SalinitypptSiteContinuous*
Habitat    
 Site habitat conditionnoneSiteContinuous*
 Site depthmSiteContinuous*
 EmbeddednessScoreSite0–20*
 Flow statusScoreSite0–20*
 Sediment Deposition scoreScoreSite0–20*
 Bank stabilityScoreSite0–20*
 % Canopy cover% CoverageSite4 categories*
 Substrate classNoneSite42 categories*
 Riparian vegetative communityNoneSite18 categories*
 Macrophyte typeNoneSite4 categories*
 Algae typeNoneSite4 categories*
 Woody debris abundanceNoneSite3 categories*
 Macrophyte abundanceNoneSite3 categories*
 Detritus abundanceNoneSite3 categories*
 Algae abundanceNoneSite3 categories*
 Substrate composition% CoverageSite% of 5 categories*
 Vegetated banks% CoverageSite0–100*
 Snags% CoverageSite0–100*
Land cover    
 Annual cropland% CoverContributing Watershed0–100*
 Native pasture% CoverContributing Watershed0–100*
 Improved pasture% CoverContributing Watershed0–100*
 Hay% CoverContributing Watershed0–100*
 Barren% CoverContributing Watershed0–100*
 Built-up% CoverContributing Watershed0–100*
 Wetlands% CoverContributing Watershed0–100*
 Water% CoverContributing Watershed0–100*
 Forest% CoverContributing Watershed0–100*

Note: *Not considered due to assumed possible relationship with anthropogenic perturbation.

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 (m3 · s−1), Peak Mean Daily Discharge (m3 · s−1), and Annual Mean Volume (dam3). Finally, Effective Drainage Area (EDA, km2) was estimated using the two-year return period contributing area described above (PFRA Hydrology Division 1983; Table 3).
Turbidity (NTUs), Dissolved Oxygen (mg · L−1), pH, and Specific Conductivity (μS · cm−1) 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 (z1; m), ½ Channel Width (z2; m), and 2/3 Channel Width (z3; m) at time of site benthic macroinvertebrate collection. We then applied these dimensions in Eqs. 14 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:
d1=w12+z12
(1)
where d1 = benthic distance from left bank to 1/3 CW (m), w1 = 1/3 CW (m), and z1 = depth of channel at w1 (m).
d2=w22+(z22z12)
(2)
where d2 = benthic distance from 1/3 CW to ½ CW (m), w2 = ½ CW - 1/3 CW (m), z1 = depth of channel at 1/3 CW (m), and z2 = depth of channel at ½ CW.
d3=w32+(z22z32)
(3)
where d3 = benthic distance from ½ CW to 2/3 CW (m), w3 = ½ CW - 1/3 CW (m), z2 = depth of channel at ½ CW (m), and z3 = depth of channel at 2/3 CW (m).
d4=w42+z32
(4)
where d4 = benthic distance from 2/3 CW to right bank of channel (m), w4 = 1/3 CW (m), and z3 = depth of channel at 2/3 CW (m).
CSHA=d1+d2+d3+d4
(5)
where dn = benthic distance (m) at each cross sectional division (Fig. 2).
Fig. 2.
Fig. 2. Cross-Sectional Habitat Availability estimation from Equation 5: CSHA = d1+d2+d3+d4, where d is the distance (m), w is width (m), and z is depth for each two-dimensional polygon across the width of the river site for benthic macroinvertebrate collection.

Reference site criteria

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 (SaskH2O 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.

Classification of reference site biological data

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.).

Reference site biological group abiotic discrimination

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).

Metric selection

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).
Table 4.
Table 4. Metrics and their correlations with Axis 1 of a non-metric multidimensional scaling ordination of the sites-by-taxa dataset from reference sites in each biological group.
  Correlation of metric with ordination structureCorrelation of metric with highest-ranked metric
Biological GroupMetricrprp
1Number of Species−0.5680.0031.000<0.001
 Total Abundance−0.5530.0040.3730.066
 % Detritivores−0.4680.0410.3470.090
2% Hyalella azteca0.734<0.0011.000<0.001
 Total Abundance−0.541<0.0010.0940.610
 Log(n+1) No. Coleoptera−0.4150.0160.0120.946
3% EPT−0.613<0.0011.000<0.001
 % Shredders0.5930.0020.1930.365
 Detritivores0.4130.0450.3610.083
All correlations for metric selection were performed using Spearman’s rank correlations in the hmisc package in R (Harrell 2016) and an α < 0.05.

Prediction of metric values at reference sites

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).

Evaluation of metric performance

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.
Table 5.
Table 5. Sediment polycyclic aromatic hydrocarbons (PAHs) from test site 8 in the Moose Jaw River at Moose Jaw, Saskatchewan, Canadian Council of Ministers of the Environment (CCME) interim sediment quality guidelines (ISQGs) and probable effect levels (PELs).
  CCME
 Current studyISQGsPEL
Parent PAHμg · kg−1μg · kg−1μg · kg−1
Naphthalene36034.6391
2-methylnaphthalene104020.2201
Acenaphthylene<190*5.9128
Acenaphthene2056.788.9
Fluorene43021.2144
Phenanthrene73041.9515
Anthracene55146.9245
Fluoranthene880111.02355
Pyrene35053.0875
Benz(a)anthracene510031.7385
Chrysene899057.1862
Benzo(a)pyrene510031.9782
Dibenz(a,h)anthracene20006.2135
*
Sample less than minimum detection limit for PAH.
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.

Results

Reference Site Criteria

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).
Fig. 3.
Fig. 3. Cluster analysis of 83 reference sites in southern Saskatchewan. Data are log(n+1)-transformed, and the cluster analysis is completed using Kulczynski distance linking metric with complete linkage. Red boxes indicate the similarity cutoff used to separate the three biological groups.

Classification of reference site biological data

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).
Fig. 4.
Fig. 4. Distribution of three biological groups resulting from non-metric multidimensional scaling ordination of reference sites based on community structure. Log(n+1)-transformed, Manhattan Distance Linking Metric. Stress = 0.18.
Table 6.
Table 6. Taxa richness (mean number of taxa [±SD] per sample), total abundance (mean abundance [±SD] per sample), Shannon’s diversity (mean ± SD), Jaccard’s evenness (mean ± SD), and taxa considered as very common among the three groups of reference sites generated from cluster analysis of 76 taxa from southern Saskatchewan rivers and streams.
Biological groupSpecies richnessTotal abundanceShannon’s diversityJaccard’s evennessCommon taxa (occuring in >60% sites in a group)%*
116.6 (4.2)1429.6 (1680.7)1.6 (0.5)0.6 (0.1)Chironomidae92
 Oligochaeta80
     Hydrachnidia80
     Hyalella azteca76
     Ceratopogonidae72
     Physidae64
     Caenis spp.60
214.6 (5.3)2296.8 (3997.3)1.4 (0.5)0.5 (0.2)Chironomidae97
 Caenis spp.91
     Hyalella azteca88
     Enallagama/Coenagrion spp.73
     Gammarus lacustris67
315.9 (6.7)326.6 (419.1)1.9 (0.4)0.7 (0.1)Chironomidae96
 Caenis spp.84
     Heptageniidae64
*
% of sites in that biological group possessing each taxa.

Underlying abiotic variable correlations

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).
Table 7.
Table 7. Ecoregion and stream order site membership for Biological Groups (number of sites) from reference river and stream sites in southern Saskatchewan.
 Biological group
Physical feature123
Ecoregion
 Cypress Hills Uplands520
 Mixed Grassland8112
 Moist Mixed Grassland121
 Aspen Parkland1077
 Boreal Transition11114
 Mid Boreal Upland001
Stream order
 11801
 2741
 30144
 401211
 5038
 6010

Metric sensitivity and evaluation of test site condition

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).
Table 8.
Table 8. Comparison of impairment at test sites known to be impacted from human activity compared against reference sites in reference Biological Group 1.
 Avonlea Creek RoleauSouris River Rafferty DamPipestone Creek Tributary Moosomin
Watershed Stressors   
 % Urban land use000
 % Crop land cover6569100
 % Pasture land use860
 No. Landfills4160
 No. Road crossings2838710
 No. Lagoons290
Site-Specific StressorReservoirReservoir100 % Crop Agriculture
 Mechanism of stressHydrologyHydrologyMultiple
D in metric-based TSA2.93.22.3
p-Values for individual metrics   
 Log(n+1) total abundance0.17*0.03*0.59*
 Number of species0.81*0.99*1.00
 % Detritivores1.001.001.00
 NMDS Axis 11.000.131.00
 Total <0.05010
*
A direction of change from reference that is lower.
Simply change from reference.
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 (R2 = 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 (R2 = 0.24, p < 0.001; Table 9). Percent EPT was also weakly correlated with forest cover in Group 1 (R2 = 0.10, p < 0.01; Table 9).
Table 9.
Table 9. R2 values of metric and land use regressions in each biological group from river and stream sites in southern Saskatchewan.
    % Cover
 Validated metrics              
Biological groupAlternative metricsLandfill/lagoons/roadsCroplandBuilt upUncultivatedForestWaterWetlands
1% Detritivores0.18**        0.06*  
 Number of species              
 Total abundance              
 NMDS Axis 1              
 % Hyalella azteca0.24***        0.07*  
 % Shredders              
 % EPT        0.10**    
 Log(n+1) No. Coleoptera              
2% Hyalella azteca0.13**0.11**          
 Total abundance      0.23***      
 Log(n+1) no. Coleoptera              
 NMDS Axis 1      0.26***0.12**  0.08*
 % Shredders      0.15***      
 % Detritivores        0.06*  0.08*
 % EPT0.10**    0.10**0.17***  0.14**
 Number of species              
3% EPT        0.14*    
 % Shredders              
 % Detritivores              
 NMDS Axis 1      0.29***0.27**    
 % Hyalella azteca      0.20**0.18*    
 Total abundance0.12*    0.25**0.16**    
 Log(n+1) No. Coleoptera      0.25**0.19**    
 Number of species              

Note: *p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001.

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).
Table 10.
Table 10. Comparison of impairment at test sites known to be impacted from human activity compared against reference sites in reference Biological Group 2.
 Wascana River ReginaMoose Jaw River Moose JawFrenchman River EastendWascana River Regina
Watershed Stressors    
 % Urban land use3003
 % Crop land cover75644175
 % Pasture land use3573
 No. landfills1225912
 No. road crossings53117391324514
 No. lagoons1414413
Site-specific stressorWaste waterPAH contaminationReservoirUrbanization
 Mechanism of stressAmmoniaPAH toxicityHydrologyMultiple
D in metric-based TSA258.03.52.02.1
p-Values for individual metrics    
 Total abundance1.001.001.001.00
 Log(n+1) number of Coleoptera0.01*0.01*1.001.00
 % Hyalella azteca1.001.001.001.00
 NMDS Axis 1<0.0010.510.990.98
 Total <0.052100
*
A direction of change from reference that is lower.
Simply change from reference.
The gradient analysis for Biological Group 2 indicated that % Hyalella azteca had a weak relationship with the number of landfill/lagoons/road crossings (R2 = 0.18, p < 0.01) and % cropland upstream of sites (R2 = 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 (R2 = 0.26, p < 0.001; Fig. 5b) and % wetlands (R2 = 0.08, p < 0.05), and positive relationships with % forest (R2 = 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 (R2 = 0.15, p < 0.001), while % detritivores had a negative relationship with % forest cover in the upstream watershed (R2 = 0.06, p < 0.05; Table 9). Finally, % EPT demonstrated a weak negative relationship with the number of landfill/lagoons/road crossings (R2 = 0.10, p < 0.01), % uncultivated (R2 = 0.15, p < 0.001) and % wetlands (R2 = 0.14, p < 0.01; Fig. 5c), and a positive relationship with the % forest in the upstream watershed (R2 = 0.17, p < 0.001; Table 9).
Fig. 5.
Fig. 5. Relationship between (a) total abundance (Log[n+1]-transformed) and % uncultivated land (arc sin [n/100]0.5-transformed), (b) NMDS Axis 1 and % uncultivated land (arc sin [n/100]0.5-transformed), and (c) % EPT (arc sin [n/100]0.5-transformed) and % wetlands (arc sin [n/100]0.5-transformed) in Biological Group 2. Relationship between (d) total abundance (Log[n+1]-transformed) and % uncultivated land (arc sin [n/100]0.5-transformed), (e) NMDS Axis 1 and % uncultivated land (arc sin [n/100]0.5-transformed), and (f) % forest cover (arc sin [n/100]0.5-transformed) in Biological Group 3. The solid center line is the regression line, while the inner bounding lines are the 95% confidence band and the outer bounding lines are the 95% prediction bands.
In summary, linear regression showed that for Group 2 the % uncultivated land in the contributing watershed had a significant relationship with total abundance (R2 = 0.23, p < 0.001; Fig. 5a) and community structure (NMDS Axis 1; Fig. 5b; R2 = 0.26, p < 0.001), while there was a significant negative relationship between % EPT and the number of wetlands upstream of a site (R2 = 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 (R2 = 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 (R2 = 0.12, p < 0.05), % uncultivated (R2 = 0.25, p < 0.01; Fig. 5d), and % forest (R2 = 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 (R2 = 0.25, p < 0.01) and negative relationship with % forest in the upstream watershed (R2 = 0.19, p < 0.01; Table 9). NMDS Axis 1 had a negative relationship with % uncultivated (R2 = 0.29, p < 0.001; Fig. 5e) and positive relationship with % forest in the upstream watershed (R2 = 0.27, p = 0.002; Fig. 5f; Table 9). In the alternative metrics, % Hyalella azteca had a positive relationship with % uncultivated cover (R2 = 0.20, p < 0.01), but a negative relationship with % forest in the upstream watershed (R2 = 0.18, p < 0.0; Table 9).
Table 11.
Table 11. Comparison of impairment at test sites known to be impacted from human activity compared against reference sites in reference Biological Group 3.
 Whitesand River Theodore Reservoir
Watershed Stressors 
 % Urban land use0
 % Crop land cover50
 % Pasture land use0
 No. landfills9
 No. road crossings271
 No. lagoons4
Site-specific stressorReservoir
 Mechanism of stressHydrology alteration
D in metric-based TSA0.9
p-Values for individual metrics 
 % EPT1.00
 % Shredders1.00
 % Detritivores1.00
 NMDS Axis 11.00
 Total <0.050

Note: *A direction of change from reference that is lower. Simply change from reference.

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).

Discussion

The first objective of this study was to determine if variation in benthic macroinvertebrate communities could be explained by minimally impacted, abiotic characteristics of rivers and streams in the highly anthropogenically modified NGPs region. We found that stream order and ecoregion predicted 68.7% of minimally impacted sites correctly using cross-validation, which is comparable to previous studies where abiotic geographical predictors have been found to be useful in explaining waterbody to waterbody variation of benthic macroinvertebrate communities (Corkum 1989; Bailey et al. 1998). Applying easily measured abiotic geographical descriptors to predict components of variation in benthic macroinvertebrate communities can greatly improve biological assessments by refining the predictions of biological communities expected at a test site in the absence of the stressors of interest. Similar to our findings, Reynoldson et al. (2001) found that stream order and ecoregion variables were predictive of biological reference site group classification in Fraser River Watershed. The Fraser River catchment in British Columbia has very different ecoregions and environmental conditions relative to that of southern Saskatchewan, suggesting these variables may have broader application as predictive of communities in the absence of significant human activity.
In Canada, the BEnthic Assessment of SedimenT (BEAST) (Reynoldson et al. 1995; Resh et al. 2000) approach formed the foundation of multivariate bioassessment models that have since expanded and fall under an umbrella of the Canadian Biomonitoring Network (CABIN). This network is the national biomonitoring program recommended by Environment and Climate Change Canada and uses the RCA as its common study design (Environment Canada 2012) but is not applicable to the NGPs portion of Saskatchewan as this region typically lacks the riffle habitat necessary for the CABIN protocols and the CABIN dataset has not expanded to represent lotic waterbodies in the prairies. However, the TSA and RCA model developed here fills the gap in coverage provided by the CABIN program.
TSA is equivalent to BEAST when only ordination metrics are used but is an extension of BEAST when other common metrics (e.g., abundance) are incorporated (Bowman and Somers 2006) as was the case in this study. The use of simple abiotic variables to link test sites and reference groups (stream order and ecoregion) combined with easy to calculate metrics of the benthic macroinvertebrate community allows for a simplified and easily usable condition assessment for land use management. Hynes (1960) clarified that water has lost its natural qualities only when pollution is severe enough to cause changes that exceed the boundaries of natural variation. Applying a large number of metrics to evaluate human impact can increase the sensitivity of this natural variation in biomonitoring models, but including metrics that do not discriminate between impaired and reference conditions will decrease sensitivity. Future experiments to refine metrics forensic of human impact will assist the biomonitoring model developed here and enable better quantification of variation from natural condition.
We evaluated the sensitivity of metrics for biological community types by performing traditional bioassessments for test sites with known human stressors as well as relationships of these metrics with various human activity gradients. For Biological Group 1, we identified that large reservoirs could impact ecosystem health, and that Log(n+1) total abundance, number of species, and NMDS Axis 1 may be promising metrics for the detection of reservoir impacts in this group (further description in Supplementary Material). Further, anthropogenic gradient analyses for Group 1 communities suggested that % detritivores or % Hyalella azteca metrics could be forensic indicators of landfill/lagoon/roads impacts (further description in Supplementary Material).
For Biological Group 2, we found that waste water and PAH contamination had significant impacts on the benthic macroinvertebrate community, but the presence of an upstream reservoir and urban center did not. The waste water effluent (likely due to high concentrations of ammonia) imparted large impacts in the ecosystem health, and the number of Coleoptera and community composition in NMDS Axis 1 provided promising responsiveness to waste water effluent impacts (further description in Supplementary Material). We also found metric responsiveness to PAH toxicity provided in an extreme example of oil contamination; the number of Coleoptera was reduced, and community composition in NMDS Axis 1 was altered relative to reference condition, similar to our results for ammonia pollution (further description in Supplementary Material). In contrast to waste water and PAH stressors, the single reservoir site in Biological Group 2 did not result in responsive metrics nor did the highest % urban cover we could locate in an upstream watershed, despite both possessing large amounts of co-occurring landfills, lagoons, and road crossings (further description in Supplementary Material). Ultimately, the metrics that show the most promise for detecting stress along a gradient for Biological Group 2 are total abundance, community composition (i.e., NMDS 1) and % EPT (further description in Supplementary Material).
Finally, we found no evidence of an ecosystem health impact of reservoirs in Biological Group 3 but did find land use relationships with some metrics of ecosystem health (further description in Supplementary Material). Stressor gradient analysis showed that the proportion of EPT can be an indicator of forest cover, and like Biological Group 2, % EPT increased as forest cover increased in the contributing watershed. Total abundance increased with increasing landfill/lagoons/roads and with increasing uncultivated as well as forest cover. Interpretation of how to apply this metric is complicated by its responsiveness to multiple variables, but combining it with a metric such as % EPT and/ or % Hyalella azteca may suggest an impact from landfill/lagoon/road stress if the latter two metrics are unaffected and total abundance is high relative to reference. As with Biological Group 2, the NMDS Axis 1 followed the similar gradient with uncultivated to forest cover.
Additional metrics for Group 3 demonstrated promise in detecting the landfill/lagoon/roads and cropland stressors. Specifically, % Hyalella azteca was positively related to the amount of cropland in a contributing watershed, but the number of Coleoptera were even more strongly positively related to the amount of uncultivated land in the upstream watershed and negatively related to forest cover. This may reflect a gradient of autochthonous to allochthonous inputs and autotrophic-based to heterotrophic-based food webs along a gradient from open pasture to forested watersheds as Hyalella azteca and Coleoptera are replaced in the benthic macroinvertebrate community by EPT taxa.
In addition to the strongest relationships with metrics, only Group 3 showed strong relationships with a number of single taxa. For example, haliplid beetles, Hyalella azteca, and elmid beetles all had negative correlations with increasing forest, while the caenid mayflies haliplid beetles and Hyalella azteca all produced positive correlations with increasing uncultivated open habitat. This would be an interesting hypothesis to pursue in the future to better understand the underpinning mechanism of these metrics and increase confidence in their forensic capabilities in the NGPs.
Regardless, the four metrics selected in each biological group are intended to be coarse summaries of the benthic macroinvertebrate community, and if they indicate disturbance relative to reference condition then their appropriate application would be to trigger further investigation on the part of water managers. In addition to more resolute water chemistry and habitat investigation would be the inclusion of additional metrics that may be more forensic of specific stressors such as the alternative metrics identified here, or emerging metrics being developed for the NGPs (e.g., Hoemsen 2015 for sediment or Collopy 2019 for PAH concentrations).

Implications of what reference means in the NGPs

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.

Acknowledgments

Benthic macroinvertebrates and habitat variable field collection from 2006 to 2009 were primarily performed by Kevin Kirkham, Deanne Schulz, Christine Markel, Alison Anton, Kirby Ebel, Amy Goodbrand and Lindsay Kosier. Janet Halpin conducted all sorting, and Dale Parker of AquaTax (Saskatoon, Saskatchewan) provided the bulk of the taxonomic identification. Assistance in mapping was provided by Dave MacDonald, acquisition of hydrology data by Curtis Hallborg, and interpretation of PAH data and comment on the manuscript provided by John-Mark Davies. Funding for the field work associated with this project was provided by the Water Security Agency of Saskatchewan’s Water Sustainability Fund, a grant from Saskatchewan Ministry of Agriculture, funding from the Fish and Wildlife Development Fund, and a grant from the Prairie Adaptation Research Collaborative (PARC). Support for manuscript preparation was provided by Natural Sciences and Engineering Research Council PGS-D and PDF scholarships to IDP.

References

Agriculture and Agri-Food Canada. 2005. Saskatchewan Soil Resource Database User’s Manual for SKIS. SKSISv2.pdf. Agriculture and Agri-Food Canada, Saskatoon Research Centre, Saskatchewan Land Resource Unit, Saskatoon, Saskatchewan, 43 pp.
Agriculture and Agri-Food Canada. 2008. 1:250,000 Quaternary Geology Map Series. [online]: Available from publications.gov.sk.ca/documents/310/97304-surficial_250k.htm
Agriculture and Agri-Food Canada. 2013. Effective Drainage Area of the AAFC Watersheds Project – 2013. [online]: Available from open.canada.ca/data/en/dataset/aeb7959a-9683-421a-be35-fb7c6100d0dc.
Agriculture Financial Services Corporation. 2006. 2006 Landcover 56m raster layer. [online]: Available from envgis.gov.sk.ca/arcgis/rest/services/Agriculture/AG_ACLMV2/MapServer/57.
Bailey RC, Kennedy MG, Dervish MZ, and Taylor RM. 1998. Biological assessment of freshwater ecosystems using a reference condition approach: comparing predicted and actual benthic invertebrate communities in Yukon streams. Freshwater Biology, 39: 765–774.
Barbour MT, Gerritsen J, Snyder BD, and Stribling JB. 1999. Rapid bioassessment protocols for use in wadeable streams and rivers: Periphyton, benthic invertebrates and fish. 2nd edition. EPA 841-B-99-002. Office of Water, US Environmental Protection Agency, Washington, DC.
Beatty JM, McDonald LE, Westcott FM, and Perrin CJ. 2006. Guidelines for Sampling Benthic Invertebrates in British Columbia Streams. British Columbia Ministry of Environment. 34 pp. [online]: Available from a100.gov.bc.ca/pub/eirs/finishDownloadDocument.do;jsessionid=EBCB5CA3593C4182940228356B2B8007?subdocumentId=6772
Boivin NL, Zeder MA, Fuller DQ, Crowther A, Larson G, Erlandson JM, Denham T, and Petraglia MD. 2016. Ecological consequences of human niche construction: Examining long-term anthropogenic shaping of global species distributions. Proceedings of the National Academy of Sciences, 113: 6388–6396.
Bouchard RW, Niemela S, Genet JA, Yoder CO, Sandberg J, Chirhart JW, et al. 2016. A novel approach for the development of tiered use biological criteria for rivers and streams in an ecologically diverse landscape. Environmental Monitoring and Assessment 188: 196.
Bowman MF, and Somers KM. 2005. Considerations when using the reference condition approach for bioassessment of freshwater ecosystems. Water Quality Research Journal of Canada, 40: 347–360.
Bowman MF, and Somers KM. 2006. Evaluating a novel test site analysis (TSA) bioassessment approach. Journal of the North American Benthological Society, 25: 712–727.
CCME (Canadian Council of Ministers of the Environment). 1999. Canadian Sediment Quality Guidelines for the Protection of Aquatic Life: Polycyclic Aromatic Hydrocarbons (PAHs). Canadian Council of Ministers for the Environment. [online]: Available from ceqg-rcqe.ccme.ca/download/en/243.
Clarke RT, Furse MT, Wright JF, and Moss D. 1996. Derivation of a biological quality index for river sites: Comparison of the observed with the expected fauna. Journal of Applied Statistics, 23: 311–332.
Clarke KR, and Warwick RM. 2001. Change in Marine Communities: an Approach to Statistical Analysis and Interpretation, 2nd edition. PRIMER-E, Plymouth, UK.
Clifford H. 1991. The Aquatic Macroinvertebrates of Alberta. University of Alberta Press, Edmonton, Alberta, Canada.
Collopy P. 2019. Developing an ecosystem health response fingerprint for freshwater oil contamination.M.SEM. thesis, School of Environment and Sustainability, University of Saskatchewan, Saskatoon, Saskatchewan. 33 pp.
Corkum LD. 1989. Patterns of benthic invertebrate assemblages in rivers of northwestern North America. Freshwater Biology, 21: 191–205.
Covich AP, Fritz SC, Lamb PJ, Marzolf RD, Mathews WJ, Poiani KA, et al. 1997. Potential effects of climate change on aquatic ecosystems of the Great Plains of North America. Hydrological Processes, 11: 993–1021.
Davies SP, and Jackson SK. 2006. The biological condition gradient: A descriptive model for interpreting change in aquatic ecosystems. Ecological Applications, 16: 1251–1266.
Dodds WK, Gido K, Whiles MR, Fritz KM, and Mathews WJ. 2004. Life on the edge: The ecology of Great Plains prairie streams. Bioscience, 54: 205–216.
Environment Canada. 2012. Canadian Aquatic Biomonitoring Network Field Manual. Freshwater Quality Monitoring and Surveillance – Atlantic Water Quality Monitoring and Surveillance Division, Science and Technology Branch, Environment Canada, Dartmouth, Nova Scotia. [online]: Available from ec.gc.ca/Publications/C183563B-CF3E-42E3-9A9E-F7CC856219E1/CABINFieldManual_EN_2012.pdf.
Evans-White M, Dodds WK, Huggins DG, and Baker DS. 2009. Thresholds in macroinvertebrate biodiversity and stoichiometry across water-quality gradients in Central Plains (USA) streams. Journal of the North American Benthological Society, 28: 855–868.
Faith DP, Minchin PR, and Belbin L. 1987. Compositional dissimilarity as a robust measure of ecological distance. Vegetatio, 69: 57–68.
Freeman PL, and Schorr MS. 2004. Influence of watershed urbanization on fine sediment and macroinvertebrate assemblage characteristics in Tennessee Ridge and Valley Streams. Journal of Freshwater Ecology, 19: 353–362.
Godwin RB, and Martin FRJ. 1975. Calculation of Gross and Effective Drainage Areas in the Prairie Provinces. In Canadian Hydrology Symposium – 1975 Proceedings, Winnipeg, Manitoba, 5 pp.
Government of Saskatchewan. 2019. Mining and Petroleum GeoAtlas, Oil and Gas GIS layer. [online]: Available from gisappl.saskatchewan.ca/Html5Ext/index.html?viewer=GeoAtlas.
Harrell FEJ. 2016. Hmisc: Harrell Miscellaneous. [online]: Available from CRAN.R-project.org/package=Hmisc.
Hawkins CP, Norris RH, Hogue JN, and Feminella JM. 2000. Development and evaluation of predictive models for measuring the biological integrity of streams. Ecological Applications 10: 1456–1477.
Hilsenhoff WL. 1988. Rapid field assessment of organic pollution with a family-level biotic index. Journal of the North American Benthological Society, 7: 65–68.
Hoemsen BM. 2015. Indicator invertebrates: Determining change in benthic macroinvertebrate communities due to deposited sediment in the Northern Great Plains. M.Sc. Thesis, University of Saskatchewan, Saskatoon, Saskatchewan.
Hynes HBN. 1960. The Biology of Polluted Water. Liverpool University Press, Liverpool, UK.
Karr JR, Fausch KD, Angermeyer L, Yant PR, and Schlosser IJ. 1986. Assessment of biological integrity in running waters: A method and its rationale. Illinois Natural History Survey Special Publication No. 5, Illinois Natural History Survey, Urbana Champaign, IL.
Kilgour BW, Somers KM, and Mathews DE. 1998. Using the normal range as a criterion for ecological significance in environmental monitoring and assessment. Ecoscience, 5: 542–550.
Klemm JJ, Blocksom KA, Fulk FA, Herlihy AT, Hughes RM, Kaufmann PR, et al. 2003. Development and evaluation of a macroinvertebrate biotic integrity index (MBII) for regionally assessing Mid-Atlantic Highland streams. Environmental Management, 31: 656–669.
Kolkowitz R, and Marsson M. 1908. Ökologie der pflanzlichen Saprobien. Berichte der Deutschen Botanischen Gesellschaft, 26A: 505–519.
Linke S, Norris RH, Faith DP, and Stockwell D. 2005. ANNA: A new prediction method for bioassessment programs. Freshwater Biology, 50: 147–158.
Mason PG, and Parker DW. 1994. Additions and corrections to the list of non-biting midges in Saskatchewan. Blue Jay, 54: 200–203.
McCormick FH, Hughes RM, Kaufmann PR, Peck DV, Stoddard JL, and Herlihy AT. 2001. Development of an index of biotic integrity for the Mid-Atlantic Highlands Region. Transactions of the American Fisheries Society, 130: 857–877.
McCune B, and Grace JB. 2002. Analysis of ecological communities. MjM Software Design, Gleneden Beach, Oregon.
Merritt RW, Cummins KW, and Berg MB. (Editors). 2008. An introduction to the aquatic insects of North America. 4th edition. Kendall Hunt Publishing Company, Dubuque, Iowa.
National Ecological Framework for Canada. 2006. Ecoregions GIS Dataset. [online]: Available from sis.agr.gc.ca/cansis/nsdb/ecostrat/gis_data.html.
Natural Resources Canada. 2016. Canadian Digital Elevation Model. [online]: Available from open.canada.ca/data/en/dataset/7f245e4d-76c2-4caa-951a-45d1d2051333.
Natural Resources Canada. 2018. National Topographic System Maps, 1:50,000. [online]: Available from open.canada.ca/data/en/dataset/28a599de-59d5-50e8-8539-f4566f7982b6.
National Resources Canada. 2019. National Road Network shapefiles. [online]: Available from open.canada.ca/data/en/dataset/3d282116-e556-400c-9306-ca1a3cada77f.
Newbury RW. 1984. Hydrologic determinants of aquatic insect habitats. In The ecology of aquatic insects. Edited by VH Resh and DM Rosenberg. Praeger, New York. pp. 323–357.
Ode PR, Rehn AC, and May JT. 2005. A quantitative tool for assessing the integrity of Southern California coastal streams. Environmental Management, 35: 493–504.
Ode PR, Rehn AC, Mazor RD, Schiff KC, Stein ED, May JT, et al. 2016. Evaluating the adequacy of a reference-site pool for ecological assessments in environmentally complex regions. Freshwater Science, 35: 237–248.
Papworth SK, Rist J, Coad L, and Milner-Gulland EJ. 2009. Evidence for shifting baseline syndrome in conservation. Conservation Letters, 2: 93–100.
Parker D, and Glozier N. 2005. First record of the non-biting midge Zavreliella marmorata (Wulp.) (Chironomidae: Diptera), from Saskatchewan. Blue Jay, 63: 200–202.
PFRA (Prairie Farm Rehabilitation Association) Hydrology Division. 1983. The Determination of Gross and Effective Drainage Areas in the Prairie Provinces. Hydrology Report #104, Agriculture Canada, PFRA Engineering Branch, Regina, Saskatchewan, 22 pp.
Phillips ID, Parker D, Hoemsen BM, Bell AJ, and Chivers DP. 2013. Biological notes and range expansion of the non-biting midge Odontomesa fulva (Kieffer) (Diptera: Chironomidae). Western North American Naturalist, 73: 244–247.
Phillips ID, Pollock MS, Bowman MF, and Chivers DP. 2015. Thermal alteration and macroinvertebrate response below a large Northern Great Plains reservoir. Journal of Great Lakes Research 41: 155–163.
Phillips ID, Davies J-M, Bowman MF, and Chivers DP. 2016. Macroinvertebrate communities in a Northern Great Plains river are strongly shaped by naturally occurring suspended sediments: implications for ecosystem health assessment. Freshwater Science, 35: 1354–1364.
Podani J. 2000. Introduction to the Exploration of Multivariate Biological Data. Backhuys Publishers, Leigen.
Poff LL, Allan JD, Bain MB, Karr JR, Prestegaard KL, Richter BD, et al. 1997. The natural flow regime. Bioscience, 47: 769–784.
Qian SS, King RS, and Richardson CJ. 2003. Two statistical methods for the detection of environmental thresholds. Ecological Modelling, 166: 87–97.
R Development Core Team. 2013. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0, [online]: Available from R-project.org.
Resh VH, and Jackson JK. 1993. Rapid assessment approaches to biomonitoring using benthic macroinvertebrates. In Freshwater biomonitoring and benthic macroinvertebrates. Edited by DM Rosenberg and VH Resh. Chapman and Hall, New York. pp. 195–223.
Resh VH, Rosenberg DM, and Reynoldson TB. 2000. Establishing reference conditions in the Fraser River catchment, British Columbia, Canada, using the BEAST (BEnthic Assessment of SedimenT) predictive model. In Assessing the biological quality of fresh waters. RIVPACS and other techniques. Edited by JF Wright, DM Sutcliffe, and MT Furse. Freshwater Biological Association, Ambleside, UK. pp. 181–194.
Reynoldson TB, Bailey RC, Day KE, and Norris RH. 1995. Biological guidelines for the freshwater sediments based on BEnthic Assessment of SedimenT (the BEAST) using a multivariate approach for predicting biological state. Australian Journal of Ecology, 20: 198–219.
Reynoldson TB, and Wright JF. 2000. The reference condition: Problems and solutions. In Assessing the Biological Quality of Fresh Waters. RIVPACS and other Techniques. Edited by JF Wright, DW Sutcliffe, and MT Furse. Freshwater Biological Association, Ambleside, UK. pp. 293–303.
Reynoldson TB, Rosenberg DM, and Resh VH. 2001. Comparison of models predicting invertebrate assemblages for biomonitoring in the Fraser River catchment, British Columbia. Canadian Journal of Fisheries and Aquatic Sciences 58: 1395–1410.
Ripley B, Venables B, Bates DM, Hornik K, Gebhardt A, Firth D, and Ripley MB. 2013. Package ‘mass’. Cran r. 538: 133–120. [online]: Available from stats.ox.ac.uk/pub/MASS4/
Roy AH, Rosemond AD, Paul MJ, Leigh DS, and Wallace JB. 2003. Stream macroinvertebrate response to catchment urbanization (Georgia, U.S.A.). Freshwater Biology, 48: 329–346.
SaskH2O. 2019. Landfills data. [online]: Available from saskh2o.ca/landfills.asp.
Saskatchewan Ministry of Environment. 2006. Sde.water_stream_network_v2. Based on National Topographic Series 1:50,000 water data.7.
Stoddard JL, Larsen DP, Hawkins CP, Johnson RK, and Norris RH. 2006. Setting expectations for the ecological condition of streams: The concept of reference condition. Ecological Applications 16: 1267–1276.
Strahler AN. 1964. Quantitative geomorphology of drainage basins and channel networks. In Handbook of Applied Hydrology. 4th edition. Edited by VT Chow. McGraw-Hill, New York. pp. 39–76.
Statistics Canada. 2016. 2016 Census – Boundary Files. [online]: Available from www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2016-eng.cfm.
Vörösmarty CJ, McIntyre PB, Gessner MO, Dudgeon D, Prusevich A, Green P, et al. 2010. Global threats to human water security and river biodiversity. Nature, 467: 555–561.
Waters CN, Salasiewicz J, Summerhayes C, Barnosky AD, Poirier C, Galuszka A, et al. 2016. The Anthropocene is functionally and stratigraphically distinct from the Holocene. Science, 351(6269): 137.
Webb JM. 2002. The mayflies of Saskatchewan. M.Sc. thesis, University of Saskatchewan, Saskatoon, Saskatchewan, Canada.
Webb JM, Parker DW, Lehmkuhl DM, and McCafferty WP. 2004. Additions and emendations to the mayfly (Ephemeroptera) fauna of Saskatchewan, Canada. Entomological News, 115: 213–218.
Williams WD. 1985. Biotic adaptations in temporary lentic waters, with special reference to those in semi-arid and arid regions. Hydrobiologia, 125: 85–110.
World Water Assessment Programme. 2009. United Nations World Water Development Report 3: Water in a Changing World.
Wright JF. 1984. A preliminary classification of running-water sites in Great Britain based on macroinvertebrate species and the prediction of community type using environmental data. Freshwater Biology, 14: 221–256.
Wright JF. 1995. Development and use of a system for predicting the macroinvertebrate fauna in flowing waters. Australian Journal of Ecology, 20: 181–197.

Supplementary material

Supplementary Material 1

Information & Authors

Information

Published In

cover image FACETS
FACETS
Volume 8Number 1January 2023
Pages: 1 - 31
Editor: Irene Gregory-Eaves

History

Received: 12 July 2022
Accepted: 28 December 2022
Version of record online: 30 March 2023

Data Availability Statement

Data generated or analyzed during this study are available in the Dryad repository, https://doi.org/10.5061/dryad.4f4qrfjgt.

Key Words

  1. macroinvertebrates
  2. stream bioassessment
  3. Northern Great Plains
  4. reference condition approach
  5. Test Site Analysis

Sections

Subjects

Plain Language Summary

A new tool for aquatic ecosystem health assessment in the Canadian Prairies

Authors

Affiliations

Iain D. Phillips [email protected]
Ecological and Habitat Assessment, Water Security Agency, 10-3904 Millar Avenue, Saskatoon, SK S7P 0B1, Canada
Department of Biology, University of Saskatchewan, 112 Science Place, Saskatoon, SK, S7N 5E2, Canada
Glen McMaster
2275 Retallack Street, Regina, SK S4T 2K7, Canada
Douglas P. Chivers
Department of Biology, University of Saskatchewan, 112 Science Place, Saskatoon, SK, S7N 5E2, Canada
Michelle F. Bowman
Forensecology, 70 Swift Crescent, Guelph, ON N1E 7J1, Canada

Author Contributions

All conceived and designed the study.
IDP performed the experiments/collected the data.
IDP analyzed and interpreted the data.
IDP, GM, and MFB contributed resources.
All drafted or revised the manuscript.

Competing Interests

The authors declare there are no competing interests.

Metrics & Citations

Metrics

Other Metrics

Citations

Cite As

Export Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

There are no citations for this item

View Options

View options

PDF

View PDF

Get Access

Media

Media

Other

Tables

Share Options

Share

Share the article link

Share on social media