Relating environmental variables with aquatic community structure in an agricultural/urban coldwater stream

Abstract

Background

Urban areas are often built along large rivers and surrounded by agricultural land. This may lead to small tributary streams that have agricultural headwaters and urbanized lower reaches. Our study objectives assessed are as follows: (1) landscape, geomorphic, and water quality variables that best explained variation in aquatic communities and their integrity in a stream system following this agricultural-to-urban land use gradient; (2) ways this land use gradient caused aquatic communities to differ from what would be expected for an idealized natural stream or other longitudinal gradients; and (3) whether the impacts of this land use gradient on aquatic communities would grow larger in a downstream direction through the agricultural and urban developments. Our study area was an impaired coldwater stream in Michigan, USA.

Results

Many factors structured the biological communities along the agricultural-to-urban land use gradient. Instream woody debris had the strongest relationship with EPT (Ephemeroptera, Plecoptera, and Trichoptera) abundance and richness and were most common in the lower, urbanized watershed. Fine streambed substrate had the strongest relationship with Diptera taxa and surface air breather macroinvertebrates and was dominant in agricultural headwaters. Fish community assemblage was influenced largely by stream flow and temperature regimes, while poor fish community integrity in lower urban reaches could be impacted by geomorphology and episodic urban pollution events. Scraping macroinvertebrates were most abundant in deforested, first-order agricultural headwaters, while EPT macroinvertebrate richness was the highest downstream of agricultural areas within the urban zone that had extensive forest buffers.

Conclusion

Environmental variables and aquatic communities would often not conform with what we would expect from an idealized natural stream. EPT richness improved downstream of agricultural areas. This shows promise for the recovery of aquatic systems using well-planned management in watersheds with this agricultural-to-urban land use pattern. Small patches of forest can be the key to conserving aquatic biodiversity in urbanized landscapes. These findings are valuable to an international audience of researchers and water resource managers who study stream systems following this common agricultural-to-urban land use gradient, the ecological communities of which may not conform with what is generally known about land use impacts to streams.

Introduction

Many large cities are located along navigable waterbodies because of historic and/or present day uses of the rivers such as commercial transport, hydroelectric production, and waste disposal (Gallup et al. 1999). Oftentimes, the soils in large river valleys are very fertile, which also encourages agricultural land use in river floodplains (Gallup et al. 1999). As a result of this land use pattern of urban areas along waterbodies surrounded by agricultural lands, numerous small, tributary stream watersheds have become characterized by predominantly agricultural land use in the headwaters and urban land use in the lower reaches. These waterbodies are surrounded by dense populations and often provide numerous recreational opportunities, thereby generating interest in their restoration (Moerke and Lamberti 2004; Alexander and Allan 2007; Schwartz and Herricks 2007). However, agricultural and urban land use can have multiple impacts on stream systems, and understanding their interactions is necessary for successful restoration efforts in these systems (Poff 1997; Allan 2004; Cooper et al. 2009).

The distribution and abundance of stream biota often respond to the interaction of multiple environmental variables (Raleigh et al. 1984; Poff 1997), which can be caused by gradients of land use (Baumgartner and Robinson 2017). We refer to a land use gradient as a change in the relative intensities of different land use types from one point to another, causing a transition from one class of environmental stressors to another. Variation in ecosystem structure can also be caused by longitudinal stream gradients in large, natural streams such as the River Continuum Concept (Vannote et al. 1980; Tornwall et al. 2015). Rivers and streams are known to exhibit longitudinal patterns of many physical properties including average stream width and depth, shading, velocity, and habitat characteristics. Such patterns are central to the River Continuum Concept and its more recent elaborations (Vannote et al. 1980; Tornwall et al. 2015). The River Continuum Concept deals mainly with patterns that occur over broad ranges of stream order (i.e., 1–12) in watersheds without significant anthropogenic impacts (Vannote et al. 1980).

For management purposes, a greater diversity of aquatic macroinvertebrates, especially those sensitive to environmental stressors, can indicate a healthier stream, whereas a greater abundance of pollution-tolerant organisms indicates poor stream quality (Hilsenhoff 1987). The feeding guilds and other traits of stream organisms can also reflect environmental variables (Usseglio-Polatera et al. 2000; Merritt et al. 2006). The composition of aquatic communities can be vastly different depending on the land use surrounding streams, for instance urban and agricultural areas (Nascimento et al. 2018), as well as longitudinally (Vannote et al. 1980). Low-agriculture areas with extensive vegetation buffers may support the greatest biological integrity for aquatic communities because of the habitats there, such as pools and woody debris (Effert-Fanta et al. 2019). Conversely, urban areas often have the lowest richness and diversity of aquatic macroinvertebrates (Nascimento et al. 2018; Fierro et al. 2019).

Relationships between these variables and biological communities can be difficult to characterize quantitatively because they are often complex, with numerous variables impacting multiple facets of biological health (Poff 1997; Johnson et al. 2007; McNair 2009) (Table S1). However, a combination of approaches can be used to analyze and understand these relationships (Baumgartner and Robinson 2017; Fierro et al. 2019). The method of assessing relationships can include water quality and pollution variables, such as electrical conductivity and herbicides (Berger et al. 2017). It can also include land use coverage variables, such as agriculture or urban coverage (Baumgartner and Robinson 2017), habitat variables (Fierro et al. 2019), or a combination of the above (Lemm and Feld 2017). Environmental gradients are useful to examine the spatial variability of natural ecosystems (Whittaker 1967) and environments that are anthropogenically altered (McDonnell and Pickett 1990). Since both agricultural and urban areas often have a highly developed core with more natural areas interspersed and also radiating outward, the steepness of the change from rural to more anthropogenically modified will determine ecosystem structure and function (McDonnell and Pickett 1990).

Our objective was to answer the questions:

• 1) Among landscape, geomorphic, and water quality/pollution variables, what best explains variation in aquatic communities and assessment of their integrity in a stream system with this agricultural-to-urban land use gradient?

• 2) Does this land use gradient cause aquatic communities to not conform to longitudinal gradients such as the idealized natural stream from the River Continuum Concept?

• 3) Do impacts to aquatic communities from this land use gradient grow larger in a downstream direction through the agricultural and urban developments, or are there other factors that can improve biotic integrity despite these land use stressors?

We hypothesized that this agricultural-to-urban land use gradient would cause aquatic communities to differ from what would be expected of an unimpaired natural stream, and that these impacts would increase in a downstream direction through the agricultural and urban landscapes. We designed a study to relate environmental variables with macroinvertebrate and fish communities in the Indian Mill Creek watershed of Michigan, USA. Based on land use, the watershed shifts from being dominated by agriculture along the east and west branches, as well as the upper reaches of the mainstem, to urban development along the lower reaches of the mainstem (Fig. 1). This is similar to other tributaries of major rivers in southern Michigan, the USA, and temperate areas of the world that originate in agricultural areas surrounding urban developments. This study’s findings contribute to our understanding of the response of aquatic communities to anthropogenic disturbance and provide critical information on the effective management of human dominated ecosystems.

Fig. 1

Environmental variable, macroinvertebrate, and fish inventory sites in the Indian Mill Creek (IMC) watershed of Michigan, USA (2017). BC Brandywine Creek

Materials and methods

Study area

Indian Mill Creek (HUC 040500060504) is a third-order tributary to the Grand River in Kent County, Michigan, USA (Fig. 1, Fig. S1). It is 18.5 km long with a 44-km² watershed. Indian Mill Creek is a coldwater stream with stable to moderate temperature fluxes, conducive to rainbow trout (Oncorhynchus mykiss), brown trout (Salmo trutta), and sculpins (Cottus spp.) in the lower reaches (Wehrly et al. 1999; Sigdel 2017). The amount of groundwater influx to the creek increases along an upstream to downstream continuum (Sigdel 2017). Land use in the watershed is predominately urban (43%) and agricultural (39%) (LGROW 2011). Commercial and residential development is mostly in the lower reaches as it flows through the Grand Rapids metropolitan area, natural and urban lands are in the middle reaches, and farmland and orchards are in the agricultural upper reaches (Fig. 1, Fig. S2) (LGROW 2011). Impervious surfaces cover 12% of the watershed and up to 25% of some lower catchments (AWRI, unpublished data; Sigdel 2017). Indian Mill Creek is designated as a coldwater trout stream by the State of Michigan; however, it is currently listed as impaired by the Michigan Department of Environment, Great Lakes, and Energy (EGLE) due to degraded fish communities (Goodwin et al. 2016).

Geologic and geomorphic features of the Indian Mill Creek watershed were formed by retreating glaciers that deposited hills of medium-textured till in the upper watershed, which contributed cobble and rock to the creek (Farrand and Bell 1982). Glacial meltwater carved the larger Grand River Valley (Larson and Schaetzl 2001). Overall, the creek descends 65 m in elevation from its headwaters to its mouth. The steepest descent is along five stream kilometers starting downstream from the present location of Interstate 96 and descending 24 m in elevation (Fig. 2) (Gesch et al. 2002). The lower watershed gently slopes in an outwash of sand and gravel with postglacial alluvium (Farrand and Bell 1982). One low-head dam (Richmond Dam) is present just upstream of site IMC7 (Fig. 1).

Fig. 2

Indian Mill Creek (IMC) of Michigan, USA, and its tributaries following a topographical gradient (elevation data from Gesch et al. 2002). Stars indicate sampling sites. BC Brandywine Creek, WAD Walker Avenue Ditch

Sampling sites

Nine sites in the Indian Mill Creek watershed were monitored for environmental variables and macroinvertebrates (Fig. 1). Seven sites were monitored for fish (Table S2). All sites were chosen in public or private locations with permitted access and to reflect the watershed’s spatial variation (Fig. 1). Environmental variable inventory sites overlapped with all macroinvertebrate and four fish assessment sites. All locations within the watershed were perennially flowing except for the Walker Avenue Ditch site, which becomes dry during long periods between storms.

Environmental variable inventory

The environmental variable inventory assessed attributes of riffle/pool habitat variability, riparian and bank structure, substrate composition, instream woody debris, suspended sediment, stream velocity, bedload sediment, water temperature, stream slope, and episodic urban pollution events (oil, grease, and industrial byproducts) (Table S3). Habitat components were surveyed in June and July 2017. Riffles, pools, and other geomorphic habitats (runs, glides, and cascades) were surveyed using a modified Basinwide Visual Estimation Technique (Dolloff et al. 1993; Tip of the Mitt Watershed Council 2015). Riparian and bank structure conditions were documented using the Great Lakes Environmental Assessment Section (GLEAS) Procedure 51 habitat survey (EGLE 2008). A riparian and bank structure score out of 60 was assigned as the sum of these scores for each site. Scores of 49–60 represented excellent riparian and bank structure, 31–48 were good, 13–30 were fair, and 0–12 were poor. Factors such as habitat variability and streambank stability led to better scores. Streambed substrate was examined using the zigzag method of Wolman pebble counts (Bevenger and King 1995), with 100 to 125 pebble samples per site. Median particle size was calculated for pebble count data (Bevenger and King 1995), as well as the percentage of the count that was fine substrate under 2 mm along the intermediate axis. Woody debris was surveyed following methods of Cordova et al. (2007), who counted all wood pieces greater than 1 m in length and 10 cm in diameter.

Suspended sediment was sampled monthly from May through September 2017, with two additional sampling events immediately after storms using grab samples. These samples were collected in 1-l, polyethylene bottles in the center of the stream at mid-depth and suspended sediment concentration was analyzed by APHA method 2540 D (Greensberg et al. 1992). Stream velocity was measured during these same events using a Marsh-McBirney Flow Mate 2000 velocity meter (Hach Company, Loveland, Colorado, USA) attached to a top-setting wading rod. Seven to twelve velocity measurements were taken across the stream at 60% depth. The maximum velocity of each monitoring event for each site was averaged to represent the water velocity for the site. Bedload sediment was sampled using a Helley-Smith Sampler (Bunte et al. 2008).

Stream temperature was recorded every 30 min in July and August 2017 using automated Tidbit loggers (Onset Computer Corporation, Bourne, Massachusetts, USA). We determined that loggers were in the water throughout the entire deployment after checking data for temperature spikes or fluctuations that would suggest otherwise. Thermal and fluctuation regimes were examined (Wehrly et al. 1999), as well as the number of hours the stream temperature was above the optimal brown trout temperature range of 12 to 19 °C (Raleigh et al. 1984). Channel slope was measured in a Geographic Information System over 30-m elevation data (Gesch et al. 2002).

Episodic urban pollution

Episodic urban pollution events occur irregularly and are typically ephemeral so directly assessing them in real time is difficult. To assess the magnitude of pollution events that have occurred in Indian Mill Creek, we enumerated the number of violation notices for pollution events that have been logged in the EGLE’s MiWaters website (https://miwaters.deq.state.mi.us/nsite/ (accessed 7/16/2020)) in the last decade. Pollution events were recorded if hazardous materials (oil, grease, industrial byproducts) were released either directly into Indian Mill Creek or into a drain that enters the creek. These data were not integrated into the quantitative analyses of environmental variables at each site because ecological impacts could not be established in real time on a site-by-site basis, and toxicity data from the time of the study was limited to three sites. Rather, we enumerated the number of discharges that had occurred and the description of contaminants, then presented it in summary as a potential explanation of the results in the watershed over a longer span of time. The boundary of this assessment is the Indian Mill Creek watershed, within which a pollution discharge would affect reaches downstream, as the compounds would flow downstream, but would not affect reaches upstream of the discharge.

Within the last decade, EGLE has cited 16 facilities either for spill incidents within the Indian Mill Creek watershed that would drain to the creek or for directly discharging contaminated wastewater into Indian Mill Creek. Illicitly discharged and spilled materials have included oil, sodium hydroxide, and metal plating wastewater effluent (EGLE MiWaters Explorer, https://miwaters.deq.state.mi.us/nsite/ (accessed 4/24/2020)). In 1998, ammonia refrigerant from a meat-packing facility spilled into Indian Mill Creek, resulting in a complete fish kill in the ~3 km reach from the discharge point to its confluence with the Grand River (Hanshue 1998). The sites included in this section are U/S dam, IMC7, and Turner Avenue.

Whole sediment toxicity tests of samples from the lower reach of Indian Mill Creek, collected in 2017, resulted in reduced 10-day growth of Chironomus dilutus and reduced 10- and 28-day growth and survival of Hyalella azteca, which may be a result of elevated sediment PAH (polycyclic aromatic hydrocarbon) concentrations (Parker 2018). Parker (2018) found a significant decrease in growth and survival of Hyalella in three lower Indian Mill Creek sediment samples tested for toxicity, suggesting that the sediment may have picked up some of the episodic pollution events. The three sites that Parker (2018) tested for Hyalella toxicity in 2017 were sites “U/S dam,” “IMC7,” and “Turner Avenue” in the present study. Also, Parker (2018) tested three streambed sediment samples for 23 different pesticides in lower Indian Mill Creek in 2017 and none were found. The three sites that Parker (2018) tested sediment at in 2017 were sites “U/S dam,” “IMC7,” and “Turner Avenue” in the present study as well. Additionally, Sigdel (2017) performed a Microtox assay at 18 sites in the Indian Mill Creek watershed and concluded that, although some seeps and drains showed low toxicity levels, routine stormwater toxicity did not likely have a major effect on aquatic communities.

Benthic macroinvertebrate and fish surveys

Macroinvertebrates and fish were surveyed in July 2017 with GLEAS Procedure 51, used by EGLE to evaluate fish and macroinvertebrate communities and habitat of wadeable streams throughout Michigan (EGLE 2008; Riseng et al. 2010). Both macroinvertebrates and fish were assessed using Procedure 51. The framework of Procedure 51 surveys is a regionally modified Index of Biotic Integrity (Karr 1991). Procedure 51 metrics successfully assess differences in stream communities based on physical variables, providing comparable systematic evaluations among sites (Haller 2010; Woodruff et al. 2010) (Table S3). Our objective for using the Procedure 51 assessment was so the results of this study can be integrated with biological assessments of streams throughout Michigan, which are assessed and managed using Procedure 51. The assessments also allow us to characterize macroinvertebrate and fish communities and how they are affected by environmental variables. Procedure 51 metrics can also be normalized with scores from other assessments (Wiley et al. 2003).

For macroinvertebrates, the Procedure 51 assessment relies on fixed-count (300 ± 60 individuals) subsampling, which is widely used to reduce costs and time for assessing impairments (Barbour and Gerritsen 1996). It also relies on multi-habitat sampling, which best represents community structure, and combines it into an overall index of stream integrity (Haller 2010; Einheuser et al. 2012). For the macroinvertebrate assessments, we used fine-mesh dip nets and hand picking. We sampled until all habitats were represented and the fixed count criteria was met, identifying macroinvertebrates to family level per Procedure 51 protocol. Family-level identification has previously been used in published literature examining impacts (Le Viol et al. 2009; Davies et al. 2010; Dallas and Rivers-Moore 2012) and in particular when associated with this technique (Riseng et al. 2010). The lengths and areas of each stream macroinvertebrate inventory varied among sites (Table S4). Additional calculated macroinvertebrate metrics included the Family Biotic Index (FBI) (Hilsenhoff 1988), total taxa richness, Shannon diversity, and Pielou’s evenness (Qu et al. 2017) (Table S3). The benefit of combining Procedure 51 with other indices such as FBI is that it can allow perspectives at different scales and a better representation of ecological condition (Ogren and Huckins 2014). Macroinvertebrate functional feeding groups and traits were also assigned to each taxon (Bouchard et al. 2004; Merritt et al. 2006, 2008) (Table S3, Table S5).

Fish were sampled with a backpack electroshocker in July 2017 at seven fish study sites (Fig. 1, Tables S2 and S3). Indices of Biotic Integrity, such as Procedure 51, have previously been used to assess fish communities in both warmwater (Einheuser et al. 2013) and coldwater (Lyons et al. 1996) streams and can be correlated with land use scenarios and anthropogenic stressors (Einheuser et al. 2013). Per Procedure 51 protocol (EGLE 2008), sites were sampled for fish communities with a single pass in a section of stream that was ten times the width of the stream. Fish were identified to species, enumerated, and released back into the stream. Each site was rated using the EGLE (2008) Procedure 51 scoring scheme. This scheme automatically considers a site to be “poor” and thus not attaining its fishery designated use if fewer than 50 fish are caught or anomalies (e.g., scoliosis or open lesions) are found on greater than 2% of fish at a site. If ≥ 50 fish are collected, the percentage of salmonids relative to total fish number needs to exceed 1% for a stream to meet its coldwater fisheries designated use.

Although the fish and macroinvertebrate surveys all occurred within the same month (and habitat surveys took place within an overlapping 2-month period), they did not all occur on the same day due to the amount of fieldwork involved. Also, the fish survey did not take place at all the macroinvertebrate sites because of limited availability of our electroshocker and dense vegetation overhanging small agricultural ditches that would make electroshocking difficult.

Data analysis

To answer our first question, we analyzed relationships between environmental variables and aquatic communities in the watershed. Spearman’s rank correlation, Canonical Correspondence Analysis (CCA), and direct comparisons of abundances and metrics were used to relate environmental variables with macroinvertebrate variables and understand the impacts of longitudinal and land use gradients. CCA analyses were performed using the Vegan package of R 3.3.2 (R Core Team 2016; Oksanen et al. 2017). This method was chosen because it is widely used in aquatic sciences (e.g., Fierro et al. 2019), which often contain both constrained and unconstrained datasets, and zero-inflated data (ter Braak and Verdonschot 1995). To confirm our interpretations of the CCA were significant, we performed Permutational Multivariate Analysis of Variance (PERMANOVA) tests with Bray-Curtis dissimilarity matrices and 100,000 permutations on each CCA’s constraining variables (Anderson 2017). Data were transformed using a log+1 transformation prior to the use in the CCA to minimize biases. The variable “Poor Function” was used to represent poor Procedure 51 macroinvertebrate scores, on a continuous scale with a stronger “Poor Function” representing biological communities that rated poorer by the Procedure 51 scoring (EGLE 2008). Woody debris and riparian condition were strongly multicollinear in the CCA (r = 0.72, p = 0.028), so riparian condition was removed from the CCA and the variable was labeled woody debris and riparian condition. Also, the bedload sediment, slope, and riffle constrainers were removed from the CCA analyses because of weak associations with the axes. Taxa richness was removed for the same reasons. Rather than relying on the CCA results alone, a Nonmetric Multidimensional Scaling (NMDS) analysis was performed alongside each CCA to assess robustness of the main results (Oksanen et al. 2017). Additionally, scatterplots were visualized to confirm that CCA results were reasonable.

Due to the low number of fish sample sites, NMDS alone was used to produce a biplot displaying fish communities. NMDS was performed on Bray-Curtis dissimilarity matrices calculated from raw species abundances and standardized by maximum abundances for 400 iterations (Faith et al. 1987; McCune et al. 2002). To verify visual interpretations of fish community groupings in the NMDS biplot, a multi-response permutation procedure (MRPP) (Mielke 1984; Zimmerman et al. 1985) was performed. Euclidean distance measures and a natural weighting, recommended by Mielke (1984), were used in the MRPP. Significance was defined as α = 0.10 for the MRPP comparisons because of the low sample size. The influence of environmental characteristics in the NMDS of fish communities was based on interpretation of characteristics of site groups along the NMDS dimensions, which represented relationships among fish communities. Abundances and metrics were also directly analyzed among fish survey locations to better understand the influence of environmental variables and gradients on the data.

To answer our second and third questions, we evaluated whether our findings conformed with what we would expect for longitudinal gradients in the watershed such as for an unimpacted natural stream, then assessed spatially whether impacts to aquatic communities increased in a downstream direction. Longitudinal gradients in environmental variables along the East Branch, West Branch, and mainstem Indian Mill Creek were assessed using Spearman’s rank correlation (r) to determine whether each variable showed a statistically significant correlation with distance downstream of the apparent stream source. This helped disentangle the impacts of the longitudinal gradients from those of land use alone. Similar to the first question, we also related abundances and metrics with environmental conditions for variables throughout the watershed using results from analyses such as the Spearman’s rank correlation and CCA. Our determination of the impacts of land use and longitudinal gradients, and therefore testing of our hypothesis, was based on the findings of these analyses as they related to specific environmental variables, aquatic communities, and locations in the watershed.

Results

We found that environmental variables and aquatic communities in the watershed were affected by longitudinal and land use gradients. In the results tables, environmental variables (Table 1), macroinvertebrate metrics (Table 2), macroinvertebrate traits (Table 3), and fish metric (Table 4) results are ordered upstream to downstream (i.e., IMC1 to IMC7), which corresponds to an agricultural to urban land cover gradient, followed by the two tributary sites. Site specific data for fish and macroinvertebrate taxa are included in the Supplementary Data (Table S6 and S7).

Table 1 Environmental variable results from the habitat surveys, water quality monitoring, and temperature loggers in Indian Mill Creek (IMC), Michigan, USA (2017) (mean ± standard deviation). SSC suspended sediment concentration, BC Brandywine Creek, WAD Walker Avenue Ditch

Table 2 Macroinvertebrate metrics results calculated using data from the Procedure 51 (P 51) surveys in Indian Mill Creek (IMC), Michigan, USA (2017). EPT Ephemeroptera, Plecoptera, and Trichoptera macroinvertebrates; FBI Family Biotic Index; BC Brandywine Creek; WAD Walker Avenue Ditch

Table 3 Macroinvertebrate trait results assessed using data from the Procedure 51 surveys in Indian Mill Creek (IMC), Michigan, USA (2017). Descriptions of these classifications and sources can be found in Table S5. BC Brandywine Creek, WAD Walker Avenue Ditch

Table 4 Fish metric results using data from the Procedure 51 (P 51) surveys in Indian Mill Creek (IMC), Michigan, USA (2017). U/S upstream of the dam

Environmental variables

Longitudinal and land use patterns were realized for many environmental variables in the watershed. Results for our longitudinal gradient analysis showed that the prevalence of riffles increased toward the mouth of Indian Mill Creek (r = 0.83, p = 0.021; Fig. S3a), although no clear pattern was observed for pools (r = 0.39, p = 0.400; Fig. S3b). Agricultural areas in the upper watershed had the least habitat variability, with virtually no riffles and limited pool habitat (Table 1). A similar pattern occurred for riparian score, with healthier riparian zones in lower reaches (r = 0.76, p = 0.049; Fig. S3c). Riparian conditions can be explained by the land use gradient throughout the watershed. The poorest scoring riparian conditions were in agricultural areas of the upper watershed (Table 1), where scores as low as 12 (poor, IMC2) were documented. Riparian conditions were fair to good in urban areas of the lower and middle watershed, where the sites often had vegetated riparian buffers (for instance, IMC4 scored 42, Table 1). The proportion of fine substrate was greater in headwaters though not significantly (r = −0.679, p = 0.110). Fine substrate proportion appeared to be related to the land use near sites, with agricultural sites IMC1 through IMC4 having finer substrate (Fig. S3d). The prevalence of woody debris tended to be greater toward the mouth of the creek (r = 0.703, p = 0.078; Fig. S3e), as did stream velocity (r = 0.75, p = 0.066; Fig. S3f). However, no sites met the representative condition for Midwest streams of 32.6 woody pieces per 100 m, which was the mean density for Midwest sites found by Cordova et al. (2007) (Table 1). Reduced large woody debris was most evident in the agricultural headwaters where few wood pieces were found (Table 1). Though still below reference conditions, woody debris were more abundant in the urbanized, lower watershed, where 14–30 pieces per 100 m were found at sites (Table 1). This could be at least partly explained by small sections of forest along the stream in the urban landscape. However, relationships between downstream distance and suspended sediment (r = 0.11, p = 0.840; Fig. S3g) and bedload sediment (r = 0.61, p = 0.167; Fig. S3h) were less clear and were affected by anomalous measurements at sites. Water temperatures in the middle watershed were above the optimal brown trout limit of 19 °C (Raleigh et al. 1984) for nearly 550 h (37.0%) throughout July and August, 2017 (Table 1).

Macroinvertebrate metrics

The integrity of stream macroinvertebrate communities was related to the environmental variables in the watershed, suggesting that there were clear impacts of land use on aquatic life, although these impacts would not always accumulate in a downstream direction. Woody debris prevalence was positively associated with abundance (r = 0.69, p = 0.038) and richness (r = 0.77, p = 0.015) of EPT macroinvertebrates (Fig. 3a, b), though this was influenced by a large proportion of EPT (30%) at the IMC7 site. This suggests that a larger study with more sample sites would aid in more precisely understanding this relationship. Fine sediment in the streambed was associated with a larger proportion of Dipteran (r = 0.55, p = 0.121) and surface air breathing (r = 0.80, p = 0.014) taxa (Fig. 3c, d), although influenced by a high percentage of Diptera (55%) at the WAD site. It is unlikely that this site had an overly large effect on the correlation test results though because we used the nonparametric Spearman’s rank correlation. Good riparian condition was positively associated with abundance (r = 0.34, p = 0.377) and richness (r = 0.55, p = 0.125) of EPT, although not significantly (Fig. 3e, f). Thus, riparian conditions could be affecting the integrity of macroinvertebrate communities in the watershed and a larger study with more sample sites could better elucidate these relationships. Also, stream slope did not have a statistically significant relationship with macroinvertebrate metrics, such as EPT richness (r = 0.40, p = 0.28) and percent Diptera (r = −0.08, p = 0.85) (Fig. 3g, h), or other variables in the study.

Fig. 3

Relationships between environmental and macroinvertebrate metric variables. EPT: Ephemeroptera, Plecoptera, and Trichoptera macroinvertebrates. a Woody debris and percent EPT, b woody debris and EPT richness, c percent fine substrate and percent air breathers, d percent fine substrate and percent Diptera, e Procedure 51 riparian score and percent EPT, f Procedure 51 riparian score and EPT richness, g slope and EPT richness, h slope and percent Diptera

The ordinations demonstrated interrelationships between macroinvertebrate metrics and environmental variables that can occur along the land use gradient in the watershed. The first two CCA axes (CCA1 and CCA2) of the macroinvertebrate metrics ordination explained 52.9% and 15.0% of the total variation in the data (Fig. 4). The first CCA axis represented a gradient of substrate size and macroinvertebrate community integrity. Positive values of CCA1 corresponded to fine substrate associated with increased numbers of Diptera and air breathers that had poor function per Procedure 51 scores, while negative values of CCA1 corresponded to a high richness of Ephemeroptera and Trichoptera taxa (Table 2, Table S6), as well as increased amounts of woody debris and riparian vegetation (Table 1). The second CCA axis largely represented a flow gradient. Negative values of CCA2 corresponded to fast flowing habitats moving more suspended sediment, while positive values of CCA2 corresponded to slow flowing habitats such as pools. The remaining constrained axes explained 15.0% of the total variation and three residual CA axes explained 17.2%. Sites in the CCA are presented in Fig. S4 for reference. PERMANOVA results showed that the constraining variables from the macroinvertebrate metrics CCA of pools, fine substrate, and suspended sediment were significant, while the woody debris and riparian condition variable and water velocity were only slightly above α=0.05, which we deemed acceptable because of the small sample size and complex ecological conditions (Table 5). PERMANOVA also showed that overall the model was significant (df=5, F=6.31, p = 0.008). NMDS also showed an association between Ephemeroptera and woody debris, Trichoptera with high velocity and suspended sediment, and Dipterans and air breathers with fine substrate and pool habitat (stress = 0.028) (Fig. S5).

Fig. 4

Canonical correspondence analysis (CCA) ordination showing the relationships between environmental stressors and macroinvertebrate metrics in the Indian Mill Creek watershed, Michigan, USA (2017). EPT Ephemeroptera, Plecoptera, Trichoptera

Table 5 Permutational multivariate analysis of variance (PERMANOVA) results for macroinvertebrate metrics and traits analyses, with significant findings (α = 0.05) in bold at Indian Mill Creek (IMC), Michigan, USA

The integrity of aquatic macroinvertebrate communities did not decrease in a downstream direction as we had expected. The greatest richness of EPT (4 taxa at IMC5 and 3 taxa at IMC4) occurred in the urban zone, downstream of agricultural impacts, in an area with forested riparian buffers and prevalent woody debris (Table 2). The greatest percent of EPT (30%) was also within the urban zone, at the IMC7 site closest to the mouth of the creek, while the greatest Procedure 51 score was also at IMC5 within the forest-buffered section of the urban area. This suggests that other factors, such as the relationship with woody debris prevalence, are improving macroinvertebrate communities despite the accumulated agricultural and urban land use impacts.

Macroinvertebrate traits

In Indian Mill Creek, linkages of environmental variables with macroinvertebrate feeding groups and traits often occurred along the agricultural-to-urban land use gradient, suggesting that relationships occur beyond those expected of an idealized natural stream. This gradient represents a transition from primarily agricultural land in the upper watershed (67% of area upstream of Interstate 96 is agricultural) to urban development in the lower watershed (78% of area downstream of Interstate 96 is urban; Fig. 1; Homer et al. 2015). Scrapers (e.g., Physidae, Heptageniidae) and herbivores (e.g., Corixidae) were most commonly found in the upper watershed, with particularly high abundance at the IMC2 site, and were associated with poor riparian condition (r = 0.57, p = 0.109 for scrapers; r = 0.64, p = 0.062 for herbivores) (Table 3, Fig. 5a, b). The particularly high abundance of scrapers at IMC2 (138) suggests that a larger study with more sites could better elucidate these relationships and supports our use of the nonparametric Spearman’s rank correlation. Scrapers and herbivores were rarely found in the lower, urbanized areas (Table 3). For instance, the IMC7 site only had one scraper (Physidae; Table S5 and S6). Collector gatherers (e.g., Baetidae, Turbellaria) and filterers (e.g., Hydropsychidae) were found throughout the watershed but were most commonly found in the middle and lower reaches (e.g., 190 collector gatherers at IMC6, 85 collector filterers at IMC7; Table 3). They were associated with high water velocity (r = 0.58, p = 0.108 for filterers) and good riparian conditions (r = 0.63, p = 0.070 for gatherers) (Fig. 5c, d), though this relationship is affected by the large number of filterers at IMC7. Collectors were not as commonly found in agricultural headwaters (e.g., only 16 collector gatherers and no collector filterers were found at IMC2; Table 3). Climbers (e.g., Aeshnidae, Coenagrionidae) and swimmers (e.g., Baetidae, Hydrophilidae) were associated with pools of the two tributary sites, for instance there were 214 swimmers in the fixed sample for the WAD site (Table 3). Additionally, predators (namely Veliidae) comprised 80% of the macroinvertebrate community at the Brandywine Creek (BC) site (Table 3, Table S6). Shredders and sprawlers (e.g., Isopoda) were found throughout the watershed but were most associated with high velocity habitats (r = 0.59, p = 0.092 for sprawlers) and coarse substrate (r = 0.60, p = 0.097 for shredders) (Fig. 5e, f). These conditions were prevalent in the urbanized middle and lower watershed (Table 1). Shredders and sprawlers were not commonly found at the agricultural site IMC2 or the tributary sites, which had a high proportion of fine streambed substrate (Tables 1 and 3).

Fig. 5

Relationships between environmental and macroinvertebrate trait variables. a Procedure 51 riparian score and scrapers count, b Procedure 51 riparian score and herbivores count, c Procedure 51 riparian score and gatherers count, d water velocity and filterers count, e water velocity and sprawlers count, and f percent fine substrate and shredders count

The ordinations demonstrated how environmental conditions can affect the traits of aquatic macroinvertebrate communities along the agricultural-to-urban land use gradient. The first two CCA axes (CCA1 and CCA2) of the macroinvertebrate feeding and traits ordination explained 52.4% and 18.8% of the total variation (Fig. S6). The first CCA axis represented an aquatic habitat structure gradient. Negative values of CCA1 corresponded to well-structured habitats with woody debris and riparian vegetation buffers, while positive values of CCA1 corresponded to poorly structured habitats with lots of fine substrate. The second CCA axis appeared to show a flow gradient. Negative values of CCA2 corresponded to faster-flowing habitats with high water velocities and suspended sediment, while positive values of CCA2 corresponded to slower-flowing habitats with more pools. The remaining constrained axes explained 11.3% of the total variation and three residual axes explained 17.6%. PERMANOVA results from the macroinvertebrate traits CCA constraining variables showed that fine substrate and suspended sediment were significant, while the woody debris and riparian condition variable was only slightly above α=0.05 (Table 5). PERMANOVA results also showed that overall the macroinvertebrate traits model was significant (df=5, F=4.07, p = 0.012). NMDS showed an association between climbers and swimmers with pools, collector filterers with water velocity and suspended sediment, and predators with fine substrate (stress = 0.041) (Fig. S7).

Fish

Fish survey results revealed how environmental conditions can affect fish communities in the watershed. We had low catch numbers at all sites (25 to 78 individuals, mean=42), though the surveys were confined to priority sites in the lower and middle watershed (Table 4, S7). Fish numbers were particularly low in the lower, urban reaches (25–26 individuals) and increased in an upstream direction to 78 individuals at the IMC4 site.

Fish community assemblage appeared to be largely structured by stream temperature and flow. Salmonids, white sucker (Catostomus commersonii), and bluegill (Lepomis macrochirus) were associated with stable, coldwater reaches, while small Cyprinids and Johnny darter (Etheostoma nigrum) were associated with higher temperature reaches (Table 1, S7). For instance, the IMC7 site had a colder average temperature of 16.2 °C and 69% of the catch salmonids, while the IMC4 site had a warmer average temperature of 18.3 °C and 53% of the catch creek chubs (Semotilus atromaculatus). The low-head dam (Fig. 1) is located between sites IMC6 and IMC7, which had the highest proportion of salmonids (37.1% and 69.2%, respectively). Each site had a poor Procedure 51 fishery score, either due to low values or insufficient catch (Table 4). Few fish were found in the middle and lower watershed sites (for instance, 25 individuals upstream of the dam), where pools occupy as low as 14–17% of the habitat area (Tables 1 and 4). The 3 Mile Road and IMC4 sites, which were the furthest upstream sites, contained the largest number of individual fish (66 and 78, respectively) and also had high proportions of small minnow and darter species (Table 4, Table S7).

The NMDS biplot of the fish community showed three distinct groups that appeared to be influenced by velocity and temperature regimes (Fig. 6). The MRPP confirmed our visual interpretation of differences among the groups (A = 0.278, p = 0.009). The MRPP comparisons between groups revealed that there was no difference between the fish communities in the fast and slow flow reaches (A = 0.288, p = 0.33). However, communities did appear to differ as a function of temperature regime with marginally significant differences between the fast velocity and warm temperature communities (A = 0.23, p = 0.10) and between the slow velocity and warm temperature communities (A = 0.22, p = 0.10).

Fig. 6

Nonmetric multidimensional scaling (NMDS) plot of fish communities in three distinct groups driven largely by stream temperature and flow regime in the Indian Mill Creek (IMC) watershed, Michigan, USA (2017; A = 0.278, p = 0.009). U/S upstream of dam

Discussion

Environmental variables and aquatic communities

Landscape

In the Indian Mill Creek watershed, longitudinal patterns are confounded with a gradient of agricultural headwaters to urban areas in the lower reaches. Indian Mill Creek is much smaller (third order) than streams typically evaluated with the River Continuum Concept and subject to significant and widespread anthropogenic impacts. We found that macroinvertebrates of the collector feeding group, such as Hydropsychid caddisflies and Turbellarians (Table S5), were associated with aquatic communities throughout the watershed, as predicted for first- to third-order streams like Indian Mill Creek by the River Continuum Concept (Vannote et al. 1980). They were most abundant in the middle and lower reaches, associated with high water velocity and good riparian conditions. Collectors were not as common in agricultural headwaters which generally had degraded riparian conditions. Additionally, in an idealized stream with no anthropogenic impacts, the River Continuum Concept predicts that shaded streams (such as headwaters) will have low numbers of macroinvertebrate scrapers because of low periphyton abundance (Vannote et al. 1980). We found that the scraping and herbivorous macroinvertebrates were most common in the upper watershed, particularly the IMC2 site, and were associated with poor riparian condition, supporting that Indian Mill Creek would not conform to the idealized natural stream. This pattern could be from fertilizer runoff entering the stream and increasing periphyton growth (Compin and Céréghino 2007), coupled with the absence of shading from riparian vegetation that would otherwise reduce periphyton abundance in these agricultural areas (Wooster and DeBano 2006). Shredding macroinvertebrates were associated with middle reaches that had extensive riparian vegetation buffers, variable habitats, and likely more inputs and retention of leaves and other organic matter. This relationship is similar to the predictions of the River Continuum Concept, although the location of these shredder-abundant reaches were rearranged from those of an idealized natural stream (where they would be in the headwaters; Vannote et al. 1980) because of the anthropogenic impacts. Also, similar to previous research (Baumgartner and Robinson 2017), macroinvertebrate taxa richness, particularly of EPT taxa which are generally more sensitive to environmental degradations, was lowest at the agricultural sites of upper Indian Mill Creek (e.g., IMC1 and IMC2).

Our slope variable did not have statistically significant relationships with other variables in the study, which aligns with previous findings that intensive land uses can homogenize aquatic communities so they no longer follow a longitudinal gradient (Delong and Brusven 1998). Our findings that aquatic macroinvertebrate communities (e.g., EPT richness) did not have a significant relationship with stream slope, but rather were strongly influenced by land use such as riparian management, support this. Fish communities appear to be associated with geomorphology (Section 4.1.2) or water quality and pollution (Section 4.1.3) as well.

Geomorphology

The distribution of streambed substrate in a watershed, such as the one in this study, can be explained by a combination of geomorphology as well as the agricultural-to-urban land cover gradient. Lane’s Balance (Lane 1955; Dust and Wohl 2012; Pollock et al. 2014) conceptualizes that water discharge and channel slope are related to sediment load and representative particle size and shows whether aggradation or degradation will occur under changing scenarios. The coarsest substrate in Indian Mill Creek was in the middle reaches as the creek descends the Grand River Valley. Here, a steep gradient (Fig. 2) tipped Lane’s Balance toward increased particle size and erosion of fine particles. The agricultural, upper reaches had more gradual slopes, tipping the balance toward finer substrate and sediment deposition.

Additionally, our results show that bedload sediment was highest in the urbanized, lower watershed that had softer elevation changes (Fig. 2). Here, channelization and high flows from storm water runoff (Walsh et al. 2005) likely counteract the low slope, increasing velocity and tipping Lane’s Balance toward coarser substrate, while also moving bedload sediment. This is connected to our objectives because the bedload sediment is likely reducing habitat variability and we found both macroinvertebrates and fish to be related to flow regimes. Excessive bedload sediment can homogenize the stream channel and create a long run with uniform depth and velocity instead of a riffle/pool series (Alexander and Hansen 1986), which would affect the low habitat variability we witnessed toward the mouth of the creek. Reduced habitat variability from high bedload sediment can be unsuitable for macroinvertebrates and fish such as brown trout, which have optimal conditions of 30–50% riffle and 50–70% pool area (Raleigh et al. 1984; Alexander and Hansen 1986; Culp et al. 1986; Sigdel 2017), and which were present in our study. Our results show that in the Indian Mill Creek watershed, the coverage of riffle and pool habitats is almost always lower than these optimal conditions, which is likely affecting the fish community as well as the macroinvertebrates. Also, the fish community in the small dam impoundment was more similar to the community in the slow-flowing, lower stream reach, near the confluence with the Grand River, suggesting that the dam may be artificially affecting the fish community (Table S7). These geomorphic patterns in streams with this agricultural-to-urban land use gradient are vital to watershed managers studying similar catchments because they can distinguish geomorphic relationships in streambed substrate from those related the management of sediment from urban stormwater, agricultural erosion, or streambank erosion (Paul and Meyer 2001; Allan 2004; Kiesel et al. 2009).

Water quality and pollution

The fish catch in urban reaches of the watershed was low in regards to other Indian Mill Creek sites and Procedure 51 assessments in nearby Michigan streams and was a cause for “poor” ratings (EGLE 2005, 2011, 2016). It is uncommon for total catch to be this low at streams in the region (EGLE 2005, 2011, 2016). For instance, Procedure 51 surveys of nearby Deer Creek and Plaster Creek in 2014 had total catch of 144 and 436 fish, respectively (EGLE 2016). Also, since the Procedure 51 fish survey tests for there being greater than 50 fish caught, all sites with total catch below 50 individuals were automatically rated as “poor” (EGLE 2008).

The low fish catch in comparison to other nearby Michigan streams studied using the Procedure 51 survey (EGLE 2005, 2011, 2016) warrants additional discussion of possible causes, including episodic urban pollution events and toxicity. In particular, the absence/low numbers of small-bodied, sedentary fishes such as sculpin, darters (Etheostoma spp.), and small minnows (Cyprinidae spp.) in the lower reach of Indian Mill Creek was peculiar and may be a result of the episodic, urban pollution events (Table S7). The episodic pollution events (Section 2.4) and toxic sediments (Parker 2018) could at least partially explain why our total fish catch was low at the lower urban sites, particularly for the small-bodied, sedentary fishes. Episodic pollution events can have major impacts to fish communities via direct mortality or by causing fish to seek refuge in unpolluted waters (Seager and Maltby 1989; Van Sickle et al. 1996), which could explain why these smaller fish with less mobility were in low abundance toward the mouth of Indian Mill Creek. Small fish such as darters (Mundahl and Ingersoll 1983), small-bodied Cyprinids (Mundahl and Ingersoll 1989), and sculpin (Breen et al. 2009) tend to have relatively small home ranges in stream systems. Full recovery of small fish populations from these events can take several years (Albanese et al. 2009; Kubach et al. 2011). Slow recovery rates can be further exacerbated when dispersal barriers are present (Albanese et al. 2009) such as the low-head dam in the lower reach of Indian Mill Creek. This is important for our objectives of understanding how environmental variables affect aquatic communities in Indian Mill Creek because it suggests that episodic pollution events could be a cause of the low fish catch numbers in the lower watershed sites.

Previous research has found aquatic macroinvertebrate communities to be related to water quality and pollution variables in streams with different land uses (Berger et al. 2017). These variables include herbicides, pesticides, electrical conductivity, dissolved oxygen, silicate, and caffeine (an indicator of untreated wastewater pollution) (Berger et al. 2017). Notable among these variables are agricultural pesticides, such as insecticides and fungicides, which widely impact macroinvertebrate communities in streams (Schäfer 2019). However, Parker (2018) did not detect these pesticides in Indian Mill Creek (Section 2.4). This suggests that studies focusing only on these water quality and pollution variables may omit stream habitat variables, which were an important part of our study, and thus could also miss relationships. It is best for future studies to include a diverse set of water quality and pollution variables (such as pesticides; Schäfer 2019) and aquatic habitat variables.

Longitudinal and land use gradients

A major way that the longitudinal processes in our study stream differ from those in idealized natural streams is deforestation in the agricultural headwaters. Our results of low woody debris abundance in the upper agricultural watershed show that deforested riparian zones and lack of mature trees in the agricultural headwaters likely limit the amount of large woody debris input to upper Indian Mill Creek. This could be because urban and agricultural streams are typically channelized and cleared of debris to enhance water conveyance from upland sources and to prevent flooding (Booth et al. 1996; Johnson et al. 2003). Woody debris showed an overall increase in prevalence from the upper to lower reaches of Indian Mill Creek but were most prevalent in the middle watershed, which had forested riparian areas along the valley slope. This pattern could be explained by recruitment of woody debris from bank erosion and retention by debris jams from fallen trees (Martin and Benda 2001). The prevalence of woody debris could also explain the formation of pools at these sites (Evans et al. 1993; Martin and Benda 2001). In these instances, the riparian management thus has a major impact on aquatic habitats. Woody debris are important because they provide a mechanism for habitat, refuge, and substrate for biofilm food production for aquatic macroinvertebrates (Hax and Golladay 1993; Johnson et al. 2003), which explains why they were so positively related to abundance and richness of EPT macroinvertebrates in our study. This demonstrates the importance of woody debris and riparian buffers as management tools for streams following this agricultural-to-urban land use gradient.

The influence of land use (described in Section 4.1.1) and riparian management leads us to suspect that the agricultural-to-urban land use gradient has a major influence on patterns of environmental variables and aquatic communities in the creek, compared with the longitudinal patterns expected of an idealized, natural stream (Vannote et al. 1980). This was also witnessed by relationships between distance downstream with riparian scores and woody debris prevalence, which are likely attributed to the changing land use patterns rather than longitudinal processes. Changes in land use can cause physical disturbance regimes that alter geomorphic and hydrologic processes in watersheds, affecting their habitat characteristics and ecological structure and function (Poff et al. 2006).

Locations of greatest impacts

Along the agricultural-to-urban land use gradient, there is a series of environmental variables that structure biological communities in the creek (diagram in Fig. S8), causing impacts to aquatic communities to not always increase in a downstream direction through the agricultural and urban developments. Understanding the relationship of these multiple variables within these stream ecosystems is important for successful assessment and restoration of watersheds around the world with this common agricultural-to-urban land use pattern. For example, while previous findings in the lowland Swiss Mönchaltorfer Aa catchment showed that EPT richness was highest in least impacted sites mostly upstream of urban and agricultural areas, our greatest EPT richness was along sites IMC4 and IMC5, which are in well-forested reaches downstream of the agricultural areas but within the urban zone (Baumgartner and Robinson 2017). Also, contrary to findings from Brazil’s Uberabinha river catchment and Chile’s Mediterranean ecoregion, we did not find urban reaches to have less richness or diversity of sensitive taxa than agricultural reaches in the Indian Mill Creek watershed (Nascimento et al. 2018; Fierro et al. 2019). This divergence from other catchments suggests that one possible reason for this increase in EPT richness is that there is a forested section of the middle Indian Mill Creek watershed within the urban zone, which provides habitats and shelter for aquatic organisms (Fig. 1). This finding, which we did not hypothesize to occur, is important for researchers and watershed managers because it shows promise for streams to recover from agricultural degradations upstream. It is reasonable because the impacts of urbanization can be weaker in streams already degraded by agriculture in the watershed (Cuffney et al. 2010). This particularly supports the idea that effective management of the riparian zone could improve aquatic communities in vulnerable areas of watersheds with this common agricultural-to-urban land use gradient and is an important take-away from the study. Notably, it supports that small patches of forest can be key to conserving aquatic biodiversity in urban landscapes.

Limitations and future research

Our findings show that abundant woody debris, riparian vegetation, and pool habitats can occur together. This is important because it suggests that relationships among environmental variables can often co-occur, confounding analyses of their influences on aquatic communities. For instance, a study in German lowland rivers found that fine sediment cover and high phosphorus concentrations often co-occur (Lemm and Feld 2017). These co-occurring relationships could cause difficulties when trying to narrow focus on a single environmental variable that is impacting aquatic communities. However, a study needing to reduce the number of habitat variables, for instance because of inventory costs, could do so knowing these relationships. Additionally, difficulties with watershed studies such as this include that sampling sites are not completely independent of each other as they are in the same catchment, which could violate the assumptions of statistical analyses. This highlights the need for investigation of spatial autocorrelations when developing models (Miralha and Kim 2018). It is also possible that there was temporal variability between sampling dates not accounted for and seasonal or annual variability in aquatic conditions that could affect the transferability of our results to other seasons or years.

Future studies could be improved by replicating our approach in a suitable reference catchment with no significant agricultural or urban landcover, permitting a clear separation of impacts of agricultural-to-urban gradients from natural longitudinal gradients in streams. Researchers should also further examine the co-occurrence of environmental variables, such as woody debris, riparian vegetation, and pool habitats and provide guidance for indicator variables that can reduce the costs and time of aquatic resource inventories. Additionally, although we identified macroinvertebrates only to the family level per the rapid bioassessment protocol, generic-level identification used in most similar studies would provide a more thorough analysis and likely additional insights and confidence into the relationships between environmental variables and aquatic macroinvertebrate communities. Further, steam size could be another environmental variable to investigate.

Conclusion

For river basins following a land use pattern of agricultural headwaters to urban lower reaches, interactions of landscape, geomorphic, and water quality/pollution variables can be influential on the structure of fish and macroinvertebrate communities. For instance, the integrity of stream macroinvertebrate communities, as well as the prevalence of macroinvertebrate scrapers and shredders, was associated with surrounding land use and riparian management. Also, lower fish abundance and community integrity in the lower urban reaches of Indian Mill Creek could be explained by a combination of geomorphology and episodic urban pollution events in addition to flow and temperature regimes. Next, we found that aquatic communities would not always conform to expectations for longitudinal gradients like the unimpacted natural stream from the River Continuum Concept. Deforestation and clearing of the riparian zone in agricultural headwaters caused there to be more scraping macroinvertebrates than we would have expected in a natural stream without these impacts. We also found that impacts would not always increase in a downstream direction through the agricultural and urban developments. This is because the abundance and richness of sensitive EPT macroinvertebrates increased within the urban zone of the Indian Mill Creek watershed, downstream of extensive agriculture, within a forested riparian buffer. This shows that recovery of degraded aquatic communities downstream of urban or agricultural impacts is possible given effective management choices, such as maintaining small patches of forest along the streams. Agricultural and urban land cover changes and their associated impacts to lotic ecosystems are prevalent and increasing worldwide. This study has provided a better understanding of these impacts along streams following an agricultural-to-urban land cover gradient, which may not conform with general expectations of land use patterns. It has also provided guidance for successful mitigation of these stressors.