Forecasting responses of benthic community structure and function to anthropogenic climate change is an emerging scientific challenge. Characterizing benthic species by biological attributes (traits) that are responsive to temperature and streamflow conditions can support a mechanistic approach for assessing the potential ecological responses to climate change. However, nonclimatic environmental factors also structure benthic communities and may mitigate transient climatic conditions, and these must be considered in evaluating potential impacts of climate change. Here we used macroinvertebrate and environmental data for 279 reference-quality sites spanning 12 states in the western US. For each sampling location, we described 45 environmental variables that spanned reach to catchment scales and that represented contemporary climate drivers, hydrologic metrics, and nonclimatic habitat features, as well as purely spatial metrics. We described benthic community composition at each site in terms of 7 species traits, including those considered sensitive to temperature increases and streamflow changes. All combined environmental variables explained 67% of the total trait variation across the sites, and catchment-scale climatic and hydrologic variables independently accounted for 19%. Sites were clustered into 3 community types based on trait composition, and a classification-tree analysis confirmed that climatic and hydrologic variables were important in partitioning these groups. Sensitivity of benthic communities to projected climate change was assessed by quantifying the proportion of taxa at sites having the traits of either cold stenothermy or obligate rheophily. Regression-tree analysis showed that temperature and hydrologic variables mostly accounted for the differences in proportion of sensitivity traits across the sites. We examined the vulnerability of sites to climate change by superimposing regional-scale projections of late-21st-century temperature and runoff change on the spatial distribution of temperature- and runoff-sensitive assemblages. Sites with high proportions of cold stenotherms and obligate rheophiles occur throughout the western US, but the degree of temperature and runoff change is projected to be greatest for reference sites in the Upper Colorado River and Great Basin. Thus, our results suggest that traits-based sensitivity coupled with intraregional variation in projected changes in temperature and runoff will cause reference sites in the western US to be differentially vulnerable to future climate change.
Forecasting responses of benthic community structure and function to anthropogenic climate change is an area of growing scientific and management interest. A challenge in doing so is to distinguish biological responses to climate change from a variety of other environmental factors that influence local species success and community membership. The composition of benthic invertebrate communities in a locality is affected by the interaction of many environmental variables acting at multiple spatial and temporal scales in both direct and indirect fashion (e.g., Poff 1997, Allan 2004, Heino et al. 2007). Therefore, to detect the effects of climate change on benthic species and communities, we need to first understand how benthic communities are currently regulated by climatic factors and then how climate interacts with other nonclimatic environmental determinants of benthic community structure. Because human activities have extensively modified natural background environmental variation (Allan 2004, Wohl 2005, Poff et al. 2006a, Humphries and Winemiller 2009), initial efforts to understand potential responses to climate change are most reliably done in least-disturbed or minimally altered reference sites (Stoddard et al. 2006).
Two key aspects of climate that strongly influence stream ecosystems are temperature and precipitation. Water temperature and the annual thermal regime directly regulate species' growth and set limits on species distributions across the landscape (Sweeney and Vannote 1978, Ward and Stanford 1982). Projected warming of the earth's atmosphere will modify thermal regimes in ways that have profound implications for stream ecosystems and include shifting the composition of macroinvertebrate communities (Sweeney et al. 1992, Meyer et al. 1999, Statzner and Bêche 2010).
Precipitation also exerts a strong control on stream structure and function by being translated via watershed processes into runoff in stream channels. Geographic variation in precipitation regimes (magnitude, seasonality, etc.), in combination with geologic and land-cover controls on rates and pathways of runoff, translate into geographic variation in streamflow regimes, which directly and indirectly regulate many aspects of benthic species performance and community composition (Poff and Ward 1989). Therefore, changes in precipitation regimes are likely to induce shifts in macroinvertebrate community structure and perhaps function (Poff et al. 2002).
Biological response metrics that are sensitive to temperature and streamflow are needed to evaluate how climatic factors influence the distribution and abundance of benthic species and shape community composition. The sensitivity of these response metrics must be evaluated relative to nonclimatic environmental factors (e.g., geomorphic setting, riparian condition) to provide a habitat context for projecting benthic responses to rapid climate change. By understanding how contemporary patterns of benthic species distribution and community structure are explained by prevailing direct climatic factors (temperature, precipitation), indirect climate-sensitive variables (flow-regime metrics), and nonclimatic habitat conditions, we can potentially gain insight for predicting future responses to changing climate in heterogeneous landscapes.
The use of species traits offers a promising approach for relating contemporary patterns of benthic species distribution and community structure to prevailing environmental factors, both climatic and nonclimatic. A traits-based approach presumes some direct or mechanistic linkage between environmental drivers and biological responses that facilitates causal inference (Poff 1997, Statzner and Bêche 2010, Webb et al. 2010). Traits have been used widely in benthic stream ecology (e.g., Townsend and Hildrew 1994, see review in Statzner and Bêche 2010). Uses include efforts to link trait composition of benthic communities to a variety of prevailing environmental drivers across local to catchment to geographic scales (e.g., Richards et al. 1997, Townsend et al. 1997, Finn and Poff 2005, Heino 2005, Dolédec et al. 2006) and to use trait composition of benthic communities as metrics of biological condition in a biomonitoring context (Statzner et al. 2005, Culp et al. 2010, Menezes et al. 2010, Pollard and Yuan 2010).
In principle, a traits-based approach can be applied to questions of climate change because changes in temperature and streamflow act as ecological filters to select for biological traits suitable to the new environment (Poff 1997). A few traits-based studies have incorporated climate variables, such as temperature (Snook and Milner 2002, Lamouroux et al. 2004, Statzner et al. 2005, Bêche and Statzner 2009); however, only Bonada et al. (2007a) have attempted to use extensive trait information to infer potential benthic responses to climate change. No previous effort has attempted to determine how benthic responses to climate change can be separated from responses to nonclimatic sources of natural environmental variation.
An important consideration in evaluating potential traits-based responses to climate change is how to define the reference condition. Because the traits approach potentially allows benthic communities of substantially different taxonomic composition to be compared directly across broad geographic scales (Townsend and Hildrew 1994, Poff 1997), a traits-based reference condition may have a correspondingly broader geographic extent than a taxonomic-based reference condition (Statzner et al. 2005). Indeed, there is some suggestion that within broad climatic zones (e.g., the Temperate Zone), a single traits-based reference condition may suffice for all minimally human-altered sites (Statzner and Bêche 2010). However, whether traits-based reference conditions are stable across strong environmental gradients (e.g., ecoregions) within major climatic zones has not been adequately examined.
We had several objectives in our study. First, we quantified the independent and interactive influence of multiple environmental variables on benthic community trait composition in the context of current climatic regimes. In doing so, we examined whether distinctive traits-based reference conditions occur and whether these can be associated with gradients of prevailing climate. Second, we evaluated the potential sensitivity of benthic communities in terms of key traits that might reasonably be expected to respond directly to increasing temperature and modified precipitation or stream flow. We wanted to understand better how such community-level sensitivity to climate change may be mediated by local or landscape context, particularly by nonclimatic environmental factors. Third, we explored the geographic distribution of benthic community vulnerability to projected changes in temperature and runoff across the western US.
We analyzed variation in traits-based community structure of benthic macroinvertebrates in reference-quality sites across the western US. These sites are presumed to have minimal human alteration, a feature that allowed us to focus on natural gradients of environmental variation. We used environmental information on climate (temperature and precipitation), catchment land cover, hydrology, geomorphology, riparian condition, channel slope, and sediment at each sampling location. We characterized traits composition of benthic communities with several species traits deemed to be sensitive or responsive to the environmental drivers. We explained variation in trait composition in terms of multiscaled environmental drivers, including those directly linked to climate, indirectly influenced by climatic variables, and essentially independent of climate. We then assessed whether projected climate change is likely to be a detectable and significant stressor for benthic communities having differing sensitivity and occurring in different habitat and climatic settings.
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 area0.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.
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 m2). 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.
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 CO2 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 CO2 emissions increase until about mid-21st 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 1st 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.
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.
The 1st 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 2nd 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 3rd 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).
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.
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.
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.
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.
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 (r2 = 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 20th 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-21st-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.
We acknowledge the modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI), and the WCRP's Working Group on Coupled Modelling (WGCM) for their roles in making available the WCRP CMIP3 multimodel data set. Support of this data set is provided by the Office of Science, US Department of Energy. We thank A. T. Herlihy for providing the WEMAP data, Dan Baker for assisting in development of catchment variables, and P. R. Kaufmann for deriving his relative bed-stability metric for our sites. This work was supported in part by USEPA Science to Achieve Results (STAR) grants SPO BS0056363 and R833833, but it does not reflect any official endorsement by the US Environmental Protection Agency.