Study area

The climate of KNP (Fig. 1) is characterized by mean temperatures that range from −6°C in July to 21°C in January. The region experiences 4 distinct seasons. Average annual rainfall recorded at the Perisher Valley gauging station within KNP is 1828 mm. Most streams in the Australian Alps are chemically similar to streams in European and North American alpine regions, with high concentrations of dissolved O₂ and low concentrations of nutrients and dissolved salts. They often have torrential flows during snow melt (Tiller 1993), but they carry low sediment loads compared to many northern-hemisphere high-altitude streams (Good 1992). Above 1500 m, frosts and blizzards can occur at any time, and snow usually accumulates between June and September. The snow line generally extends down to 1200 m. Almost continual melting of snow under the snow pack occurs during winter, and the runoff can lead to substantial transport of materials to streams (Molles and Gosz 1980). However, no snow remains in the study area throughout the year.

Figure 1

Locations of test and reference study sites in Kosciuszko National Park, Australia.

Developments, such as the Snowy Mountains Hydro-electric Scheme, have affected the hydrologic, geomorphologic, and ecological condition of many park streams (DECC 2006). Ski resorts are another major type of development within the park that can damage streams via discharge of treated sewage effluent, development activities, parking lots and associated runoff of road-deicer and salt-laden sediment, and abstraction of water for domestic supply and artificial snow making (Cullen et al. 1992, Digance and Norris 2006). In January 2003, bushfires burned 71% of the park's 673,492 ha, including the areas surrounding many of the study sites (DECC 2008).

Sixteen sites (11 test and 5 reference) were sampled on 5 streams. Test sites were positioned strategically to assess ecological responses to resort activities (Fig. 1). The sites range in elevation from 1160 m asl on Sawpit Creek (site 162) to 1760 m asl on Spencer's Creek (105). Extensive snow cover often prevents sampling at these high-altitude sites in winter. The most upstream site on all streams (11, 105, 121, 128, and 160) are minimally disturbed sites. These 5 sites are a subset of the many reference sites originally sampled to collect data to build the KNP predictive model used in our study. We resampled these 5 reference sites to provide a baseline from which to identify any regional shift from reference condition as indicated by the macroinvertebrate assemblage.

Macroinvertebrate sampling and bioindicators

Macroinvertebrates were collected quarterly in February (summer), May (autumn), August (winter), and November (spring) (quarters 1–4, respectively) each year from autumn 1994 to autumn 2008. Samples were collected from riffle habitats at each site using standard AUSRIVAS sampling techniques (Parsons and Norris 1996, Nichols et al. 2006). Each collection was subsampled to collect ∼200 animals/sample (Marchant 1989). All sites (test and reference) sampled during the study period were assessed by the Thredbo predictive model, which was developed specifically for biological assessment of the streams in the KNP region based on AUSRIVAS modeling approaches (Simpson and Norris 2000), which are similar to those described for RIVPACS (Wright 1995, Wright et al. 1998, Moss et al. 1999). The AUSRIVAS macroinvertebrate collections and analysis provided the observed/expected (O/E) taxa scores. To aid interpretation of model outputs, AUSRIVAS O/E taxa scores are allocated to different bands that represent levels of biological condition. Bandwidths are based on the distribution of O/E taxa scores for the reference sites in any particular model. Band A is centered on 1.0 and includes the middle 80% of the reference-site O/E taxa scores. For the Thredbo model, band A is between O/E taxa scores 0.87 and 1.14. Condition of test sites with an O/E taxa score within Band A is considered equivalent to reference condition. A test site with an O/E taxa score below the lower bound (<0.87 for the Thredbo model) is judged to be significantly impaired with fewer families than expected and is allocated to one of the lower bands according to its score (Coysh et al. 2000).

We calculated 4 other bioindicators for each macroinvertebrate collection: the SIGNAL score for the sample, which is scaled from 1 (presence of many tolerant taxa) to 10 (presence of many sensitive animals) (Chessman 2003), the ratio of EPT families to all families in the sample, Simpson's Diversity, and the number of Oligochaeta in the sample (log₁₀[x]-transformed). We selected these bioindicators because they represent a variety of commonly used bioindicators, including those based on richness, composition, abundance, and sensitivity.

Bureau of Meteorology data

The Australian Bureau of Meteorology supplied us with climate data up to 2008 for weather stations at Thredbo and Perisher ski resorts. Records at Perisher began in June 1976 and Thredbo in January 1969. We selected climate variables for their potential effects on flow and temperature, which in turn, can influence temporal variation (interannual and seasonal) in macroinvertebrate assemblages (Wade et al. 1989, Linke et al. 1999, Mazor et al. 2009). We used highest monthly temperature (°C), lowest monthly temperature (°C), mean 3-PM relative humidity (%), and number of days with rain (full descriptions of variables are available from:  http://www.bom.gov.au/climate/cdo/about/definitions9and3.shtml). We converted all climate data used in the analysis to quarterly series with the discontinuous piecewise constant curve method, which fits the average value (of the monthly data) to each quarterly interval and applies the most recent value to missing data points (SAS, version 9.1.3; SAS Institute, Cary, North Carolina). The meteorology variables from Perisher and Thredbo were highly correlated (all correlation coefficients >0.9). Therefore, we used the average value for the 2 weather stations when relating climate data to macroinvertebrate data. Only data from 1994 onward were included in the analysis to maintain correspondence between the 2 data sets.

Statistical analysis

We expected seasonal differences in climate data and macroinvertebrate assemblages because sites were sampled quarterly (Bêche et al. 2006). Preliminary analysis confirmed that macroinvertebrate time series and climate variables both showed a significant seasonal effect. Therefore, we ran all macroinvertebrate analyses separately for each season. We seasonally adjusted the time-series data for each climate variable with the US Bureau of the Census X-11 Seasonal Adjustment program by the standard X-11 method. The X-11 census method uses a ratio to moving-average technique (Shiskin et al. 1967, Ladiray and Quenneville 2001). We estimated the 15-y trends for each series with simple linear regression, including analysis of residuals, on the seasonally adjusted data. The residual analyses included visual confirmation that the assumptions of normality, homoscedasticity, and independence of residuals were not violated for any series. Many authors use the Mann–Kendall nonparametric approach to stream data of this type, but the parametric approach we used is as powerful as the Mann–Kendall approach when the data are not skewed (Onoz and Bayazit 2003) (confirmed by the residual analyses) and is more easily interpreted.

We determined relationships between quarterly biological indicators and climate variables with Spearman's rank correlation coefficient (Mazor et al. 2009). We were unable to model the potential autocorrelation between consecutive years in each season adequately because of small sample sizes. Therefore, we interpreted p-values conservatively and applied a Benjamini–Hochberg procedure (Benjamini and Hochberg 1995) to control the false discovery rate (FDR), rather than family-wise error rate, for multiple tests.

We analyzed the differences in bioindicators between test and reference sites among years and streams with mixed models, which we chose because they use maximum-likelihood estimates rather than least-squares estimates and are superior to traditional linear models when analyzing repeated-measures data with missing values (which was occasionally the case for macroinvertebrate data from winter). Preliminary models included season as a factor, but model complexity and interaction effects involving season and the other variables made interpretation difficult. Therefore, we fitted a separate model for each season. Models included stream (5 levels), year (15 levels), and site type (test or reference) as fixed factors with sites subjected to repeated measures. We included all interaction terms in the model and used residual analysis to assess the assumptions of normality and constant variances before interpreting the outputs for each model. We present all results, but the main tests of interest were the interactions between year and stream or site type because these effects would indicate the potential for confounding by climate-induced changes. We controlled the FDRs for the tests with the Benjamini–Hochberg method (Waite and Campbell 2006) and interpreted significant effects by comparing least squares means where they could be estimated.

Climate trends and correlation with macroinvertebrate data

Highest monthly temperature, lowest monthly temperature, and number of days with rain have increased significantly at KNP weather stations since 1994. Mean 3-PM relative humidity has decreased significantly in the same period (Table 1). These changes were all highly significant (p < 0.0001) and, therefore, cannot be dismissed as chance effects or Type I errors even though the magnitude of change is small. Some climate variables and macroinvertebrate bioindicators were correlated before correcting for false discovery (Table 2), but the relationships were inconsistent. Some climate variables were both positively and negatively correlated with bioindicator values (depending on the site) (Table 2). For example, EPT richness ratio and Oligochaeta abundance were strongly related to climate variables at 7 sites. None of these relationships were significant after Benjamini–Hochberg correction.

Table 1

Annual changes between 1994 and 2008 and significance of regression (p) after seasonal adjustment for 5 climate-related variables in Kosciuszko National Park, Australia.

Table 2

Significant (p < 0.05) Spearman rank correlation coefficients between bioindicators based on macroinvertebrates collected from sites in Kosciuszko National Park and climate variables from nearby weather stations before Benjamini–Hochberg correction (after which no correlations remained significant). Seasons: 1  =  summer, 2  =  autumn, 3  =  winter, 4  =  spring. O/E  =  Australian Rivers Assessment System observed/expected taxa scores, SIGNAL  =  Stream Invertebrate Grade Number Average Level 2 scores, EPT  =  Ephemeroptera, Plecoptera, Trichoptera.

Macroinvertebrate indicators and effects of site type, stream, and year

Reference and test sites

Mean scores for all macroinvertebrate indicators differed significantly (p < 0.05) between test and reference sites in all seasons, except for EPT ratios in autumn (Table 3). Differences between test and reference sites remained significant after Benjamini–Hochberg correction (Fig. 2A–E, Table 3, 4). Mean scores for all bioindicators were lower at test than at reference sites, except for Oligochaeta abundance, which was generally greater at test sites. Mean O/E taxa scores for resampled reference sites were assessed as similar to reference condition (AUSRIVAS Band A) (Fig. 2A–C), except for the mean O/E taxa score in spring 2007, when it dipped below Band A (Fig. 3C).

Figure 2

Mean (±1 SE) values of Australian Rivers Assessment System observed/expected (O/E) taxa scores (A), Stream Invertebrate Grade Number Average Level (SIGNAL) 2 scores (B), Simpson's Diversity Index (C), ratio of Ephemeroptera, Plecoptera, Trichoptera (EPT) family richness (D), and Oligochaeta abundance (E) by stream, site type (test vs reference), and season averaged for all years (1994–2008). For ease of interpretation, untransformed data were graphed. However, the analysis was based on least-square means. In panel A, the dotted line represents the lower limit of Band A (0.87–1.14), which represents the reference condition. Numbers on the x-axis indicate seasons (1  =  summer, 2  =  autumn, 3  =  winter, 4  =  spring).

Figure 3

Least square mean values of Australian Rivers Assessment System observed/expected (O/E) taxa scores (A, B, C), Stream Invertebrate Grade Number Average Level (SIGNAL) 2 scores (D, E, F), Simpson's Diversity Index (G, H, I), ratio of Ephemeroptera, Plecoptera, Trichoptera (EPT) family richness (J, K, L), and log₁₀Oligochaeta abundance (M, N, O) in summer (A, D, G, J, M), autumn (B, E, H, K, N), and spring (C, F, I, L, O) at test and reference sites averaged across streams in Kosciuszko National Park, Australia, from 1994 to 2008. In panels A, B, and C, the dotted line represents the lower limit of AUSRIVAS Band A (0.87–1.14), which represents the reference condition.

Table 3

Probabilities associated with mixed-model tests for differences in values of bioindicators among site types (test vs reference), years (1994–2008), and streams for different seasons (1  =  summer, 2  =  autumn, 3  =  winter, 4  =  spring). Bold indicates p-values were significant after controlling the false discovery rate using Benjamini–Hochberg correction. O/E  =  Australian Rivers Assessment System observed/expected taxa scores, SIGNAL  =  Stream Invertebrate Grade Number Average Level 2 scores, EPT  =  Ephemeroptera, Plecoptera, Trichoptera.

Table 4

Mixed-model analysis of variance table for Australian Rivers Assessment System observed/expected (O/E) taxa scores in spring as an example to show structure and degrees of freedom. * indicates a significant p-value after correction for false discovery.

Climate changes and correlations

Over the 15-y study period, the climate became slightly warmer and less humid, but with a slight increase in the number of days with rain. The slight increase in the number of days with rain has not led to more rainfall and runoff for the region (BoM 2009). Our results indicated small but highly (statistically) significant changes in climate-related variables (Table 1), but the analysis did not conclusively establish a clear, consistent relationship between macroinvertebrate bioindicators and climate variables (Table 2). However, extreme climate-related events, such as the extensive bushfires of summer 2003, corresponded to changes in the bioindicator time series. Bushfires can release nutrients that are fixed in soils and plants with the result that nutrients are transported with sediments to nearby streams (Gresswell 1999). Such inputs to streams that are normally low in nutrients and sediment are likely to alter taxonomic composition. SIGNAL (Fig. 3D–F), Simpson's Diversity values at test and reference sites (Fig. 3G–I), and O/E taxa scores for reference sites (Fig. 3A–C), were lower after the fires, and Oligochaeta abundance increased in the summer and the following autumn (Fig. 3M–O).

Further severe climatic conditions occurred from May 2006 to December 2006 when most of Australia was strongly affected by El Niño. Large regions of southeastern Australia received the lowest rainfalls on record. These dry conditions, combined with little relief from the previous 2002–2003 El Niño, contributed to another extremely dry season in 2006–2007 that resulted in long-running bushfires in Victoria and New South Wales ( http://www.bom.gov.au/climate/enso/enlist/index.shtml). In the same period, poorer stream conditions were indicated by O/E taxa ratios (Fig. 3A–C), EPT richness ratios (Fig. 3J–L), and SIGNAL scores (Fig. 3D–F). Thus, bioindicators detected the effects of extreme climate-related events.

Reference condition

Climate change can trigger major environmental changes in aquatic ecosystems, which, in turn, trigger shifts in assemblage composition and structure of benthic macroinvertebrates and the bioindicator values used to assess stream condition (Harper and Peckarsky 2006, Durance and Ormerod 2007). Understanding whether bioindicators are changing in response to climate change is particularly important where the reference-condition approach (Reynoldson et al. 1997, Bailey et al. 1998) is used to assess sites at risk from human activities. Also important are robust bioindicators that can detect effects resulting from management-generated changes designed to enhance the condition of degraded streams. All bioindicators differed significantly between test and reference sites in all seasons (Table 3). The bioindicator values showed that, on average, stream conditions at test sites were most degraded at Perisher and Piper's Creeks in winter and spring compared to other streams and seasons (Fig. 2A–E). The most likely causes of this pattern are the pressures of ski-resort use and the treated effluent discharged into streams (Cullen et al. 1992, Digance and Norris 2006).

The average values of each bioindicator differed among some years for some seasons regardless of test- or reference-site status, but only O/E taxa scores at the resampled reference sites could be used to determine that the reference condition was still valid and had not changed consequent to climate-related changes. This determination was possible because a subset of the original reference sites used to define the reference condition of the predictive model had been sampled through time (Fig. 3A–C). This result suggests that, apart from extreme climate-related events such as drought and fire, the bioindicators tested here were not sufficiently sensitive to detect significant, but subtle, changes in climate during the 15-y study (but see Lawrence et al. 2010, Poff et al. 2010, Stamp et al. 2010 for more information on sensitive indicators). However, we think that clear shifts in these bioindicators are likely if the predicted changes to the frequency of extreme weather events (DECC 2008) and snow cover in the Australian Alps occur (96% reduction in the area covered by snow for >60 d by 2050).

Differences among streams

Macroinvertebrate communities are prone to temporal and spatial variation in natural and modified streams (Linke et al. 1999, Bêche and Resh 2007). Many factors, such as differences in substrate (Quinn and Hickey 1990), stream size and position in the catchment (Minshall et al. 1985), in-stream and riparian vegetation (Read and Barmuta 1999, Quinn et al. 2004), temperature (Jacobsen et al. 1997) and flow-related variables (Boulton et al. 1992, Scarsbrook 2002, Armitage 2006, Poff and Zimmerman 2010), contribute to variation in the distribution and abundance of benthic macroinvertebrates. All macroinvertebrate bioindicators tested, except O/E taxa scores, differed among streams in all seasons (Table 3). O/E taxa scores are inherently adjusted for individual site characteristics (Reynoldson et al. 1997, Simpson and Norris 2000). The O/E taxa value is derived by comparing the taxa observed at a site with those expected based on the taxonomic composition of macroinvertebrate assemblages from groups of reference sites (least disturbed by human activity). When making an assessment, the test site is allocated a probability of belonging to each group of reference sites based on the values of model-specific environmental characteristics (Wright 1995, Simpson and Norris 2000). The reference group that is most similar to the test site will contribute most to determining the expected assemblage at each test site. Thus, streams in KNP might be in similar condition overall, but have among-stream differences in macroinvertebrate assemblages that are related simply to the situation of the streams in the landscape and site-specific differences in macroinvertebrate habitat variables (Hawkins et al. 2010). Variation among streams in the EPT richness ratio and other richness-related indicators also can be caused by stressors that exclude sensitive taxa. The tolerant taxa that remain might differ by site or stream as a function of habitat-specific requirements (Bady et al. 2005, Pollard and Yuan 2006, Hawkins et al. 2010). Our results support the view that stressed systems may have macroinvertebrate assemblages that are unique to each stream.

Concluding remark on bioindicators, study design, and climate change

The AUSRIVAS O/E taxa score (which uses the reference-condition approach and predictive modeling for bioassessment) was the only indicator in our study that was adjusted for specific characteristics of each stream, and, thus, was robust to stream differences while being sensitive to reference- and test-site variation. Furthermore, O/E taxa scores indicated that reference sites sampled through time have remained in a condition equivalent to the previously defined reference condition, and, thus, provide a valid comparison for current test sites of unknown condition. Other bioindicators (e.g., EPT richness ratio, SIGNAL, Simpson's Diversity Index) varied significantly among streams, a result indicating that care should be taken when making comparisons among sites in different streams based on indictors not adjusted for natural stream differences.

The effects of climate change can be seen in 2 main ways: 1) gradual changes in temperature, rainfall, runoff, fogs, etc., and 2) changes in the frequency and severity of extreme events, such as drought and bushfires (Hughes 2003, Campbell 2008). The bioindicators used in our study were insensitive to gradual but significant climate changes that occurred over 15 y. Generally, reference and test sites did not respond differently to gradual climate changes and the differences between them were consistent over time. However, extreme climate-related events (such as severe drought and extensive bushfires) could be detected at both reference and test sites by the bioindicators. Appropriate study designs that include standardized sampling of fixed sites (both test and reference) over long periods are needed to enable detection of climate-change effects, regardless of whether they are gradual changes or extreme events.