Many physicochemical factors potentially impair stream ecosystems in urbanizing basins, but few studies have evaluated their relative importance simultaneously, especially in different environmental settings. We used data collected in 25 to 30 streams along a gradient of urbanization in each of 6 metropolitan areas (MAs) to evaluate the relative importance of 11 physicochemical factors on the condition of algal, macroinvertebrate, and fish assemblages. For each assemblage, biological condition was quantified using 2 separate metrics, nonmetric multidimensional scaling ordination site scores and the ratio of observed/expected taxa, both derived in previous studies. Separate linear regression models with 1 or 2 factors as predictors were developed for each MA and assemblage metric. Model parsimony was evaluated based on Akaike's Information Criterion for small sample size (AICc) and Akaike weights, and variable importance was estimated by summing the Akaike weights across models containing each stressor variable. Few of the factors were strongly correlated (Pearson |r| > 0.7) within MAs. Physicochemical factors explained 17 to 81% of variance in biological condition. Most (92 of 118) of the most plausible models contained 2 predictors, and generally more variance could be explained by the additive effects of 2 factors than by any single factor alone. None of the factors evaluated was universally important for all MAs or biological assemblages. The relative importance of factors varied for different measures of biological condition, biological assemblages, and MA. Our results suggest that the suite of physicochemical factors affecting urban stream ecosystems varies across broad geographic areas, along gradients of urban intensity, and among basins within single MAs.
A growing body of research has improved our understanding of how lotic ecosystems are affected by urban development in the watershed. Several generalizations have emerged from research in streams influenced by urbanization (Paul and Meyer 2001, Walsh et al. 2005). First, ecological responses to urbanization are strongly influenced by the delivery of runoff to stream channels, which in turn is controlled by variation in impervious land cover and stormwater infrastructure. The chemical and hydrological characteristics of this runoff have profound effects on the flow regime, habitat structure, and biological processes of urban streams. Second, ecological responses to increasing urbanization are profound and include measurable changes to the chemical and physical environment. Last, dramatic losses of biological diversity or changes in assemblage composition and structure that probably impair ecosystem functioning generally occur in urban streams.
Our understanding of the potential causes of ecological impairment in urban streams is incomplete for several reasons (Paul and Meyer 2001, Walsh et al. 2005). First, factors that occur sporadically or for short periods of time (e.g., organic contaminants) are often inadequately characterized by single samples, which tend to underestimate their importance to biological impairment. Second, studies often are confined to specific geographic areas, which limits understanding of how the potential biological effects of factors vary among different environmental settings (e.g., climates). For example, salinity from road-salting operations is a potentially important factor in cold northern climates but not in warm southern climates (Kaushal et al. 2005), and broad-scale geographic differences in urban paving operations may influence the potential toxicity of runoff from asphalt surfaces (Scoggins et al. 2007). Third, biological endpoints are often limited to measures of a single assemblage, an approach that probably underestimates the ecological consequences of urbanization because different assemblages probably respond in unique ways (e.g., Carlisle et al. 2008) to the mélange of physicochemical factors present in urban streams. Last, untangling the potential causes of ecological impairment in urban streams requires experimentation at large (e.g., subbasin) scales (Walsh et al. 2005) and carefully designed observational studies that control for confounding natural factors and potential covariance among anthropogenic factors.
Our primary objective was to quantify the relative importance of several physicochemical factors in explaining variation in the condition of fish, macroinvertebrate, and algal assemblages in each of 6 different metropolitan areas (MAs) of the USA. The composition of individual biological assemblages was associated with urban development in several different MAs (Brown et al. 2009, Cuffney et al. 2010), and these associations differed among regions with different climate and antecedent land uses (Cuffney et al. 2010, Kashuba et al. 2010). Our goal was to identify which of a complex set of potential physicochemical factors were most closely related to condition of multiple biological assemblages. We identified a subset of 11 physicochemical factors (Table 1) from previous studies that were influenced by anthropogenic activities in urban settings. We used a model-selection approach with these physicochemical and associated biological data to quantify how much variation in biological condition could be explained by these factors as well as the relative importance of each factor.
Table 1.
Selected physicochemical factors measured along gradients of urbanization in 6 metropolitan areas and used as predictors. Variable codes from Giddings et al. (2009) and used in Cuffney et al. 2010, Fitzpatrick and Peppler 2010, and Steuer et al. 2010 are given in square brackets. HOC = hydrophobic organic contaminants, SPMD = semipermeable membrane device, PTI = pesticide toxicity index.
Methods
Study design
In 1999, the US Geological Survey (USGS) began an integrated assessment of the effects of urbanization on stream ecosystems in selected metropolitan areas of the USA (Tate et al. 2005). Streams were sampled along gradients of urban landuse intensity in and around 9 metropolitan areas (MAs), 6 of which are represented in this study: Atlanta, Georgia (Gregory and Calhoun 2007); Raleigh-Durham, North Carolina (Giddings et al. 2007); Denver-Fort Collins, Colorado (Sprague et al. 2006); Dallas-Fort Worth, Texas (Moring 2009); Milwaukee-Green Bay, Wisconsin; and Portland, Oregon (Waite et al. 2006). A multivariate index that represents different attributes of urbanization (e.g., infrastructure, population characteristics, land cover) was used to characterize the intensity of urbanization (McMahon and Cuffney 2000, Tate et al. 2005) for a population of candidate basins (2nd–5th-order streams) within each MA. From the candidate list, 25 to 30 streams were selected after field reconnaissance in each MA to: 1) avoid point discharges, such as sewage treatment plants or industrial discharges; 2) minimize natural variability resulting from factors such as ecoregion, slope, riparian conditions, and basin size; and 3) cover the range (gradient) of urban landuse intensity as characterized by the multivariate index.
Physicochemical factors
Several major physicochemical factors in streams, including hydrology, water temperature, channel morphology, substrate particle size, nutrients, ions, metals, pesticides, and hydrophobic organic contaminants (HOC; particularly polycyclic aromatic hydrocarbons), are influenced by urbanization (Paul and Meyer 2001, Walsh et al. 2005). We selected 11 factors (Table 1) that were quantified in the USGS urban studies and that collectively represent major categories of stressors (except metals). The derivation of each environmental variable and data for each site is reported in Giddings et al. (2009). Based on the findings of previous analyses of this data set (Giddings et al. 2009, Cuffney et al. 2010, Fitzpatrick and Peppler 2010, Steuer et al. 2010), we selected a subset of physicochemical factors that were associated with urban land cover and various measures of biological assemblage composition and structure.
Continuous stream-stage data were collected at or near the biological sampling reach for ∼1 y during the biological data collection effort. Stream-stage values collected at 15-min intervals were converted subsequently to cross-sectional area, and hydrologic metrics were calculated based on hourly changes in channel cross-sectional area (McMahon et al. 2003). Water temperature was measured continuously every 15 min at each site and converted to an hourly time step prior to calculation of annual degree days (Cuffney and Brightbill 2008). Channel geometry and streambed substrate data were collected at 11 equally spaced transects along the sample reach (Fitzpatrick et al. 1998). Morphologic indicators were used to estimate bankfull stage. Bankfull width/depth ratios were calculated for each transect based on cross-sectional profiles. Coverage of fine-particle substrates was estimated at 3 points on each of the 11 transects and averaged across the entire reach. Stream water-chemistry data were from a single collection under similar hydrologic conditions across sites within an MA (base flow) from equal-width increments across the stream channel in accordance with standard National Water Quality Assessment (NAWQA) Program protocols (Shelton 1994). Water samples were analyzed for nutrients, ions, and dissolved pesticides according to standard protocols (Fishman 1993, Zaugg et al. 1995) at the USGS National Water Quality Laboratory in Denver, Colorado. The potential additive toxicity of dissolved pesticide mixtures was characterized with a pesticide toxicity index (PTI; Munn et al. 2006). The PTI combines information on exposure of aquatic biota to pesticides (measured concentrations of pesticides in stream water) with toxicity values (from laboratory bioassays) to produce a relative index value for a water sample. Therefore, PTI provides a means of comparing the potential toxicity of pesticides across streams with varying compound mixtures and concentrations. Separate PTIs based on toxicity to macroinvertebrates or fish were calculated and used in the corresponding model. The lack of toxicity data for algal species required that we use the sum of herbicide concentrations in each sample as a simplistic measure of potential toxicity in models of algal response.
Semipermeable membrane devices (SPMDs) were used to concentrate HOCs as an alternative to the collection and analysis of water, bed sediment, or biological tissue samples (Huckins et al. 1990). SPMDs were deployed 4 to 6 wk prior to sampling biological assemblages. Detailed information about deployment, retrieval, and analysis of SPMD extracts can be found elsewhere (Bryant et al. 2007, Bryant and Goodbred 2008). SPMDs provide data relevant for assessment of exposure and potential toxicity and offer a distinct advantage over sampling other media by accumulating the bioavailable fraction of waterborne HOCs, many of which may be present at ecologically relevant levels in water or other ambient media but below the limits of detection for standard analytical methods (Gourlay et al. 2005, Huckins et al. 2006). The P450RGS assay was used to evaluate the potential and relative toxicity of the complex mixtures sequestered in the SPMDs. The P450RGS assay provides an integrated (toxicity × concentration) measure of relative toxicity from compounds that are aryl-hydrocarbon-receptor (AhR) agonists. These compounds include polychlorinated biphenyls (PCBs), polycyclic aromatic hydrocarbons (PAHs), dioxins, and furans. The P450RGS assay is used as an analytical chemistry tool to measure the levels of AhR-active compounds (expressed in toxic equivalents), but more importantly, it also is used as a model of organism response to these compounds. The P450RGS assays were conducted at the US Army Corps of Engineers Research and Development Center in Vicksburg, Mississippi, with the protocols of Ang et al. (2000).
Assessing biological condition
Algal, macroinvertebrate, and fish assemblages were sampled using the methods described in Moulton et al. (2002). Quantitative algae samples were scraped from hard substrates, which were either rocks collected from riffles or woody snags in the sampling reach. The sampling area was estimated either by using pieces of aluminum foil fitted to the scraped rocks and later traced onto graph paper or by using calipers to measure diameter and length of snags. Algae were identified and enumerated from permanent slides under 1000× power at the Patrick Center for Environmental Research (Academy of Natural Sciences, Philadelphia, Pennsylvania) following methods described in Charles et al. (2002). Macroinvertebrates were collected from 5 discrete riffle locations or woody snags with a Slack sampler (Moulton et al. 2002) and combined into a single composite sample. Samples were processed at the USGS National Water Quality Laboratory using the methods described in Moulton et al. (2000). Fish assemblages were sampled with a single backpack electrofishing unit from stream reaches that were between 150 and 300 m in length. Specimens collected from each of 2 passes were combined into a single sample. Fish generally were identified and enumerated in the field, but specimens that could not be identified were retained for determination and enumeration in the laboratory (Walsh and Meador 1998).
We characterized biological assemblages with 2 metrics: one that describes changes in relative abundances and species composition along the urban gradient within each MA (nonmetric multidimensional scaling [NMDS] ordination scores), and one (ratio of observed/expected taxa [O/E]) that quantifies deviation in each site's assemblage composition from that of regional reference sites. We chose these 2 metrics because they are derived from the entire species assemblage for each assemblage type, were associated with urbanization in many of the MAs investigated, and represent contrasting but complementary (i.e., compositional change vs deviation from regional reference) attributes of assemblage responses to anthropogenic disturbance. We did separate NMDS analyses for each MA following methods outlined in Brown et al. (2009) and detailed by Cuffney et al. (2010) and Kashuba et al. (2010). We derived 2 NMDS ordination axes (Clarke and Gorley 2006) based on 4√(x)-transformed abundance data and Bray–Curtis similarities. We rotated ordination axes (varimax) and relabeled them so axis 1 represented the greatest portion of the variance. We rescaled site scores to an origin = 0 within each MA and to a consistent direction of change with urban intensity across the MAs by subtracting from the minimum or maximum value (Cuffney et al. 2010). A potential limitation of using ordination scores as the biological response variable is that NMDS axes are arbitrary, and each is uniquely derived from the set of observations in the ordination (Clarke and Gorley 2006). In a regression context, a unit change in the axis in 1 MA is not comparable to a unit change in a different MA. Therefore, we avoided direct comparisons of responses (i.e., regression slopes) among MAs. However, in previous studies using this data set, Cuffney et al. (2010) reported that the axis scores were among the strongest indicators of change in biological assemblages (among 190 metrics investigated) along gradients of urban land cover. These findings suggest that although potentially arbitrary, these axes scores represent gradients of assemblage composition change associated with urban development, which we think justifies their use in our study.
We used an O/E index (sensu Hawkins 2006) that was developed for each assemblage in previous studies (Carlisle and Meador 2007, Carlisle and Hawkins 2008, Yuan et al. 2008, Meador and Carlisle 2009, Potapova and Carlisle 2011). The O/E index is the ratio of observed (O) and expected (E) assemblage composition, and hence, is a measure of the departure from expected natural conditions. The O/E index provides an intuitive measure of the degree to which an ecosystem supports the complement of species that should naturally occur at a site. The index is standardized by each site's potential because deviation of O from E is expressed as a ratio. Thus, the index is comparable across sites that differ naturally in their species composition and richness. With one exception (see below), we derived predictions of E for fish and macroinvertebrate assemblages from River InVertebrate Prediction and Classification System (RIVPACS)-type models (Moss et al. 1987, Hawkins 2006) that predicted the probabilities of capturing (Pcs) each species in the regional species pool at a given site. We estimated E for each site as the sum of predicted Pcs for common taxa (Pc > 0.50). We calculated O for each site as the number of expected taxa that were actually collected in the sample. We developed separate predictive models for macroinvertebrates from 338 reference sites in the eastern and central USA (Carlisle and Meador 2007), 217 reference sites in the central USA (Yuan et al. 2008), and 729 reference sites in the western USA (Carlisle and Hawkins 2008). We developed predictive models for fish assemblages from 266 reference sites in the eastern USA (Meador and Carlisle 2009).
For fish assemblages in the western USA and for diatom assemblages nationwide, we estimated E from the values of regional indices of assemblage composition and structure. We assessed fish assemblages in the western USA with an Index of Biotic Integrity (IBI) developed from 210 reference sites, where the IBI represented measures of assemblage composition other than species richness (e.g., proportion of exotic species; Whittier et al. 2007). O for fish assemblages in western sites was the observed value of the IBI calculated from the sample collected at each site (Meador et al. 2008).
We assessed algal assemblages with IBIs developed for diatoms in each of 5 geographic regions: Western Mountains, Central and Western Plains, Glaciated North, Southeastern Plains, and Eastern Highlands (Potapova and Carlisle 2011). The algal IBI was based on relative abundances of diatom taxa indicative of natural conditions at each site and on taxa indicative of disturbed conditions based on independent analyses (Potapova and Carlisle 2011). For fish in the Western Mountains and algae nationwide, we standardized the index value at each site (O) to the mean index value among reference sites ( = E) within each respective region to produce an index similar in interpretation to that from RIVPACS-type models (Hawkins 2006). We truncated all O/E values that were >1 to a value of 1 (Carlisle et al. 2008).
O/E values calculated with RIVPACS-type models and those calculated with index-based models have 2 important distinctions that may have influenced the regression models. First, index-based O/E values are ratios of variables measured on a continuous scale, whereas RIVPACS-based O/E values are ratios of common taxon counts, which are usually relatively small (e.g., <30) integers. This distinction can produce differences in the statistical behavior of the O/E values. However, our examination of the distributions of O/E ratios between the 2 approaches revealed similar ranges of values. A 2nd distinction between RIVPACS- and index-based O/E values is that the RIVPACS-type models attempt to account explicitly for natural among-site variation in E, whereas the index-based models assume a constant E for all sites within each region. As a result, index-based O/E ratios probably contain more unexplained natural variability than do RIVPACS-based O/E ratios, so use of index-based O/E ratios could lead to regression models with poorer explanatory ability.
Statistical methods
Interest in quantifying the relative importance of predictors in statistical models is growing, but currently, no method is generally accepted (Grömping 2007, Murray and Conner 2009), largely because disentangling the unique contribution of each predictor is not straightforward when predictor variables are correlated, as is often the case with empirical data. Few cases of strong collinearity (bivariate Pearson correlations on log(x)-transformed data, |r| > 0.70; Mason and Perreault 1991, Katz 2006) occurred among predictors in our data set (see Results), so we are confident that the ranks of predictor importance are interpretable, with a few exceptions. However, some authors (e.g., Murray and Conner 2009) have reported that even relatively low (e.g., |r| = 0.40) levels of correlation can obscure accurate identification of variable importance. Our goal was to use the principles of model selection to identify the most important factors by evaluating the evidence in the observed data for alternative models (Burnham and Anderson 2002) that represented different hypothesized effects of factors on biota. This approach seemed appropriate because the ecological responses to urbanization involve the action of multiple processes that are often best explained by multiple models (Hobbs and Hilborn 2006). The composite forces of measured and unmeasured factors in urbanizing streams (i.e., those responsible for observed phenomena) probably cannot be approximated by any single best model, especially given the many limitations (e.g., sample size) on model complexity. Therefore, observed data probably will provide evidence for several models, each of which represents a portion of the underlying forces at work. Thus, the process of selecting the most influential predictors requires consideration of the weight of evidence for the models in which they appear (Burnham and Anderson 2002).
We generated a separate set of regression models for each biological assemblage (× 2 indicators) and MA because we expected the factors associated with biological condition to vary among assemblages (Carlisle et al. 2008) and geographically (Cuffney et al. 2010, Kashuba et al. 2010). The number of observations (i.e., streams) within each MA ranged from 25 to 30, so we considered only models with ≤2 log(x)-transformed predictors (i.e., stressors) in simple linear regressions (Harrell 2001, Burnham and Anderson 2002). We considered all possible subsets (n = 66) of 1- and 2-variable models because no scientific basis existed for considering specific alternative sets of predictors a priori (Burnham and Anderson 2002). For example, no reason existed to propose a model that included pesticides and nutrients in lieu of models that included pesticides alone or with any other factor.
We calculated AICc (Akaike Information Criterion for small sample size) for each model within each model set (66 models for each assemblage × biological indicator × MA) and standardized it to the model with the smallest AICc value within that set. AICc is a measure of model quality that balances parsimony against quality of fit, and when standardized (i.e., ΔAICc), provides a quantitative ranking of each model. In general, models with ΔAICc values <2 are equally supported by the data (Burnham and Anderson 2002). We evaluated the plausibility of each model by calculating Akaike weights (w), which quantify the likelihood that a model is the best model, given the set of models under consideration and the data at hand (Burnham and Anderson 2002). We estimated the relative importance of each factor within each model set by summing the ws across all models in which that factor appeared (Burnham and Anderson 2002).
Results
Most physicochemical factors spanned large ranges in each MA (Table 2). MA-specific ranges in nutrient concentrations and high-flow duration varied from 10- to >100-fold. The ranges in high-flow magnitude, % fines, Cl−, and HOCs were generally >10-fold. The ranges in the frequency of high flows varied from 5- to 9-fold, whereas that of width/depth ratio varied from 2- to 5-fold. Ranges in stream temperature (degree days) were consistently <2-fold. Pesticide indices in all MAs varied from 0 (no detections) to ≥0.01. With 2 exceptions (Dallas-Fort Worth and Milwaukee-Green Bay MAs), total herbicide concentrations varied from 0 to ≥10 µg/L.
Table 2.
Ranges of factors measured along gradients of urban intensity in 6 metropolitan areas, including Atlanta, Georgia; Denver-Fort Collins, Colorado; Dallas-Fort Worth, Texas; Milwaukee-Green Bay, Wisconsin; Portland, Oregon; and Raleigh-Durham, North Carolina. Min = minimum, Max = maximum, HOC = hydrophobic organic contaminants, TEQs = toxic equivalents, herbicides = total herbicide concentration, PTI = pesticide toxicity index, invert = invertebrates.
Collinearity among factors generally was not strong (Table 3). The strongest bivariate correlations (i.e., |r| > 0.80) were between HOC and flood frequency in 2 MAs and between N and Cl− in 1 MA. Several factors were correlated less strongly (i.e., 0.80 > |r| > 0.70). HOC was correlated with hydrologic factors (2 MAs), Cl− (1 MA), and temperature (1 MA); Cl− was correlated with N (1 MA) and PTI-Invert (1 MA); and hydrologic factors were correlated amongst themselves in 2 MAs.
Table 3.
Pairs of stressor variables correlated (Pearson correlation of log[x]-transformed data) in streams within 6 metropolitan areas. HOC = hydrophobic organic contaminants, PTI = pesticide toxicity index, invert = invertebrate, – = no correlations found.
More variance in biological condition could be explained by the additive effects of 2 variables than by the effects of 1 variable. Most (92 of 118) parsimonious (ΔAICc < 2) models contained 2 predictors (Appendix; available online from: http://dx.doi.org/10.1899/10-131.1.s1). For example, 3 equally plausible models could partly explain variation in macroinvertebrate O/E in Raleigh-Durham. Duration alone explained 47% of the variation, but inclusion of either Cl− or HOCs improved model explanatory ability to >50%.
The amount of variance in biological condition explained by the physicochemical factors varied considerably among assemblages, biological indicators, and MAs (Table 4, Appendix). Overall, plausible models explained 17 to 81% of the variance in biological indicators and were generally poorest in the Dallas-Fort Worth and Milwaukee-Green Bay MAs. For algae and macroinvertebrates, O/E indices were usually less strongly related than NMDS axes to the physicochemical factors (5 of 6 MAs each for algae and macroinvertebrates). In contrast, O/E indices for fish were generally more strongly related to physicochemical factors than to NMDS axes (5 of 6 MAs). Considering either biological indicator, macroinvertebrate assemblages tended to be more strongly related than algal or fish assemblages to the physicochemical factors (5 of 6 MAs).
Table 4.
Mean % variance explained by plausible models (ΔAkaike Information Criterion for small sample size [AICc] < 2; Burnham and Anderson 2002) and total number of plausible models (in parentheses) of physicochemical factors for nonmetric multidimensional scaling scores (NMDS) and observed/expected taxa (O/E) indices of biological assemblages in 6 metropolitan areas. Spearman rank correlations |r| between NMDS and O/E also are given.
The relative importance of physicochemical factors varied among assemblage types and MAs (Fig. 1), but no single factor was universally important. Considering any assemblage type, HOCs and nutrients were most frequently (5 of 6 MAs) the most important factor. Considering any MA, nutrients were the most frequently important factor for algae (5 of 6 MAs), and HOCs were the most frequently important factor for fish (3 of 6 MAs), whereas no single factor was most frequently important for macroinvertebrate assemblages. Considering any assemblage type, the most important physicochemical factors in Atlanta were nutrients and HOCs; in Milwaukee-Green Bay was temperature; in Raleigh-Durham were HOCs and hydrology; in Dallas-Fort Worth were hydrology and nutrients; in Denver-Fort Collins were Cl− and nutrients; and in Portland were temperature, pesticides, and HOCs.
Fig.1.
Summed Akaike weights, which measure the relative importance of each predictor variable, for 1st- and 2nd-order regression models using 11 physicochemical factors to predict biological condition of algae, macroinvertebrate, and fish assemblages as measured by 2 endpoints: the 1st axis from nonmetric multidimensional scaling (NMDS-black), and the ratio of observed/expected taxa (O/E-gray) in 6 metropolitan areas of the USA. HOC = hydrophobic organic contaminants. Pesticides is a collective term for either total herbicide concentration, pesticide toxicity index-Invertebrates, or pesticide toxicity index-Fish as appropriate for each aquatic community. Annotations for each endpoint report the mean % variance explained (from Table 4) for all plausible models.
Discussion
Management practices aimed at protecting or restoring streams in urbanizing watersheds are more likely to be effective when the specific physicochemical factors associated with ecological impairment are understood. If these factors are identified, remedial actions can be targeted, which leads to more cost-effective and ecologically relevant outcomes. We identified physicochemical factors associated with the condition of 3 aquatic assemblages along gradients of urbanization in divergent geographic settings. Our findings suggest that the relative importance of factors depends on the biological assemblage of interest, and that ≥2 factors often can act additively as watersheds urbanize. In addition, we found little evidence that a single type of factor—even those thought to be widely important such as nutrients (USEPA 2006)—was important in all geographic areas.
Correlations among factors limited our ability to rank the relative importance of some factors unambiguously in some MAs, but in general, most factors were not strongly redundant. That we found relatively few strong correlations among factors was surprising because urbanization causes many physical and chemical changes to streams (Paul and Meyer 2001), and many of the factors responsible for these changes are thought to covary (Walsh et al. 2005). One explanation for this result is that different combinations of stressors occur as the intensity of urbanization increases, which is intuitive given the typical sequence of urban development and its attendant effects on stream physicochemical factors. For example, sediment loads in streams may increase as land is cleared and developed, but often decline through time as urban landscapes mature (Finkenbine et al. 2000). In contrast, HOCs typically increase as watersheds urbanize (Van Metre and Mahler 2005).
Our perception of the relative importance of some factors may be an artifact of the study design. We did not characterize all physicochemical factors known to influence streams in urban areas, and some that were measured probably were inadequately characterized. For example, pesticide toxicity was unimportant in many MAs, but samples were collected during baseflow conditions, so exposure probably was underestimated (Sprague and Nowell 2008). In fact, the maximum PTI-invertebrate values reported in our study (Table 2) were lower than values in 90% of urban streams where temporally intensive sampling was conducted across the USA (Gilliom et al. 2006). Reach-scale physical habitat factors were relatively unimportant in most MAs, but this result may be an artifact of the site-selection process, in which an attempt was made to minimize reach-scale differences in habitat conditions (Short et al. 2005).
Our results show the importance of simultaneous consideration of multiple factors when evaluating relationships between biota and physicochemical conditions in urban streams. This need is well rehearsed in conceptual models (Walsh et al. 2005), but appears to be applied rarely (Paul and Meyer 2001). In many studies of stream conditions and watershed land use, investigators use bivariate correlations to search for important factors—a strategy that is fraught with statistical pitfalls (Van Sickle 2003) and may lead to severe underestimation of the importance of factors that act additively or interactively with other anthropogenic or natural factors. No form of model selection, model averaging, or multimodel inference can really determine the relative importance or influence of candidate predictors on a response. Rather these techniques permit investigators to evaluate the relative evidence in the data that variables are predictors of the response. With this caveat, our results show the usefulness of multimodel inference. This method is useful because a single model rarely approximates fully the complex processes of natural systems (Hobbs and Hilborn 2006), especially those under intensive human influence.
Differences in model results and performance between O/E and NMDS indicators can be attributed in part to differences in their derivation. O/E (as computed in our study) measures the deviation in assemblage composition of each site from an expectation derived with a set of regional reference sites independently identified and sampled in other studies (e.g., Carlisle and Hawkins 2008). This property renders O/E (especially if based on presence/absence of common taxa—as for macroinvertebrates in this study) potentially less sensitive than indicators based on relative abundances because we would expect to observe changes in the abundance of sensitive taxa as the severity of anthropogenic factors increase before these taxa became too rare to be detected (i.e., effectively extirpated). This expectation appeared to be true for macroinvertebrates and algae in our study because physicochemical factors often explained less variation in O/E measures than did NMDS axes.
In contrast to O/E, NMDS axes portray the compositional similarity among sites based on relative abundance and species composition along the urban gradient in each MA. This property produces a dimension of change in assemblage composition that is tuned to the sites at hand, but also renders the NMDS axis vulnerable to urban gradients along which changes in assemblage composition are damped, e.g., when sites that are least urbanized are nevertheless disturbed by other anthropogenic activities, such as agriculture. We observed this phenomenon in some of the MAs in our study (Kashuba et al. 2010). When study sites with lowest disturbance along the urban gradient are more disturbed than regional reference conditions, an indicator of biological alteration benchmarked to a larger pool of regional reference sites (i.e., O/E) may be more sensitive to anthropogenic factors than NMDS axes (e.g., as for most fish assemblages in our study). Thus, NMDS and O/E provided unique information (i.e., not redundant) and increased our ability to detect and understand ecological responses to physicochemical factors associated with urbanization. Analysis of more locally derived metrics (i.e., IBIs) that are not available across all MAs would add to understanding at the MA scale.
Our results illustrate the need to consider multiple assemblages in biological assessments. Macroinvertebrate assemblages generally were more strongly associated than other assemblage types with physicochemical factors, but also responded to a different set of factors than the other assemblage types. Thus, macroinvertebrates may be more useful than fish or algae for assessing whether a stream is experiencing some level of ecological alteration caused by urbanization, but provide an incomplete assessment of the physicochemical factors harming the ecosystem. Various aquatic assemblages respond differently to natural and anthropogenic physicochemical gradients because of inherent differences in their environmental requirements and niches (Paavola et al. 2003, 2006, Passey et al. 2004, Johnson et al. 2006, Carlisle et al. 2008).
Acknowledgements
We thank the ecologists, hydrologists, and other specialists on NAWQA's Effects of Urbanization on Stream Ecosystems team, who diligently applied the protocols and, thereby, dramatically improved our ability to summarize this work at a national scale. Review comments by Bob Black and Brian Gregory improved earlier versions of this manuscript.
