Study area

Southern California has a Mediterranean climate, with hot, dry summers and cool, wet winters. Intermittent streams are typical of the region. Much of the coastal, lower-elevation areas have been converted to agricultural or urban land uses, and water importation, runoff, and effluent discharges have perennialized most streams (Roy et al. 2009, Mazor et al. 2012). However, most of the upper elevations of watersheds remain undeveloped, and have chaparral, grassland, and oak or pine forest land covers. Streams in these undeveloped portions are mostly nonperennial, although short perennial lengths can be found near surficial bedrock or spring sources. Much of the region is underlain by young, erodible sedimentary geology.

Sampling

We worked in 12 nonperennial and 3 perennial streams in southern California. We defined nonperennial streams as streams that lack surface flow for at least several days per year in most years. This definition encompasses a wide variety of streams, from ephemeral washes and headwaters that flow for only a few hours after storms, to those with sustained flows lasting nearly all year or that have perennial flow in a year with heavy rainfall. All nonperennial streams in our study had flow from groundwater discharge that persisted ≥1 mo after the last storm. By the definition of Yavercovski et al. (2004), these streams would be seasonal or near-permanent.

We selected the 12 nonperennial sites to represent a range of natural conditions (e.g., short and long flow duration, high and low gradient) and anthropogenic stressors (e.g., urban development, grazing, water diversion).We sampled 3 sites repeatedly over several years, and we deployed data loggers to measure water level and temperature continuously for 1 y at selected sites (too few loggers were available for deployment at every site or every year of sampling).

We sampled nonperennial sites during the drying phase, typically beginning in March or April and ending when flow was insufficient for sample collection (no surface water in <50% of the sampling reach). Sampling began in the eurheic phase, continued through the oligorheic (or rarely the arheic) phases, and ended when the stream was in a drier hydrologic state (hyporheic or edaphic phases) (terms follow Gallart et al. 2012). We revisited sites approximately every 2 to 4 wk.When possible, we increased the frequency of site visits to weekly near the end of the drying phase. We sampled perennial sites monthly beginning in April to overlap the normal index period for bioassessment in perennial streams in southern California (SMC 2007). Sampling of nonperennial streams was funded mostly under a program that began in 2008, was suspended in 2009 for budgetary reasons, and resumed in 2010. Sampling of perennial streams was funded mostly under a 2^nd program that began in 2009 and continued in 2010.

Overall, the study period was drier than the long-term average for the region. Normal annual rainfall for Lindbergh Airport in San Diego is 27.4 cm. In 2008, 2009, and 2010, annual rainfall was 18.3, 23.4, and 26.9 cm, respectively. Drought was even more severe in the 2 y preceding our study, with only 13.7 and 9.9 cm of rainfall in 2006 and 2007, respectively ( http://www.sdcwa.org/annual-rainfall-lindbergh-field).

Stressor assessment

We used best professional judgment to select nonperennial sites that represented a gradient of stress. We verified or modified a priori selections after evaluating each site for 50 stressors (hydrologic, physical, land use, biological) associated with the California Rapid Assessment Method (CRAM) (CWMW 2012), a measure of riverine wetland condition. A stressor was given a score of 0 if it was not observed, a score of 0.5 if the stressor was likely to have a negative effect on the stream, and a score of 1 if the effect was likely to be large. We examined the distribution of scores to identify 3 groups: low stress (scores ≤ 1), moderate stress (scores 2–7), and high stress (scores 7.5–9; higher scores were not observed). A maximum score of 50 is theoretically possible, but scores >6 are uncommon, except at highly developed sites. The highest score observed in a statewide probabilistic data set of 924 sites was 21, and only 25% of the sites in this statewide data set would be classified as high stress (data not shown). Our approach to establishing a stressor gradient by counting stressors is similar to the approach used by Sánchez-Montoya et al. (2009a) to identify reference sites. Details on stressors identified at each site were provided by Mazor et al. (2012) and are given in Appendix S1.

Assessment of hydrologic regimes

We evaluated hydrologic regimes and stressors with a combination of data sources: continuous data loggers, direct measurements, and visual observation during site visits. We deployed continuous water-level loggers (HOBO U20 Water Level Data Logger; Onset, Bourne, Massachusetts) at a subset of sites on the 1^st sampling date and retrieved them after the last sampling date. At 4 sites (Bear Canyon [BC], South Fork Santa Ana River [AN], Cedar Creek [CD], and Noble Canyon [NC]), we corrected waterlevel measurements for air pressure measured with a 2^nd logger deployed at the site above the water line. At the other sites, we used air-pressure data from nearby weather stations. We calculated discharge from water velocity measured with an electromagnetic or propeller-type velocity meter, but we used flotation time of a neutrally buoyant object when conditions were too slow or shallow for the velocity meter. Last, we noted whether streams were flowing on each sampling date. We used these data to identify periods when the reaches contained flowing water, and when they were dry (or intermittently dry). We gave visual observation the highest priority, followed by direct measurements. Long periods (>6 h) during which loggers recorded water-level readings <0 were interpreted as times when the stream was dry.

Benthic macroinvertebrate collection

We used standard state bioassessment protocols to sample benthic macroinvertebrates (Ode 2007). California protocols are derived from those developed by the US Environmental Protection Agency for national stream surveys (Peck et al. 2006). We divided a 150-m reach into 11 equidistant transects. At each transect, we used a 500-µm D-frame kick net to sample 929 cm² of stream bed at 25, 50, or 75% of the transect distance. At 3 low-gradient sites where stable habitats were restricted to the stream banks (Aqua Caliente Creek [AC], San Juan mainstem [SJ], and Pine Valley Creek [PC]), we sampled at 0, 50, and 100% of the transect distance (Mazor et al. 2010). We combined all samples into a single composite sample and preserved it in 70% ethanol. For each sample, we removed ∼600 invertebrates from the detritus and identified them to Standard Taxonomic Effort (STE) Level 2 (mostly to species with Chironomidae to genus, where possible) established by the Southwest Association of Freshwater Invertebrate Taxonomists (Richards and Rogers 2011). We sampled macroinvertebrates only when flow was sufficient for sampling, i.e., the stream was in a eurheic or oligorheic state for ≥50% of the reach.

Habitat characterization

We measured physical-habitat variables by standard state protocols (Ode 2007) with a few modifications. First, we fixed the locations of the 11 transects so that they did not vary over the course of the study. Thus, we sampled well beyond the protocol-mandated time at which a site would be rejected because of lack of wet habitat (Ode 2007). Second, we assumed that slope, gradient, bank width, and bank height were stable, and we measured these variables only once a year. Third, we added measurements of algal cover in 2009 and 2010 after standard methods of estimation were published (Fetscher et al. 2009). Where possible, we analyzed physical-habitat data by calculatingmetrics published by Kaufmann et al. (1999). We analyzed other data (e.g., data from field probes) without modification.

Landscape variables

We calculated % impervious surface and road density for the watershed area upstream of each site. We used National Landcover Data from 2006 ( http://www.epa.gov/mrlc/nlcd-2006.html) to assess imperviousness, and we derived road density from a custom geographic information system (GIS) road layer (PRO, California Department of Fish and Wildlife, unpublished data).

Data analysis

Benthic macroinvertebrate community structure We used the IBI (Ode et al. 2005) and nonmetric multidimensional scaling (NMDS) to summarize macroinvertebrate community structure. We graphed these measurements over time and compared them with environmental variables to examine changes in benthic communities over the course of the study.

We calculated metrics for the IBI as described by Ode et al. (2005). To maintain consistent sample sizes, we used random subsampling to reduce the number of individuals in each sample to 500, and then aggregated individuals from STE Level 2 to STE Level 1 (i.e., mostly to genus, with Chironomidae to family; Richards and Rogers 2011). We calculated 7 metrics: 1) Coleoptera richness; 2) Ephemeroptera, Plecoptera, and Trichoptera (EPT) richness; 3) predator richness; 4) % collectors; 5) % intolerant individuals; 6) % noninsect taxa; and 7) % tolerant taxa. These metrics are scored so that higher scores reflect lessdegraded condition. We summed metric scores and rescaled the sum to a 100-point scale, as described by Ode et al. (2005). We compared values to a threshold of 39, which is 2 standard deviations (SD) below the reference calibration mean, as described by Ode et al. (2005).

We ran NMDS on presence/absence data from all samples to examine differences in assemblage composition (PC-ORD, version 5.12; McCune and Mefford 2006). We used Bray—Curtis distance, and tested up to 4 axes via 100 Monte Carlo runs with real and randomized data. We used the default procedure in PC-ORD to select axes, i.e., the highest dimension that reduced stress by ≥5 (on a scale of 0–100) and had less stress than that of 95% of runs with randomized data. We set the maximum number of iterations to 250, the stability criterion to 0.000001, and step length to 0.2. We rotated final axes with varimax rotation, and calculated taxon scores with weighted averaging. We checked scores for association with pollution tolerance by calculating correlations with tolerance values (from the California Aquatic Macroinvertebrate Laboratory Network; CAMLnet 2003).

Table 1.

Names and characteristics of sites sampled in our study.

We used a linear, mixed-effects model to test whether flow and stress status affected IBI scores. We used the nlme package in R (version 3.1–109; Pinheiro et al. 2013) to create a mixed-effects model for each biological response variable based on 4 classes of flow and stress (perennial low stress, nonperennial low stress, nonperennial moderate stress, and nonperennial high stress) with site specified as a random effect. Differences in variability between low-stress perennial and nonperennial sites were examined with an F-test using pooled variance estimates.

Relationships between biological and environmental variables We used 2 separate analyses to examine the relationship between environmental variables and assemblage composition (i.e., IBI scores and NMDS axis scores). First, we used Spearman rank correlations between within-site means of environmental and biological variables to evaluate gradients associated with among-site differences. Second, we evaluated the gradients associated with within-site changes by calculating Spearman rank correlations on variables after subtracting the within-site mean. We used only those environmental variables expected to vary with changing hydrological conditions for the 2nd analysis, and we excluded variables expected to be constant over the course of the study (e.g., all landscape variables, bank dimensions). Our goal was to characterize the relationships between biological and environmental variables, so we did not test the statistical significance of these correlations (McCune and Grace 2002).

Stressor assessment

Perennial (AN, BC, and CD) and 3 of the 12 nonperennial sites (Temescal Creek [TE], AC, and Carney Canyon [CC]) were under low stress (Table 1). Six nonperennial sites were under moderate stress (NC, Arroyo Seco [AS], San Diego River headwaters [SR], Cañada Verde [CV], PC, and Santa Ysabel Creek [SY]), and 3 were under high stress (Ortega Falls [OF], SJ, Trabuco Creek [TC]). The most stressed site (TC) had a score of 8.5. Stressors included grazing (SY, SR, and CV), nutrients (SJ and OF), and urban runoff (PC and TC) (Table 1, Appendix S1; Mazor et al. 2012). We found evidence of nonnatural flow regimes (e.g., floods in the absence of precipitation at AS) at several sites. Communication with nearby land managers indicated that groundwater diversions affect sites SJ and OF.

Patterns in hydrologic regimes

Data-loggers showed that hydrologic patterns differed among sites, including the 3 perennial sites (Fig. 1). Hydrographs were most stable at perennial sites, particularly AN, and at certain nonperennial sites (i.e., SY, SR, and NC), where water levels were stable until they decreased abruptly. In contrast, SJ and OF showed periods of large fluctuations. Communication with the manufacturer of the loggers (Onset) indicated that these fluctuations may indicate periods of extremely low water levels, during which the logger may be partially exposed to air. During these periods, intermittent drying and rewetting of the logger are likely. Intermittent drying also was evident at AS, and this site experienced several short high-flow events, despite the lack of recent precipitation. We did not deploy loggers at low-stress nonperennial sites.

Direct observation at 3 sites with multiple years of sampling showed that hydrologic characteristics also varied among years (Fig. 2A–F). All 3 sites dried in June or July during the 1^st sampling year (2008 for SY and SR, and 2009 for NC). However, in 2010 (the wettest year of the study), only SR dried in June. SY started to dry in September, and flow never ceased at NC. Data from multiyear deployment of water-level loggers were not available.

Figure 1.

Water surface levels at a subset of sites. Shading indicates dry or fluctuating periods. Data for SJ, SY, SR, AS, and OF were obtained in 2008; BC, AN, and NC were obtained in 2009; and CD were obtained in 2010. Where water levels are >0 and background color is gray (as at sites SJ and OF), surface water is inferred to be present during a period of water-level fluctuation.

Trends and patterns in IBI scores

IBI scores were high at low-stress sites, regardless of flow status (Fig. 3A, B, Table 2) or discharge (Fig. 4). The 3 low-stress nonperennial sites had IBI scores comparable to scores at low-stress perennial sites, and no sample from these sites had a score <39 (reference threshold). Flow status had a small (6.6 ± 5.4 [SE] points) nonsignificant effect on IBI score (linear mixed model, p = 0.65; Fig. S1), but power was low (only 3 sites in each low-stress group). IBI scores responded to stress at nonperennial sites as expected (Fig. 3C, D, Table 2). Mean scores declined as site quality declined (r = -0.91). Low-stress sites had higher mean IBI scores (55.6) than moderate-stress sites (40.3), which had higher scores than high-stress sites (26.8). Nearly all samples from high-stress sites had scores <<39.

IBI scores were least temporally variable at low-stress sites (Fig. 3A, B). At low-stress sites, flow status did not affect variability (pooled SD = 3.3 at perennial sites and 4.0 at low-stress nonperennial sites, F_11,11 = 1.47, p = 0.27) and no directional trend was evident at perennial or nonperennial low-stress sites. In contrast, IBI scores declined sharply near the end of the sampling period (typically in June) at several stressed nonperennial sites (OF, AS, SJ) and underwent erratic fluctuations at others (CV, NC). Discharge and IBI score were not strongly related (Fig. 4). At 2 sites that were resampled over 2 y (SR and SY), a mid-spring decline in IBI scores was followed by an early summer increase. This pattern was similar in both years of sampling, but no individual metrics appeared to consistently influence this pattern. The 3^rd temporally replicated site (NC) did not show similar patterns in each year, perhaps because it underwent different hydrologic regimes each year (see below).

Figure 2.

Photos of sites SR (A, B), SY (C, D), and NC (E, F) sampled over multiple years. Each photo was taken in the same month of each year at approximately the same location in each reach.

Stability of IBI scores arose from patterns that varied across sites (Fig. 5A–D). In some cases, stability in the IBI score arose from stability in the underlying metrics (e.g., AN; Fig. 5A). In other cases, decreases in some metric scores were offset by increases in others (e.g., CD; Fig. 5B). However, large fluctuations in metric scores that yielded smaller fluctuations in the IBI score were more common (e.g., TE; Fig. 5C). At sites where the IBI declined (e.g., OF; Fig. 5D), metric scores declined in unison near the end of the sampling period.

Figure 3.

Trends in Index of Biotic Integrity (IBI) scores at perennial low-stress (A), nonperennial low-stress (B), nonperennial moderate-stress (C), and nonperennial high-stress (D) sites. Each trajectory represents samples from a single site in a single year of sampling. See Fig. S1 for trends for individual sites. The dotted line indicates the threshold IBI score for reference condition.

Trends and patterns in community structure

Perennial and nonperennial streams did not have distinctly different macroinvertebrate community types. Instead, they represented extremes along a continuous biological gradient (Fig. 6A, B). NMDS resulted in a 3-axis solution that represented >80% of the total variance with a final stress of 17.8. Axis 3 represented more variance (41%) than the other axes (20 and 19% for axes 1 and 2, respectively). Visual inspection showed weak segregation of sites by flow status along axes 2 and 3. This segregation was driven substantially by 1 site (BC on axis 2 and AN on axis 3), but samples from the other 2 perennial sites were dispersed among the nonperennial sites. Samples also segregated weakly along a stressor gradient, with low-stress sites clustered at the positive ends of axes 2 and 3, and the negative end of axis 1.

Taxa segregated primarily along axis 3. EPT taxa were more common at sites with high positive values on axis 3 (Fig. 6D). Coleoptera occupied a slightly lower position along axis 3, suggesting a shift from EPT taxa to beetles as sites moved down this axis. In contrast, major taxonomic groups were strongly interspersed along axes 1 and 2 (Fig. 6C). Tolerance values were more strongly correlated with species scores along axis 3 (Spearman's ρ = -0.49) than axis 1 (ρ = 0.36) or axis 2 (ρ = -0.9). Taxa did not segregate by functional feeding group along any axis.

Table 2.

Site stress scores and mean and standard deviation (SD) Index of Biotic Integrity (IBI) scores. n = number of samples, n > 39 = number of samples with scores greater than the threshold for identifying reference condition. See Table 1 for definitions of site codes.

We plotted trajectories in ordination scores as change in score from initial sampling and found that community composition at most sites was either stable or changed in unison (Fig. 7A, B). For example, all trajectories were either constant or moved in a positive direction along axes 1 and 2, or in a negative direction along axis 3. Thus, if sites changed at all, they changed toward communities more characteristic of nonperennial and high-stress sites, with fewer EPT taxa and more Coleoptera, Diptera, and noninsects. In general, trajectories were shortest for perennial and low-stress nonperennial sites, and longest for high-stress nonperennial sites.

Figure 4.

Index of Biotic Integrity (IBI) scores vs measured discharge. For clarity, only a subset of sites are shown. Trajectories connect samples within sites in order of sampling date, and move in the direction of increasing discharge (i.e., toward the left end of the x-axis). The dotted line indicates the threshold IBI score for reference condition.

Relationships between biological and environmental variables

Many environmental variables were associated with biological differences among sites, but few were associated with differences within sites (Fig. 8A–D, Table 3). For example, 11 of the 28 environmental variables evaluated had strong relationships (|ρ| > 0.5) with IBI scores, and 16 variables had strong relationships with ≥1 of the 3 ordination axes. The strongest relationships of environmental variables with IBI scores were observed for stressor scores (ρ = -0.91), a few habitat variables (e.g., % fastwater habitat: ρ = 0.75; % shading: ρ = 0.68), and specific conductance (ρ = -0.69). For ordination axes, stressor scores were less important (i.e., strongest ρ = -0.53, with axis 2), but many habitat variables, particularly those related to substrate (e.g., % cobble embeddedness: ρ = -0.85 with axis 3), hydrology (e.g., % fast-water habitat: ρ = 0.70 with axis 2), and riparian vegetation (e.g., mean riparian vegetation cover: ρ = 0.70 with axis 3) showed strong relationships with site differences. Specific conductance also was strongly related to ordination axes (e.g., ρ = 0.65 with axis 3).

Figure 5.

Trends in metric scores at AN (A), CD (B), TE (C), and OF (D). Gray lines represent trends in scores for each of the 7 metrics in the Index of Biotic Integrity (IBI). The black line represents the mean metric score.

In contrast, few environmental variables were associated with within-site biological changes, and only 2 were associated with IBI scores with |ρ| > 0.2 (temperature: -0.31, wetted width: 0.27). Relationships with ordination axes were stronger, but only 4 of the 22 variables considered had |ρ| > 0.5. The strongest relationships were for variables related to water quality (e.g., temperature: ρ = 0.66 with axis 2; dissolved O₂: ρ = 0.57 with axis 3) and hydrology (e.g., discharge: ρ = -0.57 with axis 2, wetted width: ρ = -0.67 with axis 2). Physical-habitat variables not directly related to water availability had much weaker relationships with within-site biological changes. Among these variables, the strongest relationships were observed for % sands and fines (ρ = -0.32 with axis 1) and % shading (ρ = 0.43 with axis 2). The strength, and sometimes the direction, of these relationships varied from site to site. For example, within-site correlations between IBI scores and temperature were negative for most sites, but ρ was positive at TE (ρ = 0.1) and strongly positive at NC (ρ = 1.0).

Assessment tools work in nonperennial streams

Nonperennial streams are widely described as supporting biological communities that are different from communities in perennial streams (Álvarez and Pardo 2007, Datry 2012, Bogan et al. 2013). However, fauna from the streams in our study appear to occupy different positions along a continuous gradient. For example, samples from perennial and nonperennial streams were interspersed in ordination space, emphasizing the overall similarity of these stream types. The relative similarity of perennial and nonperennial streams observed in our study contrasts with results of several other studies (e.g., Arscott et al. 2010). However, many of these studies had limited spatial (e.g., Bêche et al. 2006) or temporal (e.g., Lunde et al. 2013) replication, and these limits alone could exaggerate the apparent distinctness of nonperennial streams. Our findings are supported by those of Gallart et al. (2012), who suggested that ecological status may be assessed in the same way at temporary streams as at perennial streams if flow permanence and seasonal predictability are relatively high.

Figure 6.

Results of a nonmetric multidimensional scaling (NMDS) ordination of macroinvertebrate presence-absence data showing ordination scores for each sample in the study (A, B) and weighted-average taxon scores (C, D) on axes 1 and 2 (A, C) and 1 and 3 (B, D). Centroids for 2 sites referenced in the text (AN and BC) are annotated. Numbers in axis labels show the % total variance represented in the axis. P-L: perennial, low-stress streams. N-L: nonperennial, low-stress streams. N-M: nonperennial, moderate-stress streams. N-H: nonperennial, high-stress streams.

The comparability of bioassessments between perennial and nonperennial streams may be limited to regions like southern California, where nonperennial streams are the dominant stream type. Benthic macroinvertebrates collected at a site represent a subset of the regional fauna that can tolerate the local environmental conditions (Southwood 1977, Townsend and Hildrew 1994, Statzner et al. 1997). In predominantly arid regions, taxa that are dependent on perennial streams may be excluded from the regional fauna because of the scarcity of habitat (e.g., only 500 km of streams in the San Diego region are perennial) (NHDSPlus;  www.horizon-systems.com/nhdplus). Thus, communities in perennial streams in arid regions may be constrained to be similar to communities in nonperennial streams. In addition, temporal variability in flow status, as observed in our study and elsewhere (e.g., Gasith and Resh 1999), may further limit the regional fauna to taxa adapted to intermittent flow. The richer regional fauna of wetter regions may allow greater divergence between perennial and nonperennial streams. The findings of Arscott et al. (2010) support this idea. They observed that communities in intermittent and ephemeral streams in New Zealand were nested subsets of communities found in nearby perennial streams, and reflected a loss of desiccation-sensitive taxa, rather than a gain of desiccation-tolerant taxa.

Figure 7.

Trajectories in ordination space for each site and year plotted against time for nonperennial (A) and perennial (B) streams. In each panel, scores for the initial sampling event were subtracted from scores for all subsequent samples at a site and year. Abbreviations are the same as in Fig. 6.

Figure 8.

Among-site correlations between environmental variables and nonmetric multidimensional scaling (NMDS) ordination axes calculated from within-site means (A, B) and variables with within-site means subtracted (C, D). The length of each line represents the strength of correlation (ρ) with the ordination axes. Only selected variables are labeled. Fast = % fast-water habitat, Slow = % slow-water habitat, Shade = % shading, CPOM = % cover by coarse particulate organic matter, Vel = mean water velocity, flow = discharge, macro = % macroalgae cover, embed = mean embeddedness, SC = specific conductance, T = water temperature, stress = stressor score, RipVeg = mean riparian vegetation cover, roads = road density, dry = % dry habitat, big = % particles larger than cobble, SAFN = % sands and fines, cobble = % cobbles, width = mean wetted width, DO = dissolved O₂.

The similarities between communities in perennial and nonperennial streams in our study may be restricted to the types of nonperennial streams we studied (long-lasting, seasonal streams that flow for several months nearly every year). Less-predictable or more-ephemeral streamsmay support biological communities that differ more from those in perennial streams than what we observed. Anna et al. (2008) found large differences in bioassessment indices between ephemeral and intermittent streams, and Gallart et al. (2012) found that unpredictable hydrologic regimes decreased the measured ecological status of temporary streams. Therefore, bioassessment tools developed from perennial streams may provide comparable interpretations only at nonperennial streams with sufficient relative flow permanence. Identification of critical thresholds in flow permanence or predictability should be a focus of further research.

Table 3.

Correlations of environmental and biological variables with the Index of Biotic Integrity (IBI) and nonmetric multidimensional scaling (NMDS) ordination axes 1, 2, and 3. Among sites = Spearman rank correlation (ρ) calculated from site means, within sites = Spearman rank correlation (ρ) calculated from variable values minus within-site means, CPOM = coarse particulate organic matter, NA = not assessed.

Assessment scores are consistent over time

Multivariate analyses showed large seasonal changes in community composition at most sites, but IBI scores were relatively stable, at least at low-stress sites. This stability was particularly evident at perennial sites, where all samples ranged ∼5 points at a single site. The consistency was less pronounced at low-stress nonperennial sites (typical point ranges ≈ 10 points). However, at moderate- and high-stress nonperennial sites, IBI scores were noticeably unstable and declined steeply at several sites, particularly at the end of sampling. Some point ranges at these sites were >40 on the 100-point IBI scale. In unstressed intermittent streams in Spain, Sánchez-Montoya et al. (2009b) similarly found low intra-annual variability in multimetric indices and in certain metric types, such as those derived from life-history traits, further suggesting that benthic communities in undisturbed nonperennial streams are relatively resilient to seasonal changes.

Stability in IBI scores at nonperennial sites also extended to interannual replicates, at least at 3 of the 4 sites selected for multiyear sampling, a finding consistent with the observations of others (Morais et al. 2004, Sánchez-Montoya et al. 2009b). Despite different weather patterns over the years of study, these 2 sites had similar patterns in IBI scores and trajectories in ordination space. The stability in IBI scores at low-stress sites may be a result of the use of a broad index period (April–October) during development of the IBI (Ode et al. 2005). Thus, the IBI appears to integrate the temporal variability evident in the ordinations and creates a stable measurement of structural and functional attributes of the benthic macroinvertebrate assemblage. Our results suggest that a common assumption underlying development of many bioassessment indices—the substitution of space for time during index development—may effectively broaden their applicability (Stoddard et al. 2008, Schoolmaster et al. 2012). That is, indices calibrated from data covering many sites, but only a few years, may be resilient to temporal changes. When a large number of reference sites is used, streams are likely to be sampled at different stages in their individual phenologies, which may produce biological variability similar to what would be found in repeated visits and, thus, a biotic index with resilience to temporal variability.

Few changes in the environment were associated with trends in biology

Given the overall stability of IBI scores, it is not surprising that few environmental variables were related to withinsite changes in scores. Moreover, only a few variables were associated with changes in ordination scores. Unsurprisingly, these variables were mostly related to water quantity or quality (e.g., velocity, % fast-water habitat, temperature, dissolved O2); relationships with variables related to physical structure (e.g., substrate composition), primary productivity, and riparian vegetation were weaker. In studies of temporary streams in other Mediterranean-climate regions, investigators have found similar relationships among these variables and macroinvertebrate community structure (e.g., Morais et al. 2004, Álvarez and Pardo 2007). García-Roger et al. (2011) found that most changes in habitat were site-specific, but those changes related to quantity of aquatic habitat were strongly associated with changes in community structure. In contrast, among-site differences were stronger and associated with a larger variety of environmental variables, consistent with many studies of temporary (e.g., Álvarez and Pardo 2007) and perennial (e.g., Sandin and Johnson 2004, Yuan et al. 2008, Mazor et al. 2009) stream ecosystems.

One of our more surprising findings was the apparently high sensitivity of nonperennial streams to nonnatural flow regimes. For example, several sites that had very few disturbances apart from altered hydrology (AS, NC) had relatively low IBI scores, and the lowest scores were observed at sites with altered hydrology combined with other stressors (TC, OF, and SJ). In contrast, sites with substantially higher stress but relatively stable hydrographs (SY, SR) had much higher IBI scores. Skoulikidis et al. (2011) also found that biological assemblages were particularly sensitive to modified flow regimes and artificial drying, and recommended distinguishing naturally from artificially nonperennial streams for assessment purposes. Therefore, watershed managers should monitor alterations to hydrologic regimes in nonperennial streams as much as they do in perennial streams.

Implications for bioassessment programs

Our study was limited to a small number of sites, but it illustrates that certain nonperennial streams can be incorporated into routine bioassessment programs. Hydrologically stable and predictable nonperennial streams sampled during eurheic or oligorheic states are biologically similar to perennial streams, at least in arid regions like southern California. Expanding water-quality assessment programs to include nonperennial streams would give resource managers the ability to manage a greater extent of their streams and to address impacts to some of the most sensitive portions of their watersheds. These changes may be most profound in arid regions, like southern California, but the global ubiquity of nonperennial streams (Tooth 2000) suggests that watershed protection in both wet and dry climates could be greatly improved by including nonperennial streams in assessment programs.