Site selection

We used 279 sites from the US Environmental Protection Agency's (USEPA) Environmental Monitoring and Assessment Program – Western Pilot Study (WEMAP). Personnel working on the WEMAP study collected >1500 samples from 1340 perennial streams and rivers in the western US from 2000 to 2004 (Stoddard et al. 2005). Biological, chemical, and physical habitat data were collected from each site following the procedures and guidelines in Peck et al. (2006). Data were collected from most sites during summer (June–September) and from a few during May or October. Site condition was assessed according to local chemistry and habitat conditions. We only used WEMAP sites that were defined by Herlihy et al. (2008) as reference-condition sites. Of the 1340 total WEMAP sites, 326 were considered to meet the criteria for reference condition. We eliminated from this subset sites that lacked macroinvertebrate samples or some habitat variables, as well as a few sites for which derivation of catchment-scale variables was not feasible (e.g., because of ambiguous catchment boundaries). These restrictions resulted in 279 sites available for analysis. For sites having biological, chemical, and habitat data taken on >1 occasion, we randomly selected 1 sample for our analysis.

Data generated for each site

The complete set of local and catchment variables (Appendix; available online from:  http://dx.doi.org/10.1899/10-030.1.s1 (10.1899_10-030.1.s1.doc)) can be characterized as direct climate variables (temperature, precipitation), indirect climate-sensitive variables (hydrologic metrics), and largely climate-independent habitat or landscape variables (land cover, habitat, geomorphology, etc.). Because the sites in our study occur across a broad geographic extent, we also included purely spatial variables (latitude, longitude) in the environmental variable set to examine patterns that might be purely spatially determined.

We derived data at the reach scale from the WEMAP database (taken from 11 transects). We modeled environmental variables at the catchment scale with a variety of sources and techniques after delineating the upstream contributing area to a site from 100-m-resolution digital elevation models (DEMs) obtained from National Elevation Dataset (Gesch 2007). The 45 environmental variables used in the analysis are summarized in the Appendix.

We derived the direct climatic variables, air temperature and precipitation at the catchment scale, from the 800-m-resolution Parameter-elevation Regressions on Independent Slopes Model (PRISM) database (30-y period of record from 1971–2000; PRISM Climate Group, Oregon State University, Corvallis, Oregon;  http://www.prismclimate.org). Percent of total precipitation that falls as snow was calculated according to McCabe and Wolock (2009).

We derived indirect climate-sensitive (hydrologic) variables for the ungauged WEMAP sites using a statistical method (random forests, see Statistical analysis below) as described in Carlisle et al. (2010). For each site, catchment topography and location, land cover, soil properties, geology, and climate variables (precipitation and temperature) were used to estimate hydrologic metrics for the long-term average conditions at the sites (Appendix) based on models calibrated to regional, reference-quality US Geological Survey (USGS) stream gauges (see Carlisle et al. 2010 for details). Most of the hydrology metrics used in our paper are described in Olden and Poff (2003). The mean values of the estimated metrics are calibrated to within 0 to 3% of the observed values for the reference gauges (SD  =  15–40%). Thus, these estimated metrics provide suitable resolution for capturing differences among sites in current flow regimes at ungauged sites. Relative bed stability was calculated according to Kaufmann et al. (2008).

We derived nonclimatic environmental variables at 2 spatial scales. The reach-scale variables were mostly extracted from the WEMAP database. Samples for water chemistry were collected at one point in each site. Reach-scale geomorphology, substrate, riparian and habitat characteristics were evaluated along 11 cross-section transects spanning ∼40× average channel widths between tributaries with methods from Kaufmann et al. (1999). Substrate and wetted-width measurements were taken at an additional 10 sites between transects. We derived the specific stream power between tributaries with a geographic information system (GIS) and calculated it as valley slope × drainage area^0.4 (Cuhaciyan 2006).

We derived catchment-scale nonclimatic variables using a GIS. We derived the proportion of general geology types underlying the catchment from a map developed by Reed and Bush (2005). We derived catchment area, mean watershed elevation and slope, and relief ratio [(mean elevation – minimum elevation)/(maximum elevation – minimum elevation)] from the DEMs. The final catchment-scale geomorphic metric, catchment stream stability, describes catchment slope relative to bed mobility and was calculated as catchment slope × (catchment area/mean substrate size)^0.4. This metric represents a combination of the balance of erosive and resisting forces, dominant channel type, and disturbance regime. We derived the proportion of catchment landuse types from the 2001 National Land Cover Dataset (Homer et al. 2004). Barren landuse type consisted of codes 31 to 33, forest consisted of codes 41 to 43, and agriculture consisted of codes 81 to 85.

Species traits

At each site, the USEPA personnel collected macroinvertebrates by placing one kicknet quadrat sample downstream of each of the 11 cross-section transects used to collect the geomorphology, substrate, riparian, and habitat data (quadrat size was 0.09 m²). The 11 samples were combined and up to 500 (±50) individuals were identified in each composite sample using a fixed count procedure (Stoddard et al. 2005). Macroinvertebrates were identified to lowest feasible taxonomic level, usually genus.

We selected 7 aquatic insect traits from the database in Poff et al. (2006b) to characterize functional composition of the benthic community at each of the 279 sites. These traits reflect a range of attributes, and they are expected to vary in their frequency across the kinds of gradients captured by our environmental variables. Five of these traits (voltinism, thermal tolerance, occurrence in drift, habit, and rheophily) are relatively phylogenetically independent. The other 2 traits (female dispersal and desiccation resistance) are phylogenetically dependent, but they are attributes that should be strongly correlated with hydrologic or other types of disturbance.

Each of these traits has 2 to 5 categorical states to give a total of 21 states. Female dispersal and desiccation resistance have 2 states. Dispersal is categorized as high (>1 km flight before laying eggs) and low (< 1 km), and desiccation resistance is categorized as absent or present. Voltinism, occurrence in drift, thermal tolerance, and rheophily each have 3 states. Voltinism states are categorized as multivoltine (>1 reproductive generation/y), univoltine (1), and semivoltine (<1). Occurrence in drift states are rare (catastrophic drift only), common (typically observed in drift), and abundant (dominant in drift samples). Thermal tolerance states are cold stenothermal, cool/warm eurythermal, and warm eurythermal. Rheophily states are erosional obligate, depositional obligate, and both erosional and depositional. The final trait, habit, has 5 states: burrow, climb, sprawl, cling, and swim. An additional state for the habit trait, skate, is found in the Poff et al. (2006b) trait database, but no taxa with this trait were found at our sites.

We used the approach of Poff and Allan (1995) to transform the taxonomic composition of each site to species-trait composition. For each of the 7 traits, we tallied the proportion of total taxa (usually genera) possessing each trait state to produce a 279 (site) × 21 (trait state) matrix.

Statistical analysis

We conducted 4 sets of statistical analyses. First, we used a hierarchical partitioning of variance approach to determine how variation in trait composition across the 279 sites was explained by the environmental variables at the local vs catchment scales and in terms of climatic, hydrologic, and nonclimatic variables (see Appendix). We modified Cushman and McGarigal's (2002) hierarchical canonical correlation analysis (CCA) to account for underlying linear trait–environment relationships that are better captured with redundancy analysis (RDA; Heino et al. 2007). In the hierarchical structure, we initially divided the environmental variables into catchment- and local-scale variables with latitude and longitude as purely spatial variables. In the next level of the hierarchy, we further divided local-scale variables into chemistry, riparian, and stream-structure variables. We divided stream-structure variables into geomorphology and substrate/habitat variables. We divided catchment-scale variables into climatic (temperature, precipitation) plus climate-sensitive variables (hydrology) and nonclimatic variables (geology, land use, geomorphology). Two variables, proportion of total precipitation that falls as snow and mean relative humidity of catchment, represented interactions between precipitation and temperature, but we included them in the precipitation category to maintain just 3 climatic-sensitive categories (temperature, precipitation, hydrology). We used the vegan package in R to run the RDA analyses (Oksanen et al. 2009).

Second, we used cluster analysis to characterize similarity among sites in terms of trait composition. Using complete linkage, we identified functionally similar groups of sites in the 7-dimensional trait space (as in Poff and Allan 1995). We used the complete linkage method because it detects naturally formed groups (Hill and Lewicki 2007) and tends not to artificially produce evenly sized clusters, a desirable feature for our dataset. We used the elbow method to select the number of clusters (Ketchen and Shook 1996) by assessing the additional amount of variation explained by increasing the number of clusters. This technique revealed that 3 and 8 clusters optimized the amount of variation explained for our data. We selected the 3-group solution based on ease of interpretation of the groups. Each group of sites had communities with similar trait composition that may have different responses to environmental change, including climate change. We used the R program hclust (in the stats package) to create the clusters (version 2.10.0; R Development Core Team, Vienna, Austria).

Third, we used classification and regression tree (CART) analysis (De'ath and Fabricius 2000) to examine whether spatial variation in the 3 similar traits-based community types could be explained in terms of the multiscaled environmental variables. CART models (Breiman et al. 1984, Clark and Pregibon 1992) have the desirable features of accommodating nonlinearities and correlations between explanatory variables. A CART analysis accounts for the variation in the response variable by splitting sites into increasingly homogenous groups (nodes) until criteria are met for stopping. To avoid the problem of over-fitting by generating too many terminal nodes (Hill and Lewicki 2007), we used the cross-validation process described in De'ath and Fabricius (2000). We used the tree package in R to run the CART analysis (Ripley 2009).

CART produces only a single, static best-fit tree. Therefore, correlated environmental variables that produce similar levels of homogeneity in a split will not be expressed in the model, and variables important for defining community types could be missed. To overcome this deficit, we used a random forest procedure (Breiman 2001) to determine which environmental variables best explained the trait clusters. Random forest algorithm uses a bootstrapping technique to minimize over-fitting and determines variable importance by running thousands of CART models, each produced with a randomly selected subset of explanatory variables and 70% of the data (see Carlisle et al. 2010). The remaining 30% of the data are used to evaluate the model. Throughout this process, the random forest algorithm measures the amount of error reduced by inclusion of each variable as well as the relationships between variables to produce a measure of importance for each variable. We ran the random forest analysis based on 5000 trees with 18 variables for each split in the randomForest package in R (Liaw and Wiener 2002).

Fourth, we explored sensitivity of sites to projected climate change based on 2 particularly sensitive trait states. Cold stenothermal taxa are expected to be directly sensitive to warming of surface waters because they prefer cool waters. Taxa with a strong preference for flowing water (erosional obligate taxa) could be expected to respond indirectly to reduced precipitation, which would diminish stream flow and, thus, the availability or permanence of flowing-water microhabitats. We used these 2 trait states to characterize the sensitivity of benthic communities at each of the 279 sites to future warming and precipitation reduction in 2 ways: proportion of cold stenothermal taxa present, and proportion of erosional obligate taxa present. For each of these classifications, we used CART analysis to explore how much among-site variation in these sensitivity traits was explained by the environmental variables followed by a random forests procedure to assess overall environmental variable contribution.

We used these results to evaluate geographic distribution of benthic community vulnerability to climate change by overlaying projected future changes in temperature and stream runoff on the geographic distribution of site sensitivities. We obtained the projected change in runoff in the conterminous US from fig. 4.10 in Lettenmaier et al. (2008), which was derived from values in Milly et al. (2005). The figure in Lettenmaier et al. (2008) shows model-averaged change in median runoff as the average runoff in 2041–2060 relative to average runoff in 1901–1970 for each of the USGS hydrologic unit regions (HUC2) in the western US. We recalculated these values from Milly et al. (2005). We overlaid the proportion of erosional obligate taxa at each site on these changes in runoff to examine which sites were most likely to be sensitive to projected runoff changes.

We overlaid the proportion of cold-stenothermal taxa at each site on the change in average temperate in 2041–2060 relative to 2051–2070 to examine which sites are most likely to be sensitive to projected changes in temperature. We obtained the projected temperature data from the World Climate Research Programme's (WCRP) Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel data set. The CMIP3 data set consists of 16 distinct models, each model having from 1 to 7 different sets of initial conditions (runs). In addition, each model was calibrated for 3 distinct CO₂ scenarios. These bias-corrected and spatially downscaled climate projections derived from CMIP3 data are available from:  http://gdo-dcp.ucllnl.org/downscaled_cmip3_projections/. For our analysis, we selected the middle scenario where CO₂ emissions increase until about mid-21^st century. Because no single global circulation model (GCM) has valid a priori superiority, we followed the procedure of model-averaging of the 16 ensemble runs (cf. Milly et al. 2005) to obtain a singular temperature change for each pixel. In addition, we averaged the change in temperature within each HUC2 to make our results compatible with the change in runoff. All maps in our paper were produced using ArcMap™ 9.2 GIS software (ESRI, Redlands, California).

General relationships between trait composition and environmental variables

Hierarchical RDA revealed that 67% of the total variation in functional (trait) composition across the 279 sites was explained by the environmental variables used in the analysis. Explained variation was partitioned into local (reach) vs catchment and climate vs nonclimate variables (Fig. 1). When we accounted for the other 2 groups of variables in the 1^st level of the hierarchy, catchment variables explained 12.2% of the variation compared to reach (4.9%) and spatial (0.4%). Much variation was explained by >1 set of variables, especially catchment–local (20.4%) and catchment–local–spatial (19.5%). This latter interaction reflected the dominant contribution of streamwater conductivity, which was related to variables in all 3 groups in sites in the northern Great Plains.

Figure 1

Results of hierarchical redundancy analysis showing partitioning of explained variation in functional composition of 279 sites according to 45 environmental variables (see Appendix). Total variation explained (67%) was partitioned into spatial, local (reach), and catchment variables and their interactions. All variables were partitioned into reach-scale vs catchment-scale variables. Reach- and catchment-scaled variables explained 45% and 61% of total variation, respectively. Reach-scale variables were further partitioned into component variables (riparian, chemistry, and habitat structure [which consisted of component variables]), and catchment-scale variables were further partitioned into climate and nonclimate variables (each with their component variables). Numbers in parentheses indicate the number of component variables in each partition.

Partitioning of catchment variables indicated that climate variables explained substantially more variation (55.6%) than nonclimate variables (41.8%). Most of that variation (36.4%) could be explained as an interaction between both sets of variables. However, 19% of the variation was explained by climate variables and their interactions. After accounting for interactions, the climatic variables temperature and precipitation explained only <1% and 2.3%, respectively, and hydrologic variables explained 7.8%. The main contribution of these environmental variables was in their interaction, which accounted for 32.9% of total variation. This result was not surprising, given the strong relationships among climate-driven variables. Among the nonclimate variables, geomorphic metrics explained the most variation (15.2%), although the interaction with land use was large (14.4%).

Partitioning of local variables showed that structural habitat metrics accounted for 38.8% of total variation in trait composition, but most of that variation (26.4% of total variation) also could be explained by either chemistry or riparian characteristics. Substrate/habitat and geomorphic metrics contributed about equally to overall explained variation.

Species trait groups and environmental correlates for western WEMAP sites

Three groups of sites with similar trait composition were identified using cluster analysis. The trait-state composition of sites in these 3 groups for each of the 7 traits is shown in Fig. 2.

Figure 2

The average proportion of taxa in each trait state for the 3 community types defined by cluster analysis. The cold, stable (CS) type was present at 152 sites with higher proportions of taxa with trait states characteristic of colder and more hydrologically stable streams. The warm, unstable (WU) type was found at 33 sites with higher proportions of taxa with trait states characteristic of warmer and more unstable streams. The mixture of trait states (M) type was found at 94 sites with communities with a wide variety of trait states.

The 1^st community type was the largest and was found at 152 sites. It was characterized by the highest average proportions of taxa that were cold stenothermals (  =  0.45), univoltine (  =  0.72), erosional obligates (  =  0.24), rare in the drift (  =  0.45), and poor female dispersers (  =  0.69). This community type was dominated by clingers (  =  0.47) with a mix of other habits. These communities were generally characteristic of relatively cold mountain streams where taxa had traits expected in cool habitats (cold stenothermy, univoltine life cycle) and stable habitats (low mobility, adult dispersal, rheophily). We designated this community type as cool, stable (CS).

A 2^nd community type was found at 33 sites. It was characterized by the highest average proportions of taxa that were multivoltine (  =  0.46), cool/warm eurythermal (  =  0.77), depositional obligates (  =  0.44), abundant in the drift (  =  0.65), and strong female dispersers (  =  0.69). This community type had a high proportion of burrowers (  =  0.52) and sprawlers (  =  0.25) and was associated with warmer, lowland streams. We designated this community type as warm, unstable (WU).

A 3^rd community type was found at 94 sites. It was an intermediate community type where most traits consisted of mixtures of trait states. For example, the proportions of all 3 states of the occurrence in drift trait were substantial: rare (  =  0.35), common (  =  0.20), and abundant (  =  0.45). The functional flexibility of communities in this group of sites suggests that these sites would be found across the western US at localities ranging from mountain to lowland streams. We designated this community type as having a mixture of trait states (M).

CART analysis

We eliminated 2 sets of variables to evaluate response of the 3 community types to climate change. Spatial variables explained very little variation (see Fig. 1) and were removed. Streamwater conductivity had a very strong geographic signature (sites in upper Great Plains) and tended to be a dominant splitter in the classification tree. Because this variable was correlated with many variables associated with this region (e.g., base flow: r  =  −0.64, July temperature: r  =  0.56), we eliminated it from the analysis in favor of climatic variables that allowed us to evaluate better how a change in climate may modify these sites. Elimination of conductivity did not qualitatively change the structure or interpretation of the classification tree.

The CART analysis of the 3 traits-based community types included 2 climate-influenced hydrologic variables (Fig. 3). Mean annual flow was the initial split and base flow further split sites with low annual flow and higher stream slope. Catchment-scale geology, catchment geomorphology, and local habitat and geomorphology further divided the community types. The misclassification rate for the model was 15%. The tree placed CS community type at sites that fell under 2 sets of environmental conditions. The great majority of CS sites were associated with high mean annual flow and a high proportion of fast-water habitat (terminal node 8 in Fig. 3). A smaller group of 19 CS sites were associated with low annual flow, high basin relief, high base flow, and high reach slope (terminal node 6). A group of 14 M sites also was associated with these environmental conditions.

Figure 3

The classification tree for the 3 traits-based community types (CS  =  cold, stable, M  =  mixture of trait states, WU  =  warm, unstable; see Fig. 2). The explanatory variables displayed in the regression tree are mean annual runoff (C.H.MAF), mean catchment slope (C.M.Slp), proportion of fast-water habitat within the reach (R.H.Fast), baseflow index (C.H.BF), proportion of crystalline geologic type underlying catchment (C.G.Cry), mean slope of reach (R.M.Slp), and relative bed stability (R.S.RBS). The histogram below each terminal node shows the distribution of community types within that node.

WU sites were all associated with low mean annual flow, but with 2 sets of environmental conditions. Terminal node 1, with 23 WU sites, was defined by low mean annual flow and low basin relief. Another group of 8 WU sites in terminal node 3 was associated with low annual flow, high catchment slope, a high proportion of crystalline geology, and low relative bed stability.

The classification tree placed M community type under environmental conditions that would support various trait states, such as high mean annual flow and low proportion of fast-water habitat. The M sites were associated with a range of environmental conditions, including both high and low mean annual flow, base flow, and proportion of crystalline geology.

The random forest analysis of the community types supported the classification tree analysis by including most of the same variables, but with more climatic variables. Mean annual flow was the most important variable, followed closely by proportion of total precipitation that falls as snow, February precipitation, catchment slope, fast-water habitat, base flow, July temperature, median annual coefficient of variation of daily flows, and annual precipitation. Proportion of total precipitation that falls as snow was correlated with both proportion of fast-water habitat and base flow. February precipitation and annual precipitation are highly correlated with each other as well as with mean annual flow. These correlations would indicate that purely climatic variables are important in defining the environmental differences among the CS, WU, and M sites shown in classification model.

The geographic distribution of community types revealed distinct patterns (Fig. 4). The WU community type was distributed among arid lowland sites, with most sites in the northern Great Plains but many in the desert southwest lowlands. The CS community type occurred primarily in the Rocky Mountains, mountains of the Great Basin, and coastal mountain regions. The M community type was scattered across the region, but most sites were in mountain and foothill streams, especially the more arid mountains of southern California and Arizona.

Figure 4

Geographic distribution of the 3 community types predicted by cluster analysis (see Fig. 2 for abbreviations). The underlying ecoregions are Omernik (1987) Level III ecoregions, which have been generalized into 3 regions, as in Stoddard et al. (2005).

Sensitivity of benthic communities to climate change

CART analysis of the 2 sensitive trait states (cold stenothermal and erosional obligate) revealed the strong influence of climatic variables (Fig. 5A, B). Sites with a higher proportion of cold stenotherms (>0.30) had cold July temperatures or, if they had higher July temperatures, had lower daily flow variation with fewer low-pulse counts (Fig. 5A). The random forest analysis of cold stenotherms also revealed July temperature as the most important variable, followed by proportion of total precipitation that falls as snow. These 2 climatic variables were followed by 4 hydrology variables: high-flow pulse duration (indicative of snowmelt), mean annual flow, coefficient of variation of daily flows, and base flow.

Figure 5

Regression trees for the proportions of 2 sensitive trait states. A.—Cold stenotherms, with explanatory variables (see Appendix) mean July air temperature of catchment (C.T.JulT), median annual coefficient of variation of daily flows (C.H.CV), reach turbidity (R.C.Tur), and mean low-flow pulse count (C.H.LFPC). B.—Obligate rheophiles with explanatory variables mean February precipitation for the catchment (C.P.FebP), median annual coefficient of variation of daily flows (C.H.CV), and proportion of fast-water habitat within the reach (R.H.Fast). The value at each terminal node is the average proportion of sites grouped in that node. Percent variance explained for each tree was 69.0% for cold stenotherms and 40.3% for erosional obligates.

Sites with a higher proportion of erosional obligates (>0.25) occurred either in sites with low February precipitation regimes and a high proportion of fast-water habitat or high February precipitation and lower daily flow variation (Fig. 5B). The random forest analysis of erosional obligates also revealed February precipitation as the most important variable, followed by proportion of sand and fines, and proportion of fast-water habitat. These dominant variables were followed by mean annual flow, catchment slope, annual precipitation, and coefficient of variation of daily flows.

We overlaid the proportion of cold stenotherms on the map of projected change in temperature from 1951–1970 to 2041–2060 (Fig. 6A) and the proportion of erosional obligates over change in runoff from 1901–1970 to 2041–2060 (Fig. 6B). Changes in runoff and temperature are projected to be most severe in the upper Colorado River and Great Basin HUC2s. Sites with the greatest proportion of cold stenotherms (0.50–0.70; Fig. 6A) are presumably most sensitive to warming and tended to be in the CS community type, which were primarily in mountainous regions. Of the 50 thermally sensitive sites (proportion of thermally tolerant taxa between 0.50 and 0.70, see Fig. 6A), 44 (88%) were CS sites and 6 were M sites. Thus, only 29% of the 152 CS sites were considered highly thermally sensitive. Similarly, 7 of the 11 sites sensitive to streamflow reduction (i.e., proportion of taxa that are erosional obligates between 0.35 and 0.50, Fig. 6B) were CS community types, and the rest were M types. Only 1 site (in the Pacific Northwest HUC2) was defined as both thermally sensitive and sensitive to streamflow reduction.

Figure 6

Proportion of taxa at each site that were cold stenotherms (A) or obligate erosional taxa (B) overlaid on site maps showing projected changes from 1951–1970 to 2041–2060 in temperature (A) and 1901–1970 to 2041–2060 in runoff (B). Changes in temperature and runoff were averaged across US Geological Survey hydrologic unit regions (see text for discussion). The traits-based community types were cold, stable (CS), warm, unstable (WU), and mixture of trait states (M) (Fig. 2). The proportion of taxa with the sensitive trait in the community is represented by size of the symbol.

Traits-based benthic communities and environmental gradients

Our analysis is the first attempt to relate species trait composition of benthic macroinvertebrate assemblages in undisturbed (reference) streams to multiscaled environmental drivers with the goal of distinguishing the relative contributions of climatic and nonclimatic environmental correlates of trait composition as a basis for examining sensitivity to projected climate change. As expected, we found that climatic variables contribute to a substantial portion of among-site variation in trait composition at the broad geographic extent of the western US. This result may reflect, in part, the greater number of climatic and hydrologic variables than nonclimatic and local variables (see Appendix). However, many nonclimatic variables at both reach and catchment scales also were important in explaining community variation.

Previous studies of environment–macroinvertebrate trait relationships using constrained ordination (Feld and Hering 2007, Heino et al. 2007) found that local-scale and landuse characteristics accounted for much of the variation in trait composition, but they did not incorporate catchment-scale climate or climate-sensitive (hydrologic) variables. This omission may explain why their analyses accounted for less geographic variation in trait structure than ours did. In addition to the climate-based variables, habitat (substrate and geomorphic variables) explained substantial portions of variation in our data set. However, most of the explained variation in trait composition at all sites could be accounted for by >1 set of variables, a result indicating a high degree of variable interaction and a complex response across natural gradients of environmental variation for these reference sites.

We found reasonably strong separation among sites grouped by trait composition along dominant environmental gradients (Fig. 3). In particular, hydrologic metrics and both stream and basin slope were primary environmental factors discriminating the groups. Despite their relative weight in the hierarchical RDA analysis, climatic variables did not appear in the CART as primary splitting variables in identifying traits-based community types, but they were very heavily weighted in generating multiple possible trees. The climatic variables, proportion of total precipitation that falls as snow, February precipitation, and annual precipitation, were the top 3 variables, and these had strong correlations (r²  =  0.62–0.97) with the hydrologic variables, mean annual flow and base flow, which were identified by the CART analysis as primary splitting variables. This result indicates that climatic factors are important in distinguishing among community types, which is not surprising given the large subregional (ecoregional) differences in annual precipitation driven by prevailing climatic patterns and steep elevation gradients in the western US. Together, these variables influence volume and seasonality of runoff and extent and duration of late summer base flow. This hydrological influence presumably allows the community types to be associated with different climatic settings. For example, most sites that supported CS assemblage type (terminal node 8 in Fig. 3) were high elevation, snowmelt-driven sites with high velocities characteristic of steep mountain streams throughout the region. By contrast, WU sites, which occurred at lower elevations in flatter terrain, possessed mostly opposite trait states and were distinguished by inverse correlations with these same climatic and landscape variables. Similar distinctions between CS and WU community types also were seen along elevation gradients in the Colorado River Plateau, an area not represented in our study (Moline 2007).

Catchment-scale factors were important descriptors of benthic assemblages at a regional scale, a result that has been observed in other studies. For example, using functional feeding groups and habit traits, Heino et al. (2007) found that catchment, local, and spatial variables accounted for 43.8% of total variation in trait-based community structure of macroinvertebrates in Finland. Catchment-scale variables accounted for 22% of total variation, local-scale variables accounted for 27% of total variation, and spatial variables accounted for 22% of total variation. They found that ∼½ of the explained variation (20%) was involved in interactions between the 3 sets of variables. Heino et al. (2007) found that spatial variables were much more important in explaining variation than we did. However, they did not include climatic variables per se in their analysis, and it may be that latitude is correlated with climate conditions that influence environmental conditions in Finnish streams. In the western US, environmental conditions of streams are more influenced by steep elevation gradients, which are not necessarily correlated with latitude or longitude, except for the easternmost plains sites included in the data set.

The differences in trait composition of our 3 community types were substantial, especially for the CS and WU types (Fig. 2). Interestingly, the suite of traits that define the CS and WU types tend to be associated across the North American lotic insect taxa (see table 3 in Poff et al. 2006b), a result indicating that these trait groups can be expected to co-occur often as trait complexes in stream benthic communities. The catchment- and local-scale environmental correlates of these trait groups clearly differ, a result suggesting that a single traits-based definition of reference condition might not be appropriate for the entire western US.

In contrast, researchers studying variation in trait composition at large, regional scales in Europe have found only small differences in trait composition among reference-quality streams. Archaimbault et al. (2005) found fairly small but statistically significant differences in some traits between streams with different underlying geology at 12 sites in France and Belgium. Bonada et al. (2007a) found small differences in traits-based community structure between temperate and Mediterranean streams in Europe and North Africa, and Johnson et al. (2004) found few differences in functional feeding groups across 6 ecoregions in Sweden. Statzner and Bêche (2010) argued that these small differences are statistically but not ecologically significant and that the small observed ecological distinctions in trait structure imply a universal trait-based community structure for reference condition streams within major climatic zones (e.g., the Temperate Zone). Statzner and Bêche (2010) reinforced this point by demonstrating highly correlated trait compositions between 527 reference-condition sites across Europe and 435 reference-condition sites across the conterminous US.

Our contrasting findings may reflect the fact that our study occurred in a region with steep climatic and nonclimatic environmental gradients and that we were able to use a set of derived and measured metrics to capture this environmental heterogeneity effectively. Archaimbault et al. (2005) included both lowland and mountainous sites with various geology types, but all sites were located within a single climatic region. Bonada et al. (2007a) and Johnson et al. (2004) generated average trait composition within biogeographic regions or ecoregions but in doing so averaged over mesoscale factors (sensu Bonada et al. 2007a), such as the effects of mountains and lowlands within an ecoregion. Similarly, the intercontinental-scale analysis of Statzner and Bêche (2010) incorporated sites in the US that were mainly in humid or perennial flow regions (the Pacific Northwest, Rocky Mountain, or Mid-Atlantic Highlands), and their study may have been statistically insensitive to trait variation associated with the few sites located in the xeric West or plains, where harsher hydrologic conditions and warmer climatic conditions exist.

Our study sites also spanned several ecoregions (Fig. 4), but we incorporated climatic, climate-influenced, and nonclimatic variables for every site and, thus, expressed trait composition as average values for an entire ecoregion. This coupling of at-a-site trait composition with multiscaled and multitype environmental data probably allowed us to detect arguably different trait reference conditions.

The question of how to define traits-based reference condition(s) for assessing benthic community response to environmental change requires more attention because the answer may depend on how the question is addressed. Not only is the issue of capturing landscape (mesoscale) variation in trait structure important, but the sensitivity of reference type to the set of traits selected also is unresolved. We used a subset of traits from Poff et al. (2006b) to define trait groups, and these traits differed qualitatively from those used by Bonada et al. (2007a) and Bêche and Statzner (2009). Others have used still different subsets of traits to infer the need for multiple reference trait groups (e.g., functional feeding groups along pH gradients, Heino 2005). We may have obtained different results with a different set of traits, although we have not examined that possibility. Thus, our results are suggestive of the need to consider more than a single traits-based reference condition within minimally human-altered streams. They are not conclusive, but they do point to a need for further analysis.

Risk and vulnerability to climate change

The 3 assemblage types (Fig. 2) varied substantially in trait composition in ways that reflect direct and indirect sensitivity to climatic and other nonclimatic environmental features. For example, number of generations per year is a thermally sensitive trait for many species (Vieira et al. 2006), and thermal tolerance is defined in terms of temperature. However, other traits are difficult to relate directly to temperature change, and these traits may be more easily interpreted in terms of disturbance resistance (e.g., juvenile mobility in the drift and adult female dispersal). Therefore, we used only the cold stenothermal and obligate erosional (rheophilic) trait states to examine sensitivity to climate change.

Our vulnerability analysis for western WEMAP reference sites showed clear subregional differences. The Upper Colorado and Great Basin HUC2s are projected to have the most dramatic increases in temperature and reductions in runoff for the western US (see also Seager et al. 2007). Therefore, sites with high proportions of cold stenotherms in the Upper Colorado and Great Basin HUC2s may experience relatively large shifts in community composition caused by selection against this trait state. This region is where sensitive CS community types (as well as CS and M community types with a moderate proportion of thermally sensitive taxa) are most likely to experience the greatest stress from increasing temperatures. By contrast, thermally sensitive CS types in other areas, such as the Pacific Northwest, are expected to experience less exposure to thermal warming and so should be relatively less vulnerable to changing traits-based community structure driven by a warming climate.

A similar geographically distributed vulnerability is associated with projected declines in runoff. Communities with high proportions of obligate erosional taxa are presumably at risk because decreased runoff can cause intermittency or prolonged periods of low flow (and low velocities). The highest proportions of obligate erosional taxa occur in CS community types. Therefore, CS community types in the Upper Colorado, Lower Colorado, and Great Basin HUC2s are projected to be most vulnerable to community alteration. Other regions where communities are runoff sensitive are not projected to experience reductions in annual runoff (e.g., the Pacific Northwest), and benthic assemblages in these regions presumably would not experience selection against runoff-sensitive taxa. However, other components of runoff might be sensitive to change, even if annual flows remain similar. For example, thermally driven reductions in snowpack can modify the timing and duration of high and low flows and potentially could affect benthic invertebrates via indirect pathways mediated by habitat. A pattern of earlier snowmelt in the western US has already been documented in the late 20^th century (Stewart et al. 2005).

WU and possibly some M-type communities appear most robust to increased temperature and decreased runoff. The WU group in particular is composed of taxa having traits characteristic of warm and intermittent-flow conditions. These community types probably have already been filtered by current (mesoscale) climatic and hydrologic factors for taxa that are able to tolerate environmentally harsh conditions (see Poff and Ward 1990, Bonada et al. 2007b, Hering et al. 2009). Taxa at these sites might be more resilient to additional warming, although additional information on lethal tolerances clearly is needed. At the regional scale, we might expect to see taxa characteristic of these groups becoming more widespread and dominant across the region as climate change progresses.

Implications of climate-driven changes for benthic community composition

As stream temperatures and runoff patterns change to some substantial degree, we would expect to see species replacements and changes in the traits-based community composition of sites. Given the general correspondence between thermally sensitive sites and the 3 trait-based community types, we would expect to see a shift in communities from CS toward M and WU types. However, linkage among species traits combined with slow adjustment of nonclimatic environmental conditions could complicate this pattern. For example, high-gradient (mostly CS) sites with high local stream velocities would not be expected to provide the kinds of microhabitats available in lowland WU sites. Thermal conditions might shift in these more montane sites and favor cool/warm eurythermal taxa. However, such taxa are largely depositional obligates and burrowers, and they might not be able to colonize effectively and thrive in warming but still erosional sites. About 62% of taxa at CS sites in the Upper Colorado HUC2 are cold stenotherms, and loss of these taxa because of a strong climate-change filter would greatly alter community composition at these sites. Consider the following hypothetical illustration. Suppose that all cold stenotherms were removed, and warm/cool eurythermal taxa randomly colonized these sites. If we assume a 50% probability of successful colonization for those potential colonizing taxa that are also depositional obligates or burrowers, then, on average, a 12% reduction in genus richness would occur at these sites (from 36 to 32). Obviously, many scenarios of species loss and colonization success are possible, and more focused modeling efforts will be needed to examine this question rigorously.

Species losses and gains in response to warming could also affect components of ecosystem function, such as energy flow through detritus- and algal-based food chains. For example, the 297 sites in our study contain 24 taxa that are defined as shredders, and 16 of these are cold stenotherms. Were warming temperatures to eliminate cold stenotherms, then we would expect to see a loss of total shredder richness across many sites. We might expect to see the dominant shredder genera shift from cold stenotherms like Pteronarcys and Zapada to cool/warm eurytherms like Tipula and Malenka, which might have different production and litter-processing efficiencies. A similar kind of species substitution might be seen among grazers. Of the 51 herbivore genera across the 297 sites, 25 are cold stenotherms, 24 are cool/warm eurytherms, and only 2 are warm eurytherms. Cold stenotherms, such as Heptagenia or Neothremma, could be replaced by cool/warm eurytherms, such as Drunella or Glossosoma. Grazer control of algal production in high-elevation western streams depends on grazer species and their performance under near-bed flow (Wellnitz and Poff 2001, Poff et al. 2003), so shifting grazer species composition or modifying runoff could alter the environmental context of grazer–algal interactions and stream energy flow.

However, we wish to emphasize that, based on our analysis, the response of local communities to climate change probably will depend on nonclimatic environmental factors. For example, we observed a high proportion of thermally sensitive taxa in sites with warmer air temperatures but relatively stable flows (low coefficient of variation and few low-flow pulses:  =  0.315 in Fig. 5A). This combination of environmental factors suggests the possibility that local flow regime (and temperature regime) may be influenced by groundwater inputs, a variable that cannot be reliably modeled. Indeed, previous work by Hawkins et al. (1997) has shown that local controls on stream temperature can obscure expected differences based on air temperature (e.g., elevation) in montane streams in the western US. Such local controls provide an important context (and challenge) to predicting stream benthic responses to climate change in heterogeneous landscapes (see also Hering et al. 2009, Hawkins et al. 2010).

This study allowed us to begin exploring how specific reference sites may be affected by climate change, as a function of both site-specific sensitivity and subregional degree of climate change. The WU reference sites in the northern Great Plains appear to be intrinsically robust against increased temperature or drying, but this region is not expected to experience much warming or drying, so little biotic change would be expected in any case. However, the M and CS sites in the higher-elevation ecoregions of the western US are more likely to respond to climate change because these sites contain more taxa sensitive to warming (cold stenotherms) and drying (erosional obligates). CS sites are generally the most vulnerable, but they have different vulnerabilities depending on geographic location relative to projected late-21^st-century warming and drying. Sites in the humid, cooler Pacific Northwest have relatively low vulnerability compared to CS sites in the Upper Colorado and Great Basin region. However, sites with only moderate proportions of thermally sensitive or flow-sensitive taxa may still be vulnerable in regions with large projected climate warming or drying. Thus, moderately sensitive sites in regions with large climate excursions might have a level of vulnerability similar to that of highly sensitive sites in other regions with less warming or drying. These kinds of considerations must be evaluated more carefully with future research to understand better the relative vulnerabilities of WEMAP reference sites and to project the future taxonomic and trait composition of these sites. These sites will remain the least human-impacted sites in the western US and will continue to provide an important future reference for disentangling the effects of climate change from other forms of anthropogenic disturbance.