Multiple stressor effects on stream health in the Lake Simcoe Watershed

JGLR-01095; No. of pages: 12; 4C: Journal of Great Lakes Research xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

Journal of Great Lakes Research journal homepage: www.elsevier.com/locate/jglr

Multiple stressor effects on stream health in the Lake Simcoe Watershed Noreen E. Kelly a,⁎,1, Eavan M. O'Connor b, Robert F. Wilson b, Joelle D. Young c, Jennifer G. Winter c, Lewis A. Molot a a b c

Faculty of Environmental Studies, York University, 4700 Keele St., Toronto, ON, Canada M3J 1P3 Lake Simcoe Region Conservation Authority, 120 Bayview Parkway, Newmarket, Ontario, Canada L3Y 3W3 Ontario Ministry of the Environment and Climate Change, 125 Resources Rd., Toronto, ON, Canada M9P 3V6

a r t i c l e

i n f o

Article history: Received 24 January 2016 Accepted 25 June 2016 Available online xxxx Communicated by Joseph Makarewicz Index words: Benthic invertebrates Fish Hydrology Temperature Water quality Multiple stressors

a b s t r a c t The Lake Simcoe Watershed (LSW) is an important natural resource, supporting varied agricultural, recreational, and tourism activities for millions of people. Historically, major alterations to land uses have occurred throughout the watershed, resulting in increases of nutrients and contaminants, and alteration of hydrological and thermal regimes. The combination of these stressors has the potential to elicit greater impacts on ecosystem health than any one stressor in isolation, yet no studies have examined their effects on stream health and aquatic biota in the LSW. In this study, we quantified the impacts of multiple stressors on indices of biotic integrity, as well as assessments of benthic invertebrate and fish community composition from 2004 to 2012. Using a suite of multivariate analyses, we examined stressors across three categories of environmental variables: water quality, temperature, and hydrology. Water quality explained the largest amount of variation in stream health and biological community composition, followed by water temperature and hydrological variables, respectively. Total phosphorus, dissolved oxygen, and iron concentrations, as well as the 7-day low flow and the mean summer temperature, were identified as predictors of both fish and benthic invertebrate variance. Multiple stressor interactions were detected for fish communities, although these interactions were of lower relative importance than for individual stressors. In contrast, no significant stressor interactions were found to influence benthic invertebrate communities. Overall, our results reinforce the importance of ecosystem monitoring, and the need to consider the influence of multiple aspects of ecosystem function, when examining the health of aquatic systems. © 2016 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Introduction The cumulative effect of a variety of human activities profoundly alters the chemical, physical, and biological processes occurring in streams across varying spatial and temporal scales and thereby impacts the ecological integrity of river ecosystems. Land-use alteration is a major driver of anthropogenic change in river systems, resulting in a complex cascade of changes that impacts water quality, flow regime, habitat availability, and ultimately aquatic biota (Allan, 2004). For example, the conversion of natural areas for agricultural or urban land uses often results in increases of nonpoint sources of nutrients and pollutants, degradation of riparian and stream channel habitat, and altered flow regimes (Allan, 2004). Eutrophication, caused by excessive inputs of phosphorus and nitrogen, leads to overgrowth of algae and macrophytes that alters the distribution of dissolved oxygen, impacting ⁎ Corresponding author. E-mail address: [email protected] (N.E. Kelly). 1 Current address: Department of Physical and Environmental Sciences, University of Toronto Scarborough, Toronto, ON, Canada M1C 1A4.

habitat availability and aquatic biodiversity (Carpenter et al., 1998). Modified flow regimes can influence the frequency and spatial extent of physical disturbances in a stream ecosystem and can alter the transport and distribution of nutrients and sediments, which in turn directly affects species abundance, composition, and diversity of aquatic communities (Bunn and Arthington, 2002). The loss of forest cover and riparian vegetation to accommodate urban and agricultural activities is associated with increased bank and channel erosion (Gregory et al., 1991), as well as alterations to the thermal regime of streams, which in turn influences primary production and the physiology and distribution of aquatic species (Bourque and Pomeroy, 2001; Poole and Berman, 2001). Most landscape-scale studies have contrasted the varying extent of agricultural, urban, and forested land to examine the impacts of human-dominated land use on stream health and community composition. Such studies have generally observed shifts to pollution-tolerant taxa or warm water-tolerant species, decreases in diversity and overall abundance of stream biota, and changes to aquatic food webs (e.g., Bazinet et al., 2010; Lenat and Crawford, 1994; Moore and Palmer, 2005; Stanfield and Kilgour, 2006, 2012). Collectively, these

http://dx.doi.org/10.1016/j.jglr.2016.07.007 0380-1330/© 2016 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

2

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

studies have demonstrated the importance of the surrounding landscape and human activities on the ecological integrity of streams across multiple spatial scales. Multiple environmental stressors are now a common problem facing many freshwater ecosystems around the globe (Ormerod et al., 2010). Interactions among multiple stressors can result in “ecological surprises” that can generate complex effects that lessen or amplify the individual impacts of each stressor (Christensen et al., 2006; Folt et al., 1999). Evidence is mounting that multiple stressors, acting at different spatial and temporal scales, can interact to affect water quality, biodiversity, and ecosystem function (e.g., Merovich and Petty, 2007; Piggott et al., 2012; Townsend et al., 2008). For example, in a mined Appalachian watershed in the USA, Merovich and Petty (2007) attributed extensive ecological loss and poor benthic macroinvertebrate survival to the combined effects of acid mine drainage and heat inputs, compared to either stressor in isolation. Piggott et al. (2012) manipulated nutrients, sediment, and riparian vegetation loss in agricultural streamside channels in New Zealand. Their experiments demonstrated that interactions between sediment concentrations and stream temperature altered algal and benthic invertebrate richness, with elevated temperatures strengthening the negative impact of added sediment, further reducing biodiversity. By combining survey and experimental approaches in grassland streams in New Zealand, Townsend et al. (2008) provided evidence for the antagonistic impacts of high nutrient concentrations and fine sediment cover on benthic invertebrate populations and diversity, with the largest reductions in these metrics occurring when high nutrient and sediment concentrations coincided. However, an enhanced understanding of where the combined action of multiple stressors produces complex versus simple responses, across a wide array of stream habitats and trophic levels, is still needed (Townsend et al., 2008). In addition, the net effects of multiple stressors pose significant challenges for the management of stream ecosystems, as land-use change typically alters more than one ecosystem component simultaneously. Thus, scientific studies that can directly identify anthropogenic stressors, and their subsequent impacts and interactions, provide valuable information to assist management actions aimed at ameliorating human-use impacts. The Lake Simcoe Watershed (LSW), in south-central Ontario, Canada, is an important natural resource, supporting varied agricultural, recreational, and tourism activities for millions of people (LSRCA, 2013a). Collectively, the ecological goods and services provided by the watershed are valued between $700 and $975 million annually (LSEMS, 2008; Wilson, 2008). Since European settlement in the 1800s, major changes have occurred within the LSW, as land was cleared for agriculture, roads, and urban and industrial development. Dramatic environmental changes, such as the destruction of different types of natural habitats; increased nutrient inputs to the watershed through increased erosion, fertilization, and waste disposal; and increased natural resource utilization (e.g., gravel extraction, firewood harvest, fishery pressure), have directly affected the quality of the LSW ecosystem (Evans et al., 1996). While the percentage of agricultural land has decreased slightly in recent years (Oni et al., 2014), the percentage of urban areas has increased (O'Connor et al., 2013), which could continue, as the population is projected to grow by ~30% by 2031 (MOI, 2006). Since the 1970s, excessive phosphorus concentrations have been identified as one of the primary drivers of anthropogenic stress throughout the LSW, and reductions of TP loads from tributaries have been a key management goal for several decades (LSEMS, 2008). Consequently, much of the recent research in the LSW has been aimed at quantifying the spatial and temporal dynamics of nutrient concentrations and exports to the lake (Baulch et al., 2013; Miles et al., 2013; O'Connor et al., 2011; Oni et al., 2014), as well as investigating their relationships with wetland, agricultural, and urban land uses (North et al., 2013; Oni et al., 2015), or understanding controls on nutrient dynamics across multiple subcatchments through modelling applications (Crossman et al., 2013; Oni et al., 2011; Whitehead et al., 2011).

Collectively, these empirical and modelling studies suggest that combinations of multiple stressors (e.g., changes in nutrient concentrations, stream flow, and stream temperature that are consequences of changes in land use and climate) are now prevalent in portions of the LSW and have the potential to impart greater impacts on ecosystem health than any one of these single stressors in isolation. Yet, to date, no study has explicitly quantified either the individual or the combined impacts of these potential stressors on stream health and aquatic biota in the LSW. The LSW is an ideal setting for examining the impacts of multiple stressors on tributary ecological health, due to widespread gradients in environmental characteristics (e.g., land use, water chemistry, hydrology, temperature, biology) among catchments. As a result of expanded environmental monitoring over the last decade in multiple areas of the watershed by multiple agencies (LSRCA, 2013a; Young and La Rose, 2014), a large quantitative dataset containing physical, chemical, hydrological, and biological parameters is presently available, providing the means for the analysis of multiple stressor impacts. In this study, we employ a suite of multivariate analyses to identify the parameters influencing the ecological health of Lake Simcoe tributaries. Specifically, we examine the impacts of multiple stressors across these three categories on stream biological condition and community composition. We focus on data collected from six major catchments (or subcatchments), which are spatially distributed throughout the LSW which represent gradients in water quality, temperature, and hydrology while being differentially impacted by human activity. This study provides the first quantitative assessment of the impacts of multiple environmental stressors on aquatic biota across the LSW. The findings presented here will inform management actions in the region, as well as contribute to a broader understanding of multiple stressors as drivers of ecological change in stream ecosystems.

Methods Study area and time period The LSW has a total land area of 2899 km2, which is drained by 35 tributaries within 18 catchments. As of 2008–2009, agriculture was the dominant land use within the LSW, as 45% of the entire watershed was used for hay and pasture, cropland, and turf and sod operations, while 7% was urban development and roads, and 35% was natural vegetative cover (e.g., woodlands, wetlands, and riparian areas) (LSRCA, 2013a; MOECC, 2015). A subset of sampling sites and years was selected from all available data, which comprised sampling sites that had all biological and environmental variables sampled at approximately the same geographical location within a catchment over the same period of record. This left us with a dataset containing information from one long-term sampling station located within each of six catchments (from the total 18) (Fig. 1) for a period of 2004–2012 for fish, and 2004–2007 and 2009–2011 for benthic invertebrates. The Beaver, Black, Pefferlaw, and Upper Schomberg River catchments are located in the southern portion of the LSW, represent 40% of the total area draining into the lake, and collectively have a large influence on the quality of Lake Simcoe's aquatic ecosystem (LSRCA, 2013a). These catchments are located within the most intensively farmed portion of the watershed (LSRCA, 2013a). In contrast, the Lovers Creek and Hawkestone Creek catchments are located in the west and northwest portion of the watershed, respectively. The Lovers Creek catchment contains a higher degree of urbanization (LSRCA, 2012a), while the Hawkestone Creek catchment contains one of the highest levels of natural vegetative cover (e.g., woodlands, wetlands, and riparian areas) in the LSW (LSRCA, 2013b). Collectively, these six catchments vary considerably in their land use, hydrological regime, thermal regime, and water quality (Table 1). Further information on the catchments used in this study can be found in LSRCA (2010a, 2010b, 2012a, 2012b, 2012c, 2013b).

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

3

Fig. 1. Map of Lake Simcoe and its watershed. Inset shows the location of Lake Simcoe in relation to the Laurentian Great Lakes. Squares are catchment sampling stations; circles are locations of major urban centres.

Simcoe Protection Plan) program and the PWQMN (Provincial Water Quality Monitoring Network) program. Under the LSPP program, water samples were collected year round, every 2 weeks during the spring, summer, and fall and every 3 weeks in the winter months, and

Sample collection We analyzed selected chemical parameters from stream water samples collected as part of two monitoring programs: the LSPP (Lake

Table 1 Land uses (% of total catchment area) for six Lake Simcoe catchments examined in this study. Catchment

Beaver Black Hawkestone Lovers Pefferlaw Upper Schomberg

Abbreviation

BV BL HS LV PF US

Catchment area (km2)

327 375 48 60 425 44

Land use (%) AGR

FOR

TRA

URB

WET

64.8 39.8 36.1 35.7 49.5 58.9

4.3 13.1 23.0 13.7 11.5 11.7

1.3 2.4 0.0 9.7 1.6 1.7

2.3 4.7 3.1 16.3 5.7 3.8

16.8 21.1 20.9 13.7 12.5 4.6

Q (m3/s)

T (°C)

TP (μg/L)

TN (mg/L)

Cl (mg/L)

3.45 2.76 0.55 0.75 3.97 0.35

21.4 22.5 20.2 19.5 21.7 20.1

36.8 22.0 11.5 20.2 28.8 59.5

1.07 1.27 0.94 1.91 1.18 0.93

40.9 30.0 21.5 109.1 28.6 49.3

Major categories of land-use data include agriculture (AGR; row crop, sod farm, hay, pasture), forest (FOR; coniferous, deciduous, mixed), transport (TRA; roads, rail), urban (URB; commercial, estate, industrial, golf course), and wetlands (WET; bog, fen, swamp, marsh). Streamflow (Q), summer temperature (T), and total phosphorus (TP), total nitrogen (TN), and chloride (Cl) concentrations represent mean values for the 2004–2012 period of record. Data provided by the LSRCA (see also LSRCA, 2005; Scott et al., 2006).

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

4

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

supplemented with storm/rain event sampling. The PWQMN samples were collected 8 times a year on a monthly basis during the ice-free period. Datasets from each sampling program (LSPP and PWQMN) were combined (averaged) when sampled on the same date. Chemical parameters from both programs were analyzed by the Laboratory Services Branch of the Ministry of Environment and Climate Change following standard methods (OMOE, 1983) and included total phosphorus (TP), total ammonia (NHT), total nitrates (NO2 + NO3), total Kjeldahl nitrogen (TKN), and total nitrogen (TN = total nitrates + TKN), calcium (Ca), chloride (Cl), conductivity (Cond), total suspended solids (TSS), pH, and a suite of metals (Al, Ba, Cu, Fe, Mn, Sr, Ti, V, Zn; collected from the PWQMN sampling program only). Dissolved oxygen (DO) was measured instantaneously at each site using a YSI sonde. For further details, see LSRCA (2013a). Temperature was measured in Lake Simcoe tributaries using Hobo temperature loggers installed in shallow (wading depth) streams ~ 10 m downstream from fish monitoring stations. Temperature was measured on an hourly basis from May to October following the Ontario Stream Assessment Protocol manual (Stanfield, 2010). For further details, see LSRCA (2013a). We used hydrological data collected from hydrometric gauges maintained by the Water Survey of Canada (wateroffice.ec.gc.ca) for the Black, Beaver, Hawkestone, Pefferlaw, and Upper Schomberg River catchments and by the LSRCA for Lovers Creek catchment. At each gauge, water level (stage) was recorded continuously using either a float recorder or constant flow bubbler and data logger combination. The continuous stage record was converted to continuous discharge (volume of water per unit time) using stage–discharge relationships developed for each site (LSRCA, 2013a). For continuity with other measured tributary variables (e.g., water quality, biology), the discharge time series from three tributaries (Black, Beaver, and Pefferlaw Rivers) were prorated to the location of water quality and biology sampling locations (O'Connor et al., 2009). Benthic invertebrates were sampled in the fall between September 1 and November 30. Sampling of benthic invertebrate communities followed the transect-travelling kick and sweep method from the Ontario Benthos Biomonitoring Network protocol (Jones et al., 2004); exceptions were sampling one riffle instead of two and enumerating 300 (instead of 100) individuals (to highest taxonomic resolution possible, typically genera or family level). Ten meters of each riffle was kicked for 3 min so that effort could be measured and standardized. Fish were sampled between June 1 and September 30 (with the majority of sampling occurring in July and August) primarily using backpack electrofishers, to sample all habitats within at least one pool riffle sequence following the Ontario Stream Assessment Protocol (Stanfield, 2010). Each sampling run was at least 40 m long. Shocker seconds, sampling time, wetted widths, and channel lengths were recorded so that effort could be assessed and standardized. Blocker nets were not utilized for this study. Fish were identified to species level and enumerated on site and released alive; some specimens were retained for confirmation of identification. For further sampling details, see LSRCA (2013a). Data analysis Our main objective was to examine the impact of multiple environmental stressors on the aquatic biota of Lake Simcoe tributaries. To do this, we performed four separate multivariate analyses on the benthic invertebrate dataset and the fish dataset: linear discriminant analysis (LDA), multiple linear regressions, and redundancy (RDA) and variance partitioning (VPA) analysis. The benthic invertebrate or fish response metrics for the LDAs and multiple linear regressions were the rankings and scores of biological health indices, respectively, and community composition for the RDAs and VPAs. All analyses were conducted using the statistical software R (R Core Team, 2015). The index of stream health for benthic invertebrates was the modified family-level Hilsenhoff Index of Biotic Integrity (HBI) (Hilsenhoff, 1988; LSRCA,

2013a). The HBI estimates the overall tolerance of the benthic invertebrate assemblage toward organic (nutrient) enrichment, with higher scores (e.g., N 5) indicating a higher abundance of tolerant taxa, suggesting a higher level of pollution. For fish, we used an Index of Biotic Integrity (IBI) (Karr, 1991), modified by the Toronto Region Conservation Authority to be relevant for Oak Ridges Moraine watercourses (OMNR and TRCA, 2005). The IBI is a multi-metric measure of stream quality that uses fish fauna as biological indicators and is used to rate the overall health of a stream on a scale of 9 (poor) to 45 (very good). To evaluate fish and benthic invertebrate community composition, we created a community matrix based on relative abundances of fish or benthic invertebrates, identified to species or families (e.g., lowest possible taxonomic resolution), respectively. To remove the impact of rare species in subsequent analyses, we reduced the community matrices following the rank-abundance plotting procedure outlined in Kelly et al. (2013), retaining those taxa found at ≥ 1% relative abundance and ≥ 5% occurrence across all sites. Taxon abundances were log (x + 1) transformed prior to analysis to improve normality and lessen the influence of dominant taxa. The independent variables in the multivariate analyses came from three categories of environmental stressors: water quality, temperature, and hydrology. For water quality stressors, we used the annual median concentrations of water quality parameters and metals. We adjusted the total nitrates, total ammonia, and metals datasets to account for nondetectable values (Helsel, 2012, chap. 12) using the R package NADA (Lee, 2012). Summertime stream temperature metrics were calculated from a time-series of mean daily temperature (°C) for the summer period (June 21 to September 20) and included the mean (MST), standard deviation (SDST), and coefficient of variation (CVST) of temperature, as well as the mean daily maximum (MDMax) and mean daily minimum (MDMin) temperature. For the hydrological variables, we used daily average discharge (m3/s) time series to calculate four different streamflow statistics that represent the storm and base flow patterns over annual timescales that can influence the biological conditions of streams (Konrad and Booth, 2002; Konrad et al., 2005): (1) Q mean, the annual mean streamflow; (2) Q max, the maximum annual streamflow; (3) Q min, 7-day low flow, or the minimum mean daily discharge for 7 consecutive days; and (4) TQmean, the fraction of a year that streamflow exceeds the mean annual streamflow (e.g., Q daily/Q mean N 1). Q mean is an indicator of the annual runoff volume in a stream, while TQmean provides a measure of the distribution of annual runoff between stormflow and recessional and baseflow periods. Q max and Q min are indicators of ecological disturbance, providing measures of the magnitude of the largest annual flood and the severity of the annual drought, respectively (Konrad and Booth, 2002). Further, these hydrological statistics are related to the level of urban development and thus represent mechanisms by which urbanization can influence biological health in a stream basin (Konrad and Booth, 2002; Konrad et al., 2005). All metrics were calculated for hydrological years beginning June 1 of the current year to May 31 of the following year. As several years of flow data were missing for Hawkestone Creek, flow estimates were made using the relationship between measured flows at that station and another gauge station (index station) in the LSW (O'Connor et al., 2009). Linear discriminant analysis (LDA) (also known as discriminant function analysis) was used to assess the relationship between HBI or IBI rankings (e.g., poor, good, very good) and environmental variables (Joy and Death, 2004). LDA is a multivariate procedure that generates functions of variables that best distinguish among previously defined groups of samples by maximizing the separation among those groups (Borcard et al., 2011). We used LDA (function lda in R package MASS; Venables and Ripley, 2002) to find discriminant functions of the environmental variables that best distinguished among the HBI or IBI rankings. For the available data, there were six different HBI rankings (excellent, very good, good, fair, fairly poor, poor) and three IBI rankings (good, fair, poor). Because many of the physicochemical variables within each environmental predictor category were highly correlated, we

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

used principal component analysis (PCA) to identify the primary gradients of change and reduced the number of input variables from each category for subsequent analyses. Variables selected after the PCA screening included Cl, DO, pH, TN, TP, Fe (water quality), CVST, MST (temperature), Q mean, Q max, and Q min (hydrology). The final LDA model was checked for homogeneity of multivariate dispersion (betadisper function in R package vegan; Oksanen et al., 2015), and MANOVA (Wilk's test) was used to determine the significance of environmental predictors (at α = 0.05) (Borcard et al., 2011). Stepwise multiple linear regressions were used to assess multiple stressor effects on HBI or IBI scores. To further reduce the number of environmental predictors, we used the regsubsets function in the R package leaps (Lumley, 2009) to choose a maximum of three variables from each predictor category that contributed the most to the R2adj value of a multiple regression with the variables from a single predictor category. All significant terms from each predictor category were then input into an overall regression model. Predictor variables were meancentered to reduce multicollinearity. Backwards stepwise selection was used to determine the most parsimonious model, based on minimizing the Akaike's information criterion (AIC). In addition, individual predictor variables were retained in the final models only if significant at p ≤ 0.05 and variance inflation factors were b 10. Model residuals were examined for normality, homoscedasticity, and independence; outliers were removed as necessary to meet linear regression requirements. The relative importance of all retained variables (i.e., the contribution of each predictor to the overall variance) was calculated in the final models (function calc.relimp in R package relaimpo; Grömping, 2006). Regression analyses were first performed with only the main effects and subsequently with both the main effects and their first-order interactions. Analysis of variance (ANOVA) between the linear models was used to assess whether interactions improved model fit. The final regression models were then interpreted to identify the presence of additive and multiplicative (synergistic or antagonistic) effects (Thrush et al., 2008). Additive effects were identified when variables were considered independent of others in the model, while multiplicative effects were identified by the presence of significant interaction terms. The sign of the parameter estimate for the interaction terms indicated whether synergistic or antagonistic effects (i.e., positive vs. negative parameter estimates) worked to increase or decrease the main effect identified by the model. Variance partitioning analysis (VPA) (Borcard et al., 1992) was used to explore the relationship between the three categories of environmental variables (water quality, temperature, hydrology) and benthic invertebrate or fish community composition. Prior to analysis, the benthic invertebrate (family level) and fish (species level) assemblage matrices were Hellinger-transformed (Legendre and Gallagher, 2001), and the environmental predictors mean-centered to reduce multicollinearity. Given the number of highly correlated explanatory variables in each predictor category and that the percent variance explained by each category is sensitive to the number of variables included in the analysis (Borcard et al., 1992), the complexity of the dataset was further reduced (following Paterson et al., 2008). Considering each predictor category separately, (1) we conducted constrained ordinations (RDA) to identify individual variables that explained significant amounts of the benthic invertebrate or fish variance; (2) we applied a forward selection procedure (function forward.sel in the R package packfor; Dray et al., 2013) to each RDA, using the double-stopping criteria of Blanchet et al. (2008); and (3) using the significant variables selected in step 2, we conducted a subsequent RDA and used permutational ANOVAs (function anova.cca in the R package vegan) to check the significance of each RDA model and the individual predictor variables. Significant variables were removed if their variance inflation factors were N 10. Pre-screening reduced the number of significant explanatory variables from two to three depending on the predictor category. A three-category VPA was then performed for each of the benthic invertebrate and fish assemblages in order to quantify the variation that can be

5

explained uniquely by each predictor category and by co-variation among categories (using the varpart function in the R package vegan). Permutational ANOVA was used to examine the significance of all testable fractions of variance in the VPA. Results Benthic invertebrates and multiple stressors LDA was used to discriminate among sites with scores in the six quantitative HBI rankings (Fig. 2). For HBI rankings, two LDA axes captured ~ 80% of the variation (LDA 1 = 64.2% of variance; LDA 2 = 17.1% of variance) among catchments (MANOVA Wilk's test: W5,36 = 0.031, approx F55,134 = 2.51, p b 0.0001). LDA 1 separated sites with excellent and very good rankings from the fair, fairly poor, and poor rankings, while LDA 2 separated the good rankings (Fig. 2). Five significant environmental predictors discriminated the sites: TP (F5,36 = 12.1, p b 0.0001), Fe (F5,36 = 6.68, p = 0.0002), TN (F5,36 = 3.91, p = 0.006), DO (F5,36 = 3.22, p = 0.017), and CVST (F5,36 = 4.17, p = 0.004). Sites with fair, fairly poor, and poor rankings were associated with high TP concentrations and low DO and summer temperature variability. Conversely, sites with very good and excellent rankings were associated with high DO concentrations and high summer temperature variability and low TP concentrations. Sites with good rankings were associated with higher TN concentrations and lower Fe concentrations (Fig. 2). Stepwise multiple linear regression analysis was used to assess multiple stressor effects and their interactions on HBI scores (Table 2a). Input variables from the three environmental categories consisted of Fe, pH, DO, Q mean, Q min, MST, MDMax, and CVST. The most parsimonious model for predicting differences in HBI scores among sites included DO, pH, Q mean, Q min, MST, MDMin (F6,32 = 16.84, p b 0.001, R2adj = 0.71; Table 2a). DO, pH, Q min, and MST were significantly and negatively correlated, and Q mean and MDMin significantly and positively correlated, with HBI score. DO and MDMin had the largest, while Q mean and Q min had the smallest, relative importance in explaining the total response variance (Table 2a). Based on our selection criteria, there were no

Fig. 2. Linear discriminant analysis (LDA) ordination for 6 classes of benthic invertebrate HBI rankings. Dots indicate the centroids of the ranking scores, while arrows indicate significant environmental variables that discriminate among the rankings: TP = total phosphorus concentration (μg/L), DO = dissolved oxygen concentration (mg/L), Fe = iron concentration (μg/L), TN = total nitrogen concentration (mg/L), and CVST = coefficient of variation of summer temperature (%).

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

6

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

Table 2 Regression models relating (a) HBI and (b) IBI scores to water quality, temperature, and hydrology variables. Response metric

Main effects model R2adj (p)

Main + 1st order interactions R2adj

F

p

(a) HBI scores

0.71 (b0.001)

–

–

–

(b) IBI scores

0.45 (b0.001)

0.49

7.80

b0.001

Parameter

Standardized coefficient

p

Relative importance (%)

DO pH Qmean Qmin MST MDMin TSS Fe Ti Qmin MST TSS*Qmin TSS*MST

−0.67 −0.45 0.60 −0.39 −0.77 0.87 0.59 0.38 −0.51 0.46 −0.39 0.73 −0.36

b0.001 b0.001 0.006 0.04 0.002 b0.001 0.010 0.013 b0.001 b0.001 0.003 0.044 0.037

39.6 17.4 5.0 3.2 11.4 23.5 19.4 10.9 25.4 25.2 4.7 4.5 9.9

DO = dissolved oxygen concentration (mg/L); pH; TSS = total suspended solids concentration (mg/L); Fe = iron concentration (μg/L); Ti = titanium concentration (μg/L); Q mean = average annual streamflow (m3/s); Q min = annual 7-day low flow (m3/s); MST = mean summer temperature (°C); MDMin = mean daily minimum temperature (°C).

interaction terms that significantly improved the amount of explained total response variance compared to the main effects-only model; thus, only additive effects were important determinants of the differences in HBI scores among catchments. Negative effects dominated, with higher values of DO, pH, Q min, and MST contributing to lower HBI scores (representing low anthropogenic disturbance and improved stream health), while the positive effects of Q mean and MDMin contributed to higher HBI scores (Table 2a). Variance partitioning analysis (VPA) was used to explore the relationship between three categories of environmental variables (water quality, hydrology, temperature) and benthic invertebrate assemblage composition. Significant input variables were: Fe, TP, and pH (water quality), MST, MDMax, and CVST (temperature), and Q mean and Q min (hydrology). An RDA biplot of the six catchments analyzed using VPA revealed strong primary gradients when constrained to the significant water quality, temperature, and hydrological variables (Fig. 3). RDA axis 1 (explaining 39.4% of the constrained variance) separated catchments based on their variability in water quality, while RDA axis 2 (explaining 19.4% of the constrained variance) further separated catchments based on the variability in their temperature and flow patterns. With negative RDA 1 scores, the Upper Schomberg and Beaver River catchments are characterized by high concentrations of TP and Fe, or higher pH, respectively, and are dominated by families tolerant of organic pollution. In contrast, the Black River, Pefferlaw River (most years), Lovers Creek and Hawkestone Creek catchments have positive RDA 1 scores, are characterized generally by low Fe and TP concentrations and low pH, and are dominated by families more sensitive to anthropogenic disturbance. The Beaver and Black River catchments are further characterized by higher values of temperature (MDMax, MST) and hydrological (Q mean, Q min) variables, compared to Lovers Creek, Hawkestone Creek, and Upper Schomberg River catchments. A marked shift among years was observed for the Pefferlaw River catchment, moving from negative to positive RDA 1 scores and positive to negative RDA 2 scores between 2004 and 2011, associated with decreases in TP and Fe concentrations and pH, and increases in Q min, and shifts in community composition from tolerant (e.g., Gammaridae, Tubificidae, Planorbidae, Pisidiidae) to more pollution-sensitive families (e.g., Isonychiidae, Heptageniidae, Pholopotamidae, Baetidae, Perlidae) (Fig. 3). For the three-category VPA, the total variance of the benthic invertebrate community composition explained by all predictor categories was 32.7% (p = 0.001) (Fig. 4). The largest component of explained variance was by water quality independent of hydrology and temperature (20.1%). The shared variance between hydrology and temperature explained a smaller portion (5.9%) of the total variance in the benthic invertebrate family-level data; however, this was larger than the variance

explained by either hydrology (3.8%) or temperature (4.1%) independently. This result suggests that the co-variation of hydrology and temperature, while minor, can impact the structuring of these communities at the family level. The complex co-variation among water quality, hydrology, and temperature explained a small portion (1.5%) of the variation, and there was no shared variance of water quality and temperature or water quality and hydrology (Fig. 4).

Fig. 3. Overall redundancy analysis (RDA) ordination biplot for benthic invertebrate communities from six Lake Simcoe catchments (BL = Black River; BV = Beaver River; HS = Hawkestone Creek; LV = Lovers Creek; PF = Pefferlaw River; US = Upper Schomberg River) across 7 years (2004–2007 and 2009–2011). Arrows indicate significant environmental variables that best predict differences among the benthic invertebrate communities: TP = total phosphorus concentration (μg/L); Fe = iron concentration (μg/L); pH = power of hydrogen; CVST = coefficient of variation of summer temperature (%); MDMax = mean daily maximum summer temperature (°C); MST = mean summer temperature (°C); Q mean = average annual streamflow (m3/s); Q min = average annual 7-day low flow (m3/s). Benthic invertebrate families are represented by numbers: 1, Asellidae; 2, Baetidae; 3, Caenidae; 4, Capniidae; 5, Chironomidae; 6, Coenagrionidae; 7, Corydalidae; 8, Elmidae; 9, Empididae; 10, Ephemerellidae; 11, Gammaridae; 12, Helicopsychidae; 13, Heptageniidae; 14, Hyalellidae; 15, Hydropsychidae; 16, Isonychiidae; 17, Naididae; 18, Oligochaeta; 19, Perlidae; 20, Philopotamidae; 21, Physidae; 22, Pisidiidae; 23, Planorbidae; 24, Psephenidae; 25, Taeniopterygidae; 26, Tipulidae; 27, Tubificidae.

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

7

Fig. 5. Linear discriminant analysis (LDA) ordination for three classes of fish IBI rankings. Dots indicate the centroids of the ranking scores, while arrows indicate significant environmental variables that discriminate among the rankings: TP = total phosphorus concentration (μg/L); DO = dissolved oxygen concentration (mg/L); Q max = annual maximum streamflow (m3/s); Qmin = average annual 7-day low flow (m3/s).

Fig. 4. Three-category variance partitioning analysis (VPA) displaying the amount of variation in the benthic invertebrate (BI) and fish (F) communities explained by the unique and shared effects of water quality (Wq), temperature (T), and hydrology (H) variables in six Lake Simcoe catchments from 2004 to 2007 and 2009 to 2011 (BI) or 2004 to 2012 (F). U = unexplained variance.

Fish and multiple stressors LDA was used to discriminate among sites with scores distributed across three quantitative IBI rankings (Fig. 5). Only two LDA axes were necessary to capture all of the variation (LDA 1 = 87.3% of variance; LDA 2 = 12.7% of variance) in IBI rankings among catchments (MANOVA Wilk's test: W2,51 = 0.394, approx F22,82 = 2.21, p = 0.005) (Fig. 5). LDA 1 separated sites with fair from good rankings, while LDA 2 separated the poor rankings (n = 3) from the fair and good rankings. Four significant environmental predictors discriminated among the sites: TP (F2,51 = 5.14, p = 0.009), DO (F2,51 = 3.98, p = 0.025), Q max (F2,51 = 3.19, p = 0.049), and Q min (F2,51 = 4.21, p = 0.020). Sites with fair rankings were associated with higher DO concentrations, while good rankings were associated with higher values of Qmax and Qmin. Poor rankings were associated with low TP concentrations and possibly lower DO concentrations (Fig. 5). Stepwise multiple linear regression analysis was used to assess multiple stressor effects and their interactions on IBI scores (Table 2b). Input variables for multiple linear regression from the three environmental categories consisted of TSS, Fe, Ti, Q min, Q mean, MST, and MDMax. The most parsimonious main effects model for predicting differences in IBI scores among sites included TSS, Fe, Ti, Q min, and MST (F7,41 = 7.80, p b 0.001, R2adj = 0.45; Table 2b). Two interaction terms, TSS*Q min and TSS*MST, slightly improved the amount of explained variance (R2adj = 0.49), although the difference was marginally significant (F2,41 = 3.17, p = 0.053). Q min and Ti had the largest, while MST and the TSS*Q min interaction had the smallest, relative importance in explaining the total response variance. The additive effects of TSS, Fe, and Q min were positive, with higher values contributing to higher IBI scores, while the additive effects of Ti and MST were negative, with

higher values contributing to lower IBI scores. Although the interaction terms had low relative importance (e.g., b 10%), their inclusion in the model demonstrates the impact of multiplicative effects on IBI scores. The antagonistic effect of the TSS*MST interaction demonstrates how the positive effect of TSS can be reduced as MST increases, while the synergistic effect of the TSS*Q min interaction amplifies the independent effects of TSS and Q min. Variance partitioning analysis (VPA) was used to explore the relationship between three categories of environmental variables (water quality, hydrology, temperature) and fish assemblage composition. Significant input variables were: TN, TP, Cl (water quality), MST, MDMin, MDMax (temperature), and Q mean, Q min (hydrology). An RDA biplot of the six catchments analyzed using VPA revealed strong primary gradients when constrained to the significant water quality, temperature, and hydrological variables (Fig. 6). RDA axis 1 (explaining 38.5% of the constrained variance) separated catchments based on temperature (MDMax) and hydrology (Q min, Q mean), while RDA 2 (explaining 26.9% of the constrained variance) separated catchments based on water quality, hydrology, and temperature variables. Hawkestone Creek and Upper Schomberg River catchments loaded positively along RDA 1 and were associated with lower values of Q mean, Q min, and MST and higher abundances of Rhinichthys atratulus, Nocomis micropogon, Etheostoma nigrum, and Semotilus atromaculatus. Lovers Creek catchment and some years in Beaver River catchment loaded positively along RDA 2, with communities dominated by Rhinichthys cataractae and Cottus bairdi, and were associated with high concentrations of TN and Cl and low concentrations of TP. The Black and Pefferlaw River catchments loaded negatively along RDA 2 and were associated with higher concentrations of TP, higher values of temperature (MDMax, MST, MDMin) and hydrology (Q mean, Q min) metrics, and higher abundances of Perca flavescens, Neogobius melanostomus, Micropterus salmoides, Lepomis gibbosus, and Lepomis macrochirus. For the three-category VPA, the total variance in the fish community composition explained by all predictor categories was 47.3% (p = 0.001) (Fig. 4). Water quality, independent of hydrology and temperature, explained the largest amount of the total variance (20.8%),

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

8

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

(LSEMS, 2008). In the LSW, high tributary TP concentrations have been linked to agricultural activities and urban stormwater runoff (LSRCA, 2013a; O'Connor et al., 2013), although they may also be indicative of low levels of natural heritage cover resulting in poor runoff filtering capacity. High concentrations of TP were significant predictors of fair, fairly poor, and poor stream health (HBI rankings) and were associated with benthic invertebrate families tolerant of organic pollution (e.g., Asellidae, Helicopsychidae, Capniidae, Tubificidae), particularly in the Upper Schomberg River catchment (see Fig. 7a). In the Pefferlaw and Black River catchments, high TP concentrations were associated with fish species that have generalist feeding strategies and habitat preferences and warm water tolerances (e.g., Percina caprodes, M. salmoides, L. gibbosus, L. macrochirus) (Coker et al., 2001). In wadeable streams in Wisconsin, Wang et al. (2007) not only found that HBI and IBI measures became more degraded as concentrations of TP

Fig. 6. Overall redundancy analysis (RDA) ordination biplot for fish communities from six Lake Simcoe catchments (BL = Black River; BV = Beaver River; HS = Hawkestone Creek; LV = Lovers Creek; PF = Pefferlaw River; US = Upper Schomberg River) across 9 years (2004–2012). Arrows indicate significant environmental variables that best predict differences among the fish communities: TP = total phosphorus concentration (μg/L); Cl = chloride concentration (mg/L); TN = total nitrogen concentration (mg/L); MDMax = mean daily maximum summer temperature (°C); MST = mean summer temperature (°C); MDMin = mean daily minimum summer temperature (°C); Q mean = average annual streamflow (m3/s); Q min = average annual 7-day low flow (m3/s). Fish species are represented by numbers: 1, Ambloplites rupestris; 2, Ameiurus natalis; 5, Cottus bairdi; 8, Etheostoma caeruleum; 9, Etheostoma nigrum; 10, Lepomis gibbosus; 11, Lepomis macrochirus; 12, Luxilus cornutus; 13, Micropterus dolomieu; 14, Micropterus salmoides; 15, Neogobius melanostomus; 17, Nocomis micropogon; 23, Perca flavescens; 28, Rhinichthys atratulus; 29, Rhinichthys cataractae; 30, Semotilus atromaculatus. The dotted circle indicates the plot region that contains species with overlapping low scores (b0.1) on both axes: 3, Ameiurus nebulosus; 4, Catostomus commersoni; 6, Culaea inconstans; 7, Cyprinus carpio; 16, Nocomis biguttatus; 18, Notropis atherinoides; 19, Notemigonus crysoleucas; 20, Notropis heterodon; 21, Notropis rubellus; 22, Percina caprodes; 24, Percina maculata; 25, Phoxinus eos; 26, Pimephales notatus; 27, Pimephales promelas; 31, Umbra limi.

followed by hydrology (8.8%) and temperature (7.4%). Smaller portions of the fish variance were explained by the co-variation of water quality and temperature (3.4%), hydrology and temperature (7%), and the complex co-variation of water quality, biology, and temperature (0.7%). There was no shared variance of water quality and hydrology independent of temperature (Fig. 4). Discussion Over the period from 2004 to 2012, multiple sets of environmental variables have influenced the ecological health of Lake Simcoe tributaries, including the structure of the benthic invertebrate and fish communities. Water quality was the primary category influencing stream health and biological community composition, followed by temperature and hydrology, respectively. TP, DO, and Fe concentrations, as well as Q min and MST, were consistently identified as important predictors of both fish and benthic invertebrate variance in multivariate analyses. Multiple stressor interactions had both synergistic and antagonistic effects on stream quality as determined by IBI scores. However, overall, stressor interactions explained small portions of the total variance in biological composition. It is not surprising that TP was identified as a significant anthropogenic stressor in our analyses, as excessive phosphorus concentrations have been problematic throughout the watershed since the 1970s

Fig. 7. Total phosphorus concentrations (μg/L) vs. HBI (a) and IBI scores (b). Horizontal dashed line indicates the PWQO of 30 μg/L. Vertical dashed lines separate classes of HBI or IBI rankings. For HBI: E = excellent, VG = very good, G = good, F = fair, FP = fairly poor, P = poor. For IBI: G = good, F = fair, P = poor.

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

increased but also demonstrated a threshold response where these biotic indices changed drastically beyond TP concentrations of 90 μg/L and 70 μg/L, respectively. While almost all TP concentrations measured in our study were lower than these thresholds, similar patterns of stream quality degradation and shifts in aquatic communities with TP enrichment have been documented elsewhere (e.g., Frey et al., 2011; Justus et al., 2010). In contrast, the IBI rankings indicating poor stream quality were associated with low TP concentrations in our study (see Fig. 7b). With fish species, low TP concentrations may indicate low productivity habitats (Dodds, 2007), potentially leading to fewer food resources, and thus lower fish diversity and/or abundance, or alternately, an unmeasured variable correlated with low TP concentrations. However, there were few poor ranked sites (n = 3) in our LDA analysis, suggesting that further exploration of the relationship between variations in TP concentrations and productivity of stream habitats is warranted. DO concentrations were another significant predictor of ecological stream health. High DO concentrations were associated with very good and excellent HBI and fair IBI rankings. DO is a vital component of river ecosystems, influenced by community metabolism (primary production and respiration) (Wang et al., 2003), while influencing water quality through impacts on nutrient cycling and nitrification processes (Kemp and Dodds, 2002). The dynamics of DO are complex and vary depending on rates of photosynthesis and respiration, air–water exchange, temperature, and groundwater inflow (Allan, 1995). Anthropogenic nutrient enrichment (e.g., excessive inputs of N and P) that promotes growth of algal biomass can deplete oxygen concentrations during night-time respiration, resulting in periods of hypoxia or anoxia that can impair biotic integrity (Miltner and Rankin, 1998). In our multivariate ordinations, gradients in DO concentrations generally opposed those of TP concentrations, which is consistent with an indirect effect of nutrient loading on dissolved oxygen. Fe concentrations had different impacts for fish and benthic invertebrates, with a positive effect on IBI scores (indicating healthier stream quality), but were associated with benthic invertebrate communities dominated by families tolerant to organic pollution (e.g., Asellidae, Tubificidae, Gammaridae). In the streams we examined, ~25% of Fe concentrations exceeded the Provincial Water Quality Objectives (PWQO) of 300 μg/L (OMOEE, 1994). While Fe is an important nutrient for algal growth, it may also have indirect toxic effects on aquatic biota. For example, the formation of Fe precipitates on biological surfaces (e.g., gills, eggs, shells) can affect the reproduction and feeding behavior of aquatic animals or may alter the physical characteristics and quality of benthic habitats (Vuori, 1995). Fe may be a more detrimental parameter to benthic invertebrates, given their limited spatial mobility. Fe was occasionally associated with TP in our analyses. Fe is a major component of soils (Holtan et al., 1988), and high Fe and TP concentrations may be an indicator of increased soil erosion, from agricultural sources or organically rich areas such as those dominated by wetlands. Such nutrient-rich streams may provide large amounts of food (e.g., invertebrates, periphyton) to local fish populations (e.g., Lenat and Crawford, 1994), enhancing IBI scores. In the RDA of the fish community, sites from the Lovers Creek catchment, and for a few years from the Beaver River catchment were associated with high concentrations of TN and Cl. High Cl concentrations in streams are typically the result of road salt application in winter months, while nitrogen is often associated with sewage inputs and urban runoff in urban areas and with fertilizer use in agricultural areas (LSRCA, 2013a; Winter et al., 2011). Throughout the LSW, Cl concentrations are generally highest in December to March (Winter et al., 2011). Oni et al. (2015) documented a 16% increase in impervious surface area in the Lovers Creek catchment between 1994 and 2008, leading to ~46% increase in runoff. In 5 of 7 years in Lovers Creek catchment, Cl concentrations approached or exceeded the Canadian Council of Ministers of the Environment (CCME) water quality guideline of 120 mg/L (CCME, 2011). High Cl concentrations are likely to exert negative effects on aquatic life, potentially affecting osmogregulation, reproduction,

9

mortality, and exclusion of sensitive species (CCME, 2011; Winter et al., 2011), while high nitrogen concentrations may contribute to eutrophication (Lewis et al., 2011). In Lovers Creek catchment, fish diversity was generally low and dominated by two species R. cataractae and C. bairdi that are benthic insectivores and have cool- and cold-water thermal preferences (Coker et al., 2001). Lovers Creek catchment has multiple stream reaches classified as cool- and cold-water thermal regimes (LSRCA, 2013a), so the presence of these species may be more related to thermal regime than to water quality. However, Lovers Creek and Beaver River catchments also recorded an overall fair stream quality (mean IBI score for Lovers Creek and Beaver River catchments was 24.6 and 25.1, respectively), suggesting water quality is impacting stream health. Winter et al. (2011) found that Lovers Creek catchment had among the highest Cl concentrations among catchments tested in the LSW, frequently exceeding guideline concentrations (120 mg/L; CCME, 2011). Regional increases in Cl concentrations (Kelly et al. unpubl. data) suggest that urban impacts on stream health will increase as urban areas expand in future. In the VPA analysis, hydrology independently explained little of the total variance in benthic invertebrate and fish assemblages (3.8% and 8.8%, respectively). However, despite the low overall contribution, two hydrology metrics Q min and Q mean were frequently identified as significant predictors of stream health and community change among catchments, particularly for the fish data. For example, in multiple linear regressions, a greater Q min, the 7-day low flow, implied better stream health, although its relative importance was higher for IBI than HBI scores (25.2% vs. 3.2%). Q min also predicted good IBI rankings in the LDA. The productivity and diversity of aquatic communities generally vary with habitat area (Allan, 1995); thus, biological conditions of a stream are expected to be impacted by variations in metrics that influence habitat availability. Variations in Q mean result in changes in stream depth and wetted surface area (Konrad and Booth, 2002), influencing overall habitat volume. Q min provides a measure of the lowest magnitude of stream flow during summer, and higher values may signify an increase in the available habitat for aquatic organisms (Konrad and Booth, 2002). Changes in habitat volumes may thus be more critical for larger-bodied organisms such as fish. Although temperature explained a small overall proportion of benthic invertebrate and fish community variance (4.1% and 7.4%, respectively), several temperature metrics were frequently identified as significant predictors of stream health. Mean summer temperature (MST) had a negative effect on IBI scores, predicting greater stream health with lower values of MST. This pattern likely reflects the thermal tolerances of fish, as IBI scores increase with a greater proportion of coldwater species such as Salvelinus fontinalis (OMNR and TRCA, 2005). For benthic invertebrates, mean daily minimum temperature (MDMin) had a higher relative importance than MST in the multiple regressions, such that increasing MDMin led to higher HBI scores and poor stream health. Increases in the minimum daily temperature could limit the survival of taxa (or developmental stages) sensitive to higher temperatures (e.g., cold water-acclimated species) (Cox and Rutherford, 2000). Further, the coefficient of variation in summer temperature (CVST) was associated with excellent and very good HBI rankings, suggesting variability in stream temperatures could lead to the creation of more niches, supporting higher benthic diversity and thus higher stream quality (Vinson and Hawkins, 1998). Stream temperature is affected by both natural conditions (e.g., precipitation, groundwater inputs, wind speed) and anthropogenic alterations (e.g., removal of bank vegetation, presence of dams, thermal pollution from water treatment plants) (Poole and Berman, 2001). From 2003 to 2011, LSRCA (2013a) detected increases in average daily maximum temperature for sites in the Pefferlaw River catchment. Combined with our results, these trends imply that future changes to stream temperatures due to a changing climate or land uses will have consequences for stream health in the LSW. Based on fish IBI scores, we found evidence that multiple stressor interactions are impacting tributary health in the LSW. Significant

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

10

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

interactions between water quality and temperature (TSS*MST) and water quality and hydrology (TSS*Q min) had differential effects on IBI scores and thus on stream health. Although a pollutant of concern at high levels, the TSS concentrations we observed were below CCME guidelines (≤ 25 mg/L increase over background levels; CCME, 2002) and may be an indicator of more abundant food resources (seston, periphyton), promoting fish abundance and increasing IBI metrics (Lenat and Crawford, 1994). The antagonistic effect of the TSS*MST interaction demonstrates how the apparent positive effect of TSS on IBI scores can be reduced as mean summer temperatures increase. In contrast, the TSS*Q min interaction had a positive (i.e., synergistic) impact on stream health, such that IBI scores would be expected to demonstrate a greater relative increase when both higher TSS concentrations and 7-day low flow values increased in combination than compared to any increases observed in the individual components alone. Further, for both benthic invertebrate and fish communities, the co-variation of hydrology and temperature explained similar or greater proportions of the total variance than their independent categories, suggesting that metrics in these two explanatory categories are changing in concert over the period of record, and/or that temperature and hydrology explain similar aspects of biotic community composition. This was also illustrated in the RDA ordinations, as temperature metrics (MST, MDmax) occurred together with hydrology metrics (Q mean and Q min), characterizing the Beaver and Black River catchments. Together, these results illustrate the complex nature of stressor interactions and the difficulty in predicting how the LSW responds to stress. Overall, multiple stressor interactions explained a small amount of the variation in stream health and community composition, and this pattern was consistent across the different multivariate analyses used. There were no significant interactions predicting HBI scores, while the interaction terms included in the model predicting IBI scores had low relative importance (e.g., both b10%). Variation in stressor effects among trophic levels has been documented elsewhere; Piggott et al. (2012) observed that stream algal communities responded to stressors acting individually, while pair-wise stressor interactions were detected more frequently for benthic invertebrate communities. In our study, the lack of interactions or their low relative importance among stressors was surprising, given that many components of river ecosystems are interrelated and may often have cascading consequences when a single metric is altered. For example, poor stream quality would be expected during periods of elevated water temperatures and lower base flows leading to low DO concentrations. It is possible that such scenarios did not occur frequently within the narrow spatio-temporal range of our dataset, leading to few interactions among stressors. Alternately, the stressor that invokes the largest stress response may overwhelm any smaller responses to another stressor, leading to the detection of few interactions (e.g., Folt et al., 1999). The LSW is an important economic and ecological asset, providing many ecosystem services to the surrounding human population. While significant progress has been made in recent years to improve the ecological health of the LSW (MOECC, 2015), future challenges concerning continued land-use change and climate change are projected to impact many aspects of the ecosystem (Crossman et al., 2013; Oni et al., 2014). Urban development is expected to continue, fueled by a projected population growth of 30% by 2031, while agricultural activates may decline slightly (MOI, 2006). Thus, changes to nutrient and contaminant inputs could continue to impact stream health directly through alterations to TP and Fe inputs, indirectly via changes to DO concentrations, or through alterations of flow regimes. Climate change will probably modify the thermal regime of Lake Simcoe tributaries (Crossman et al., 2013), potentially causing loss of cold-water habitats for sensitive fish and benthic invertebrate species. Multiple stressor interactions may also become more important as watershed components change in the future due to increased nutrient loading, inputs of contaminants, or changes to hydrological or thermal regimes, as a result of future climate change or urban development throughout the

watershed. While the challenge of separating the independent effects of altered flow and thermal regimes makes future projections difficult, our results suggest that maintaining ecologically relevant flow levels (e.g., maximizing Q min), particularly during summer months, could mitigate the expected increases in stream temperature with changing climate. Acknowledgements We wish to thank the invaluable assistance of the Lake Simcoe Region Conservation Authority, particularly Lance Aspden, for providing hydrology data. Katie Stammler (University of Waterloo) provided valuable insights on fish ecology and the treatment of non-detects in chemistry data. Krista Chomicki (TRCA) assisted with calculations of hydrology time series, and Hamdi Jarjanazi (MOECC) provided a map of the watershed. We thank Michelle Palmer (MOECC) and members of the Lake Simcoe Science Committee for their insight and helpful discussions. Finally, we thank two anonymous reviewers whose comments improved the manuscript. This research was supported by an Ontario Ministry of the Environment and Climate Change grant to L.M. References [CCME] Canadian Council of Ministers of the Environment, 2002. Canadian Water Quality Guidelines for the Protection of Aquatic Life: Total Particulate Matter. Canadian Environmental Quality Guidelines, 1999. Canadian Council of Ministers of the Environment, Winnipeg. [CCME] Canadian Council of Ministers of the Environment, 2011. Canadian Water Quality Guidelines for the Protection of Aquatic Life: Chloride. Canadian Environmental Quality Guidelines, 1999. Canadian Council of Ministers of the Environment, Winnipeg. [LSEMS] Lake Simcoe Environmental Management Strategy, 2008. Lake Simcoe basinwide report. Available from: http://www.lsrca.on.ca/pdf/reports/lsems/basin_wide_ report.pdf. [LSRCA] Lake Simcoe Region Conservation Authority, 2005. Lake Simcoe Watershed 2005 Environmental Monitoring Report. Available from: http://www.lsrca.on.ca/reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2010a. Black River Subwatershed Plan. Available from: http://www.lsrca.on.ca/reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2010b. West Holland River Subwatershed Management Plan. Available from: http://www.lsrca.on.ca/reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2012a. Barrie Creeks, Lovers Creek, and Hewitt's Creek Subwatershed Plan. Available from: http://www.lsrca.on.ca/ reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2012b. Pefferlaw River Subwatershed Plan. Available from: http://www.lsrca.on.ca/reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2012c. Beaver River Subwatershed Plan. Available from: http://www.lsrca.on.ca/reports/. [LSRCA] Lake Simcoe Region Conservation Authority, 2013a. Lake Simcoe Watershed 2013 Environmental Monitoring Report (2007–2011 data). Available from: http:// www.lsrca.on.ca/pdf/reports/monitoring_report_2013.pdf. [LSRCA] Lake Simcoe Region Conservation Authority, 2013b. The Oro and Hawkestone Creeks Subwatershed Plan. Available from: http://www.lsrca.on.ca/reports/. [MOECC] Ministry of the Environment and Climate Change, 2015. Lake Simcoe Monitoring Report 2014. Queen's Printer for Ontario (PIBS 9892E). [MOI] Ministry of Infrastructure, 2006. Growth Plan for the Greater Golden Horseshoe. Queen's Printer for Ontario. [OMOE] Ontario Ministry of Environment, 1983. Handbook of Analytical Methods for Environmental Samples. Queen's Printer for Ontario. [OMOEE] Ontario Ministry of Environment and Energy, 1994. Water Management: Policies, Guidelines, Provincial Water Quality Objectives of the Ministry of Environment and Energy. Queen's Printer for Ontario. Allan, J.D., 1995. Stream Ecology: Structure and Function of Running Waters. Kluwer, Dordrecht, Neth. Allan, J.D., 2004. Landscapes and Riverscapes: The Influence of Land Use on Stream Ecosystems. Annu. Rev. Ecol. Evol. Syst. 35, 257–284. http://dx.doi.org/10.1146/ annurev.ecolsys.35.120202.110122. Baulch, H.M., Futter, M.N., Jin, L., Whitehead, P.G., Woods, D.T., Dillon, P.J., Butterfield, D.A., Oni, S.K., Aspden, L.P., O'Connor, E.M., Crossman, J., 2013. Phosphorus dynamics across intensively monitored subcatchments in the Beaver River. Inl. Waters 3, 187–206. http://dx.doi.org/10.5268/iw-3.2.530. Bazinet, N.L., Gilbert, B.M., Wallace, A.M., 2010. A comparison of urbanization effects on stream benthic macroinvertebrates and water chemistry in an urban and an urbanizing basin in Southern Ontario, Canada. Water Qual. Res. J. Can. 45, 327–341. Blanchet, G., Legendre, P., Borcard, D., 2008. Forward selection of spatial explanatory variables. Ecology 89, 2623–2632. http://dx.doi.org/10.1890/07-0986.1. Borcard, D., Gillet, F., Legendre, P., 2011. Numerical Ecology with R, Use R! Springer http:// dx.doi.org/10.1007/9781-4419-7976-6. Borcard, D., Legendre, P., Drapeau, P., 1992. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055.

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx Bourque, C.P.A., Pomeroy, J.H., 2001. Effects of forest harvesting on summer stream temperatures in New Brunswick, Canada: an inter-catchment, multiple-year comparison. Hydrol. Earth Syst. Sci. 5, 599–614. http://dx.doi.org/10.5194/hess-5-599-2001. Bunn, S.E., Arthington, A.H., 2002. Basic principles and ecological consequences of altered flow regimes for aquatic biodiversity. Environ. Manag. 30, 492–507. http://dx.doi.org/ 10.1007/s00267-002-2737-0. Carpenter, S., Caraco, N., Correll, D.L., Howarth, R.W., Sharpley, A.N., Smith, V.H., 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecol. Appl. 8, 559–568. http://dx.doi.org/10.1890/1051-0761(1998)008[0559:NPOSWW]2.0.CO;2. Christensen, M.R., Graham, M.D., Vinebrooke, R.D., Findlay, D.L., Paterson, M.J., Turner, M.A., 2006. Multiple anthropogenic stressors cause ecological surprises in boreal lakes. Glob. Chang. Biol. 12, 2316–2322. http://dx.doi.org/10.1111/j.1365-2486. 2006.01257.x. Coker, G.A., Portt, C.B., Minns, C.K., 2001. Morphological and ecological characteristics of Canadian freshwater fishes. Can. MS Rpt. Fish. Aquat. Sci. 2554 (iv + 89p). Cox, T.J., Rutherford, J.C., 2000. Thermal tolerances of two stream invertebrates exposed to diumally varying temperature. N. Z. J. Mar. Freshw. Res. 34, 203–208. http://dx.doi. org/10.1080/00288330.2000.9516926. Crossman, J., Futter, M.N., Oni, S.K., Whitehead, P.G., Jin, L., Butterfield, D., Baulch, H.M., Dillon, P.J., 2013. Impacts of climate change on hydrology and water quality: future proofing management strategies in the Lake Simcoe Watershed, Canada. J. Great Lakes Res. 39, 19–32. http://dx.doi.org/10.1016/j.jglr.2012.11.003. Dodds, W.K., 2007. Trophic state, eutrophication and nutrient criteria in streams. Trends Ecol. Evol. 22, 669–676. http://dx.doi.org/10.1016/j.tree.2007.07.010. Dray, S., Legendre, P., Blanchet, G., 2013. packfor: Forward Selection with Permutation (Canoco p.46). R package version 0.0–8/r109. http://R-Forge.R-project.org/projects/ sedar/. Evans, D.O., Nicholls, K.H., Allen, Y.C., McMurtry, M.J., 1996. Historical land use, phosphorus loading, and loss of fish habitat in Lake Simcoe, Canada. Can. J. Fish. Aquat. Sci. 53, 194–218. http://dx.doi.org/10.1139/f96-012. Folt, C.L., Chen, C.Y., Moore, M.V., Burnaford, J., 1999. Synergism and antagonism among multiple stressors. Limnol. Oceanogr. 44, 864–877. Frey, J.W., Bell, A.H., Hambrook Berkman, J.A., Lorenz, D.L., 2011. Assessment of Nutrient Enrichment by Use of Algal-, Invertebrate-, and Fish-Community Attributes in Wadeable Streams in Ecoregions surrounding the Great Lakes. U.S. Geological Survey Scientific Investigations Report 2011-5009. Gregory, S.V., Swanson, F.J., McKee, W.A., Cummins, K.W., 1991. An ecosystem perspective of riparian zones. Bioscience 41, 540–551. http://dx.doi.org/10.2307/1311607. Grömping, U., 2006. Relative importance for linear regression in R: the package relaimpo. J. Stat. Softw. 17, 1–27. Helsel, D.R., 2012. Statistics for Censored Environmental Data Using Minitab and R. second ed. Wiley. Hilsenhoff, W.L., 1988. Rapid field assessment of organic pollution with a family-level biotic index. J. N. Am. Benthol. Soc. 7, 65–68. http://dx.doi.org/10.2307/1467832. Holtan, H., Kamp-Nielsen, L., Stuanes, A.O., 1988. Phosphorus in soil, water and sediment: an overview. Hydrobiologia 170, 19–34. http://dx.doi.org/10.1007/BF00024896. Jones, C., Somers, K.M., Craig, B., Reynoldson, T.B., 2004. Ontario Benthos Biomonitoring Network Protocol Manual, Version 1.0. Joy, M.K., Death, R.G., 2004. Application of the index of biotic integrity methodology to New Zealand freshwater fish communities. Environ. Manag. 34, 415–428. http://dx. doi.org/10.1007/s00267-004-0083-0. Justus, B.G., Petersen, J.C., Femmer, S.R., Davis, J.V., Wallace, J.E., 2010. A comparison of algal, macroinvertebrate, and fish assemblage indices for assessing low-level nutrient enrichment in wadeable Ozark streams. Ecol. Indic. 10, 627–638. Karr, J.R., 1991. Biological integrity: a long-neglected aspect of water resource management. Ecol. Appl. 1, 66–84. Kelly, N.E., Yan, N.D., Walseng, B., Hessen, D.O., 2013. Differential short- and long-term effects of an invertebrate predator on zooplankton communities in invaded and native lakes. Divers. Distrib. 19, 396–410. http://dx.doi.org/10.1111/j.1472-4642.2012.00946.x. Kemp, M.J., Dodds, W.K., 2002. The influence of ammonium, nitrate, and dissolved oxygen concentrations on uptake, nitrification, and denitrification rates associated with prairie stream substrata. Limnol. Oceanogr. 47, 1380–1393. http://dx.doi.org/10.4319/lo. 2002.47.5.1380. Konrad, C.P., Booth, D.B., 2002. Hydrologic Trends Associated with Urban Development for Selected Streams in the Puget Sound Basin, Western Washington. U.S. Geological Survey Water-Resources Investigations Report 02-4040. U.S. Geological Survey, Tacoma, Washington. Konrad, C.P., Booth, D.B., Burges, S.J., 2005. Effects of urban development in the Puget Lowland, Washington, on interannual streamflow patterns: consequences for channel form and streambed disturbance. Water Resour. Res. 41, 1–15. http://dx.doi. org/10.1029/2005WR004097. Lee, L., 2012. NADA: Nondetects and Data Analysis for Environmental Data. R package version 1.5–4 (http://cran.r-project.org/package=NADA). Legendre, P., Gallagher, E.D., 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129, 271–280. http://dx.doi.org/10.1007/ s004420100716. Lenat, D.R., Crawford, J.K., 1994. Effects of land use on water quality and aquatic biota of three North Carolina Piedmont streams. Hydrobiologia 294, 185–199. http://dx.doi. org/10.1007/BF00021291. Lewis, W.M., Wurtsbaugh, W.A., Paerl, H.W., 2011. Rationale for control of anthropogenic nitrogen and phosphorus in inland waters. Environ. Sci. Technol. 45, 10030–10035. Lumley, T., 2009. leaps: regression subset selection R package version 2.9, using Fortran code by Alan Miller. http://CRAN.R-project.org/package=leaps. Merovich, G.T., Petty, J.T., 2007. Interactive effects of multiple stressors and restoration priorities in a mined Appalachian watershed. Hydrobiologia 575, 13–31. http://dx. doi.org/10.1007/s10750-006-0277-y.

11

Miles, J.J., Eimers, M.C., North, R.L., Dillon, P.J., 2013. Spatial distribution and temporal variability in the forms of phosphorus in the Beaver River subwatershed of Lake Simcoe, Ontario, Canada. Inl. Waters 3, 179–186. http://dx.doi.org/10.5268/IW-3.2.531. Miltner, R.J., Rankin, E.T., 1998. Primary nutrients and the biotic integrity of rivers and streams. Freshw. Biol. 40, 145–158. http://dx.doi.org/10.1046/j.1365-2427.1998. 00324.x. Moore, A., Palmer, M., 2005. Invertebrate biodiversity in agricultural and urban headwater streams: implications for conservation and management. Ecol. Appl. 15, 1169–1177. North, R.L., Winter, J.G., Dillon, P.J., 2013. Nutrient indicators of agricultural impacts in the tributaries of a large lake. Inl. Waters 3, 221–234. http://dx.doi.org/10.5268/IW-3. 2.565. O'Connor, E.M.O., Aspden, L., Lembcke, D., Young, J., Stainsby, E.A., Winter, J.G., 2013. Annual Water Balances and Total Phosphorus Loads to Lake Simcoe (2007–2009). Lake Simcoe Region Conservation Authority (Available from: http://www.lsrca.on.ca/pdf/ reports/phosphorus_etal_2007_2009.pdf). O'Connor, E.M., McConnell, C., Lembcke, D., Winter, J.G., 2011. Estimation of total phosphorus loads for a large, flashy river of a highly developed watershed-seasonal and hysteresis effects. J. Great Lakes Res. 37, 26–35. http://dx.doi.org/10.1016/j.jglr. 2011.04.004. O'Connor, E.M., McConnell, C., Winter, J.G., Lembcke, D., Stainsby, E., 2009. Annual Water Balances, Total Phosphorus Budgets and Total Nitrogen and Chloride Loads to Lake Simcoe (2004–2007). Lake Simcoe Region Conservation Authority (Available from: http://www.lsrca.on.ca/pdf/reports/phosphorus_etal_2004_2007.pdf). Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R., O'Hara, R.B., Simpson, G.L., Solymos, P., Stevens, M.H.H., Wagner, H., 2015. vegan: Community Ecology Package. R package version 2.2–1 (http://CRAN.R-project.org/package=vegan). OMNR, TRCA, 2005. Humber River Fisheries Management Plan. Published by the Ontario Ministry of Natural Resources and the Toronto and Region Conservation Authority. Queens Printer for Ontario. Oni, S.K., Futter, M.N., Dillon, P.J., 2011. Landscape-scale control of carbon budget of Lake Simcoe: a process-based modelling approach. J. Great Lakes Res. 37, 160–165. http:// dx.doi.org/10.1016/j.jglr.2010.05.003. Oni, S.K., Futter, M.N., Dillon, P.J., 2015. Hydrological consequences of urban developments in Lake Simcoe watershed, Canada: a combined paired-watershed and change detection approach. Hydrol. Process. 29, 1829–1843. http://dx.doi.org/10.1002/hyp.10290. Oni, S.K., Futter, M.N., Molot, L.A., Dillon, P.J., 2014. Adjacent catchments with similar patterns of land use and climate have markedly different dissolved organic carbon concentration and runoff dynamics. Hydrol. Process. 28, 1436–1449. http://dx.doi.org/10. 1002/hyp. Ormerod, S.J., Dobson, M., Hildrew, A.G., Townsend, C.R., 2010. Multiple stressors in freshwater ecosystems. Freshw. Biol. 55, 1–4. http://dx.doi.org/10.1111/j.1365-2427.2009. 02395.x. Paterson, A.M., Winter, J.G., Nicholls, K.H., Clark, B.J., Ramcharan, C.W., Yan, N.D., Somers, K.M., 2008. Long-term changes in phytoplankton composition in seven Canadian Shield lakes in response to multiple anthropogenic stressors. Can. J. Fish. Aquat. Sci. 65, 846–861. http://dx.doi.org/10.1139/F08-022. Piggott, J.J., Lange, K., Townsend, C.R., Matthaei, C.D., 2012. Multiple stressors in agricultural streams: a mesocosm study of interactions among raised water temperature, sediment addition and nutrient enrichment. PLoS One 7, e49873. http://dx.doi.org/ 10.1371/journal.pone.0049873. Poole, G.C., Berman, C.H., 2001. An ecological perspective on in-stream temperature: natural heat dynamics and mechanisms of human-caused thermal degradation. Environ. Manag. 27, 787–802. http://dx.doi.org/10.1007/s002670010188. R Core Team, 2015. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (https://www.R-project.org/). Scott, L.D., Winter, J.G., Girard, R.E., 2006. Annual water balances, total phosphorus budgets and total nitrogen and chloride loads for Lake Simcoe (1998–2004). Lake Simcoe Environmental Management Strategy Implementation Phase III Technical Report No. Imp. A. 6 (Available from: http://www.lsrca.on.ca/pdf/reports/lsems/water_balance_ phos_load.pdf). Stanfield, L., 2010. Ontario Stream Assessment Protocol. Version 8.0. Fisheries Policy Section. Ontario Ministry of Natural Resources, Peterborough, Ontario (376 pp.). Stanfield, L., Kilgour, B., 2006. Effects of percent impervious cover on fish and benthos assemblages and instream habitats in Lake Ontario tributaries. Am. Fish. Soc. Symp. 577–599. Stanfield, L.W., Kilgour, B.W., 2012. How proximity of land use affects stream fish and habitat. River Res. Appl. 7, 189. http://dx.doi.org/10.1002/rra. Thrush, S.F., Hewitt, J.E., Hickey, C.W., Kelly, S., 2008. Multiple stressor effects identified from species abundance distributions: interactions between urban contaminants and species habitat relationships. J. Exp. Mar. Biol. Ecol. 366, 160–168. http://dx.doi. org/10.1016/j.jembe.2008.07.020. Townsend, C.R., Uhlmann, S.S., Matthaei, C.D., 2008. Individual and combined responses of stream ecosystems to multiple stressors. J. Appl. Ecol. 45, 1810–1819. http://dx. doi.org/10.1111/j.1365-2664.2007.0. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics with S. fourth ed. Springer, New York, USA. Vinson, M.R., Hawkins, C.P., 1998. Biodiversity of stream insects: variation at local, basin, and regional scales. Annu. Rev. Entomol. 43, 271–293. Vuori, K.M., 1995. Direct and indirect effects of iron on river ecosystems. Ann. Zool. Fenn. 32, 317–329. Wang, H., Hondzo, M., Xu, C., Poole, V., Spacie, A., 2003. Dissolved oxygen dynamics of streams draining an urbanized and an agricultural catchment. Ecol. Model. 160, 145–161. http://dx.doi.org/10.1016/S0304-3800(02)00324-1. Wang, L., Robertson, D.M., Garrison, P.J., 2007. Linkages between nutrients and assemblages of macroinvertebrates and fish in wadeable streams: implication to nutrient criteria development. Environ. Manag. 39, 194–212. http://dx.doi.org/10.1007/ s00267-006-0135-8.

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007

12

N.E. Kelly et al. / Journal of Great Lakes Research xxx (2016) xxx–xxx

Whitehead, P.G., Jin, L., Baulch, H.M., Butterfield, D.A., Oni, S.K., Dillon, P.J., Futter, M., Wade, A.J., North, R., O'Connor, E.M., Jarvie, H.P., 2011. Modelling phosphorus dynamics in multi-branch river systems: a study of the Black River, Lake Simcoe, Ontario, Canada. Sci. Total Environ. 412, 315–323. http://dx.doi.org/10.1016/j.scitotenv.2011.09.073. Wilson, S., 2008. Lake Simcoe Basin's Natural Capital: The Value of the Watershed's Ecosystem Services. Friends of the Greenbelt Foundation Occasional Paper Series, p. 45.

Winter, J.G., Landre, A., Lembcke, D., O'Connor, E.M., Young, J.D., 2011. Increasing chloride concentrations in Lake Simcoe and its tributaries. Water Qual. Res. J. Can. 46, 157–165. http://dx.doi.org/10.2166/wqrjc.2011.124. Young, J.D., La Rose, J.K.L., 2014. Lake Simcoe Comprehensive Monitoring Strategy. Ontario Ministry of Environment and Climate Change and Ministry of Natural Resources and Forestry (Available from: http://www.lsrca.on.ca/pdf/reports/moecc-cms.pdf).

Please cite this article as: Kelly, N.E., et al., Multiple stressor effects on stream health in the Lake Simcoe Watershed, J. Great Lakes Res. (2016), http://dx.doi.org/10.1016/j.jglr.2016.07.007