Littoral and benthic macroinvertebrate community responses to contrasting stressors in a large New Zealand river

Abstract

Large rivers present major challenges for biomonitoring and the development of biocriteria because of difficulties obtaining samples representative of the range of habitats present and a paucity of information on taxa responses to environmental conditions. We sampled littoral and deepwater benthic habitats to investigate macroinvertebrate community composition associated with organic, sediment and heated water discharges at 22 transects along a 73 km section of the lower Waikato River, North Island, New Zealand. Significant community differences in near-field exposure versus no exposure sites were found only in littoral habitats, partly reflecting changes in relative abundances of Crustacea and Oligochaeta. Littoral and benthic habitats supported different macroinvertebrate communities, particularly downstream of the sediment input and following the onset of a river-floodplain morphology which may have enhanced deposition of suspended sediment. Distance down river was a significant covariate for most metrics derived from benthic samples, but did not influence the relative abundance of most macroinvertebrate groups in littoral samples. Collectively, these results suggest that littoral habitats are likely to be more suitable than benthic habitats for biomonitoring since macroinvertebrate faunas were responsive to near-field stressor exposure and were less likely to be influenced by longitudinal changes in river conditions.

Keywords:

• biomonitoring
• metric
• organic
• sediment
• thermal
• Waikato

Introduction

Compared to wadeable streams, little is known of how macroinvertebrate communities in large rivers respond to human pressure or how they can be effectively used as biomonitoring indicators (Blocksom & Johnson 2009). Overseas studies have shown variable responses by river macroinvertebrates to anthropogenic pressures, partly reflecting the physiographic setting of sites, and the range of habitats and methods used to collect samples (e.g. Royer et al. 2001; Angradi et al. 2009; Schletterer et al. 2010). Furthermore, as large rivers are at the lower ends of catchments, they often integrate a complex array of upstream stressors which can interact in multiple ways making cause–effect relationships hard to disentangle (see Mathaei et al. 2010; Piggott et al. 2012). The absence of clearly defined stressor gradients and the general homogenisation of benthic communities are significant challenges for the development of large river biocriteria (Blocksom & Johnson 2009). Nevertheless, possible response relationships between indicators and stressors need to be quantified so that the results of monitoring programmes can be appropriately interpreted (Wessell et al. 2008).

Studies carried out in non-wadeable New Zealand rivers have indicated that: (1) taxonomic richness of caddisfly and dipteran larvae and the relative abundances of dipterans reflected catchment land use differences in moderate-gradient boatable rivers (Collier et al. 2013); (2) sampling of both deepwater and littoral habitats provided the most consistent response in predominantly hard substrate rivers (Collier et al. 2013); and (3) the number of mayfly, stonefly and caddisfly taxa recorded in littoral samples was related to catchment land cover in predominantly soft-bottomed non-wadeable rivers of the Waikato region (Collier & Hamer 2013). However, the latter authors found that macroinvertebrate land cover relationships were not evident in the Waikato River where multiple pressures were thought to have obscured responses, highlighting the need for more detailed analyses of macroinvertebrate community responses to individual stressors.

In the present study, we investigated whether community-level responses to stressors differed among habitats likely to be targeted for biomonitoring, acknowledging the clinal changes and progressive cumulative effects of human pressures that occur down river systems. Aquatic macroinvertebrates were sampled in littoral and deepwater benthic habitats and relevant physicochemical parameters were measured across 22 transects spanning a 73 km section of the Waikato River between Hamilton and Mercer, northern New Zealand, above and below discharges associated with municipal, dairy factory and meat-work waste treatment plants (‘organic’), suspended sediment from a large tributary (‘sediment’), and heated discharge of cooling water from a power station (‘thermal’). Our objectives were to: (1) determine the influences of expected stressor exposure on community composition on two occasions characteristic of different seasonal conditions; and (2) test whether macroinvertebrate communities differed among habitat types and therefore could exhibit differential responses to stressors in this large river. Poulton et al. (2003) recommended a similar approach comprising longitudinal multihabitat sampling at sites related to potential sources of impairment for bioassessment and biocriteria development on the lower Missouri River, USA.

Materials and methods

Study area and sampling sites

The study spanned a river distance of 73 km along the lower Waikato River, New Zealand's longest river at 442 km. Around 39% of the catchment area is forested, although agriculture dominates on the lower river (74% of catchment area downstream of the last dam) and introduced willows (Salix spp.) are the dominant riparian tree (Collier et al. 2010). Flow passes through eight hydroelectric dams, the last being some 30 km upstream of the first sampling transect within the city limits of Hamilton (population 132,000). At this point, mean annual discharge is 259 m³ s⁻¹, increasing to 422 m³ s⁻¹ near the mouth to the Tasman Sea (Brown 2010), above the threshold of 350 m³ s⁻¹ defining a ‘large river’ according to Nilsson et al. (2005). The area around Hamilton has a mild temperate and damp climate with mean annual rainfall of 1190 mm, most of it falling between June and September. Consequently, river flows tend to be higher in winter but vary throughout the year on a daily basis due to hydropeaking (mean annual diurnal water level range 0.6 m; Pingram et al. 2014). Mean annual air temperature is 13.7 °C, with daily maxima typically ranging from 22–26 °C in January to February, and 10–15 °C in July to August. River water temperatures in Hamilton can decline to below 10 °C in winter and frequently exceed 20 °C in summer.

Between Hamilton and Ngaruawahia, treated organic discharges enter the river from municipal, dairy factory and meat processing waste water treatment plants which collectively have consent to discharge up to 225 kg of phosphorus per day to the river (NIWA 2010). At Ngaruawahia, the Waipa River contributes between half to two-thirds of the flow to the lower river during peak discharge (Brown 2010), and supplies two-thirds of the suspended sediment load mainly derived from a large landslip and bank erosion (Hicks & Hill 2010). Consequently, average turbidity is around three times higher downstream of the Waipa confluence than above (Vant 2010), and the difference is most pronounced when high flows in the Waipa coincide with low flows in the Waikato River. Average water temperature increases by 1 °C between Hamilton and downstream of Mercer due to a combination of natural heating, warm tributary inflows and discharge of cooling water from the Huntly power station equivalent to 38 m³ s⁻¹ at full load (Rutherford 2010). Cooling water may be discharged into the river at no more than 10 °C above ambient. Maximum river temperature as a result of cooling water discharge is not to exceed 25 °C 1 km downstream, and temperatures within 20 m of the discharge are not to exceed 30 °C. The discharge was designed to mix cooling water with 40%–50% of river flow, enabling fish to migrate up the non-heated side of the river, thereby creating a heat plume that gradually mixes across the river.

Sampling was conducted across 14 transects in spring when collection focused above and below discharges, and across 22 transects in autumn when several extra downstream sites (T17–T22) were included to assess clinal changes (see Fig. 1). The organic reach comprised two contiguous sections of river approximately 10 km long. The upstream section received Hamilton City stormwater inputs which have short-term localised effects on the river (Hickey et al. 2001) and were unlikely to have strongly influenced the upstream sampling transects (T1–T3). Point source organic discharges occurred 1.5–3.0 km upstream of T4–T6. The sediment reach comprised one transect upstream of the Waipa confluence (T7) and three transects downstream at approximately 500 m intervals below the confluence (T8–T10). The thermal reach comprised one (spring) or three (autumn) upstream transects (T11–T13), and three transects downstream of the cooling water discharge (T14–T16) over a river length of < 1 km. Full transverse mixing of the thermal plume is thought to occur somewhere downstream of Meremere (Rutherford 2010).

Figure 1 Location of sampling transects (T1–T22) where benthic and littoral sampling was conducted around three stressor types on the Waikato River. The lower river (T17–T22) transects were sampled only in autumn.

Physicochemical measurements

Spot physicochemical measurements obtained during sampling were supported by water quality data collected as part of regular monitoring (Vant 2010; Tulagi 2012) or by modelling of water quality changes down river (NIWA 2010). For upstream and downstream sections of the organic reach, water samples (three samples per transect) were collected in 50 ml clean Vulcan™ centrifuge tubes for analysis of phosphorus and nitrogen during autumn when low flows prevailed and dilution of inputs would therefore have been less. All samples were sealed, placed on ice, and frozen upon return to the laboratory. Concentrations of total phosphorus (TP), dissolved reactive phosphorus (PO₄-P), total nitrogen (TN), ammonium (NH₄-N) and nitrate (NO₃-N) were analysed spectrophotometrically with a Lachat Instruments flow injection analyser (model QuikChem 8000 FIA+) following thawing. For total nutrients, a persulphate digestion step was used to oxidise nitrogen and phosphorus species prior to analysis (Ebina et al. 1983).

For the sediment reach, water samples were collected on two dates during autumn and spring across upstream and downstream transects (three samples per transect: near left bank; mid-river; near right bank) for measurement of turbidity (Hach 2100N meter, APHA method 2130B). Water temperatures above and below the Huntly power station were measured at: (1) sites not influenced by the thermal discharge (above and below a tributary draining a lowland lake, and the true right bank at similar locations to littoral sampling); and (2) three impact sites from the true left bank below the discharge near macroinvertebrate sampling transects. Temperatures were measured every 30 min over 10 February to 28 March (autumn) and 15 September to 1 November (spring) using Hobo TidbiT temperature loggers (Onset Computer Corporation, Bourne, MA, USA); at impact sites these were deployed at 10 m intervals along a 30 m long chain set perpendicular to the true left bank to provide spatial information on water temperatures prior to invertebrate sampling, whereas at control sites single loggers were tethered to weights. All loggers were calibrated before deployment over 3.5 days at room temperature; the largest difference in maximum temperature was < 0.3 °C and coefficient of variation for mean temperature was < 3%.

Macroinvertebrate sampling

Invertebrate samples were collected using similar effort from littoral habitats on true left and true right banks at each transect by sweeping and brushing accessible substrates into a 0.5 mm mesh D-frame net during seasonal baseflows. These near-shore samples were collected in water < 1 m deep and predominantly from macrophytes (Ceratophyllum demersum and Egeria densa), wood and submerged roots (mainly Salix fragilis but also Alnus glutinosa). Littoral samples were typically collected from at least three locations over c. 100 m along each bank and were pooled for processing. Deepwater benthic samples (n = 3 per transect) were collected at three equidistant points across the channel using an air-lift sampler (Duncan & Associates, Cumbria, UK) deployed from a boat (Pehofer 1998). This sampler consisted of sections of tube (10.5 cm internal diameter) linked together so that length approximated measured depth. Compressed air was forced through hoses which vented at the bottom of the tube that had been lowered onto the riverbed. The tube was moved up and down using similar effort at all sites, with the help of an operator holding the upper section on the boat. Dislodged bottom material was caught by the rising bubbles within the tube and belched into a Nytex sample bag (0.5 mm mesh) at the surface. This method sampled sand-gravel and pumice (with associated macrophytes if present) substrates effectively at water depths of > 1.5 m and up to c. 5.0 m (mean sampling depth at each transect ranged from 1.6 to 3.7 m). Air-lift samples were not collected from the thermal reach in spring due to equipment failure.

Invertebrate samples were preserved in c. 70% isopropyl alcohol. Processing consisted of spreading the sample across a white tray and randomly selecting grid squares from which invertebrates were picked by eye. Grid squares were processed sequentially until ≥ 200 invertebrates were obtained or the entire sample had been processed to provide taxa relative abundances (see protocol P2 in Stark et al. 2001). Unprocessed parts of the sample were searched for additional unrecorded taxa which were noted separately and allocated a value of 0.5 so that their presence was included in analyses but had only a minor influence on relative abundance estimates. Invertebrate identifications were based on Winterbourn et al. (2000; insects), Winterbourn (1973; molluscs) and Chapman & Lewis (1976; crustaceans). The level of taxonomic resolution was to genus for most insects and molluscs, and ranged from family to phylum for other groups. Sponges and bryozoans were allocated values of 0.5 to indicate presence but they were not enumerated due to likely fragmentation during sampling and processing.

Statistical analysis

To determine patterns in invertebrate community composition non-metric multidimensional scaling (nMDS) analyses were conducted on 4th-root transformed percent abundance data using Bray-Curtis similarity (Primer v.6.1.13; Primer-E Ltd 2009). PERMANOVA was performed to test for effects of sampling habitat (littoral, benthic), season (autumn, spring) and stressor exposure ranked from 0 (no exposure) to 3 (near-field exposure). No exposure sites were either upstream or on the opposite bank to the discharge, while near-field exposure sites were those most proximate to the discharge. Intermediate exposure levels were assigned ranks (1 or 2) based on knowledge of discharge dynamics and measurements of nutrients, turbidity or temperature across the channel downstream of discharges, under the premise that distance longitudinally and laterally from the discharge plume reduced exposure due to mixing. Samples with no invertebrates were excluded from these analyses, and one air-lift sample from T9 comprising 100% Naididae was also excluded because it was a major outlier in nMDS plots.

Pairwise comparisons by sampling habitat were used to test for significant differences between lowest and highest stressor exposure (0, 3) which was taken to infer an ecologically significant stressor response if relative abundances changed in relation to increasing exposure. SIMPER was used to identify taxa contributing most to the discrimination between no exposure (0) and near-field exposure (3) sites. We analysed the following response variables within habitats for T1–T16 (Fig. 1) using ANOVA with stressor exposure and season as factors and river distance as a covariate, following rank transformation (Conover & Iman 1981): total taxa and Diptera richness, and percentage abundance of Crustacea, Oligochaeta, Diptera, Trichoptera, and Ephemeroptera + Plecoptera + Trichoptera (excluding Hydroptilidae which can proliferate on filamentous algal growths) (%EPT-H).

Results

Physicochemical characteristics

Nutrient concentration, turbidity and temperature recorded at control sites (upstream and outside discharge plume influence) and impact sites (within influence of plume) are shown in Table 1. With the exception of NH₄-N, nutrient levels were significantly elevated when measured in autumn, and in the impact reach TP exceeded the concentration considered ‘satisfactory’ for limiting nuisance aquatic plant growth (Tulagi 2012). Elevated nutrient concentrations in the downstream organic reach were consistent with long-term monitoring results at Hamilton and Horotiu (Tulagi 2012). Turbidity was around two to three times higher at the sediment impact sites compared to control sites (Table 1) and values were below ‘satisfactory’ levels (see Tulagi 2012) in spring when flows were higher. These differences are consistent with more intensive measurements summarised in Vant (2010). In the thermal reach, water temperatures were on average 1–2 °C higher downstream of the powerstation discharge (see Collier & Hamer 2014 for further details). Continuous temperature records indicated maxima were higher in autumn, and that differences between ambient river temperatures and temperatures below the discharge were greater in spring (Table 1). Daily variability in water temperature was much greater in the thermal plume compared to the natural water temperature regime.

Table 1 Mean (SEM) conditions characterising stressors at sites influenced by organic, sediment and thermal discharges, and at control sites in one to two seasons.

	Control^1			Impact^2			Z
	n	Mean	SEM	n	Mean	SEM	Statistic
Organic (autumn)
Total P (mg L^−1)	9	0.035	0.002	9	0.058	0.004	3.83***
PO4-P (mg L^−1)	9	0.024	0.002	9	0.044	0.003	4.00***
Total N (mg L^−1)	9	0.275	0.003	9	0.285	0.003	1.67*
NH4-N (mg L^−1)	9	0.020	0.003	9	0.034	0.007	1.40 n.s.
NO3-N (mg L^−1)	9	0.134	0.010	9	0.181	0.008	2.61**
Sediment
Turbidity (NTU; autumn)	6	1.59	0.15	6	2.80	0.44	2.05*
Turbidity (NTU; spring)	5	3.66	0.19	6	10.58	1.50	4.09***
Thermal^3
Mean temp. (°C; autumn)	4	21.50	0.11	9	22.41	0.08	4.79***
90th percentile temp. (°C)		22.87	0.18		23.93	0.11	3.66***
Max. temp. (°C)		23.91	0.24		25.62	0.32	3.05**
Mean temp. (°C; spring)	5	14.78	0.07	8	16.43	0.14	7.86***
90th percentile temp. (°C)		16.40	0.15		18.10	0.12	6.30***
Max. temp. (°C)		17.84	0.13		19.84	0.13	7.69***

*P < 0.05; **P < 0.01; ***P < 0.001; n.s., not significant.

¹Organic reach: T1–T3; Sediment reach: T7 and T8–T10 true right; Thermal reach: T13 and T14–T16 true right.

²Organic reach: T4–T6; Sediment reach: T8–T10 inside plume; Thermal reach T14–T16 true left.

³Values refer to number of loggers recording at 30-min intervals.

Macroinvertebrate community composition

A total of 89 invertebrate taxa was recorded in the lower Waikato River in this study, with 62 found in benthic samples and 80 in littoral samples. Few taxa were specific to a particular habitat, with 7 and 11 taxa having more than one individual in only benthic or littoral samples, respectively. Benthic taxa richness was highest above the Waipa River confluence where substrates appeared generally coarser (mean per sample = 12 taxa), and declined to an average of four taxa per sample below this confluence. Overall mean taxa richness for all transects was 7 ± 0.6 (SEM) in air-lift samples and 17 ± 0.4 in littoral samples. Five of the nine introduced invertebrate species known to be present in the Waikato River (Collier & Hogg 2010) were recorded during this study: the snails Physa acuta (70 samples), keeled Gyraulus sp. (18 samples) and Planorbella sp. (one sample); the worm Branchiura sowerbyi (four samples); and the leech Helobdella europaea (one sample). The sampling also highlighted the occurrence of sensitive leptophlebiid mayfly, gripopterygid stonefly and conoesucid and hydropsychid caddisfly taxa along the entire study section, although their occurrence below the Waipa confluence was patchy. Of the 96 air-lift samples collected, five around Huntly and one near Ngaruawahia contained no invertebrates; of the remaining air-lift samples 27 contained less than 20 individuals. All of the 72 littoral samples collected contained more than 20 individuals.

Crustacea was the dominant macroinvertebrate group (42% of total numbers) for all littoral samples collected in autumn and spring combined, followed by Trichoptera and Diptera, respectively (Table 2). The amphipod Paracalliope fluviatilis made up three-quarters of the Crustacea collected and dominated littoral samples. Naididae made up 82% of Oligochaeta. Chironomidae made up 85% of Diptera, and Oxyethira 64% of Trichoptera, although the caddis Hydropsyche and Pycnocentrodes tended to be relatively more common upstream of the Waipa confluence (3% and 7%, respectively). Differences between littoral and benthic habitats were consistent when comparing samples collected in both seasons (10 transects) and all samples collected within seasons (22 transects in autumn and 14 transects in spring).

Table 2 Mean percent abundance of macroinvertebrate groups in littoral and benthic habitats sampled in autumn (22 transects) and spring (benthic 10 transects, littoral 14 transects).

	Autumn		Spring
	Benthic	Littoral	Benthic	Littoral
Ephemeroptera	0.01	0.01	1.59	0.04
Plecoptera	0.12	0.04	0.73	0.37
Trichoptera	8.22	12.64	5.02	13.09
Odonata	0.05	2.05	0.01	0.58
Hemiptera	0.00	0.00	0.00	0.07
Coleoptera	0.00	0.00	0.02	0.10
Diptera	15.48	3.90	15.95	14.57
Lepidoptera	0.01	0.10	0.00	0.00
Crustacea	10.89	60.22	36.22	59.18
Mollusca	2.70	7.60	4.22	3.19
Oligochaeta	55.61	4.96	35.20	6.17
Acari	0.67	0.46	0.28	0.34
Other	6.23	8.02	0.75	2.30

Non-metric multidimensional scaling analysis of taxa percent abundances in autumn (when all transects were sampled) highlighted the distinction between littoral and benthic samples (Fig. 2A). Although there was some evidence of a longitudinal shift in community composition, the greatest separation occurred between sites around Hamilton (T1–T6) and Huntly (T11–T15) along axis 2 of the ordination, with most other transects variably distributed between those two locations in two-dimensional space. This separation reflects the taxonomic differences among the groups mentioned above and also taxa-specific differences shown in Fig. 2B. A Gyraulus species, the flatworm Cura, cyclopoid crustaceans, the damselfly Xanthocnemis, the shrimp Paratya and the snail Physa tended to be more indicative of littoral habitats in the lower river, whereas the caddisflies Hydropsyche and Pycnocentrodes, the flatworm Prostoma, indeterminate worm species, orthoclad midges and Sinelobus shrimp tended to be more commonly associated with benthic samples from the upper part of the study section. The amphipod Paracalliope, a keeled Gyraulus species, and the hydroptilid caddis Paroxyethria were strongly associated with littoral samples, while Naididae worms tended to be more indicative of benthic habitats, particularly downstream of the Waipa confluence where sediments retrieved during air-lift sampling indicated bottom materials were finer. Littoral samples were more similar to each other (average Bray-Curtis similarity = 55.5) than benthic samples (38.8), while similarity between littoral and benthic samples averaged 30.3.

Figure 2 Non-metric multidimensional scaling plot of autumn percent abundance data for 22 transects sampled between Hamilton City and Mercer on the Waikato River. A, Transect numbers; B, taxa correlated with ordination axes (r_s > 0.5). Air-lift samples across transects were combined and samples with no taxa (n = 6) or only one taxon (n = 1) were excluded from this analysis. Littoral samples on the true right of the river were not directly exposed to the focal discharges.

Stressor responses

PERMANOVA indicated statistically significant effects on invertebrate community composition associated with sampling habitat and season for the organic and sediment reaches, but no effect of stressor exposure (Table 3). For the thermal reach there was a significant effect of habitat in autumn and a significant effect of stressor exposure and season for littoral samples (a fully factorial PERMANOVA was not possible because air-lift samples were not collected in spring) (Table 3). Interactions between exposure × habitat and exposure × season were not statistically significant. Pairwise comparisons of no exposure versus near-field exposure samples indicated no significant differences for benthic samples, but littoral samples from the organic (t = 1.54, P < 0.05), sediment (t = 1.80, P < 0.01) and thermal (t = 1.54, P < 0.05) reaches were all significantly different between stressor extremes. Four taxa in littoral samples were identified by SIMPER to contribute > 8% to the discrimination between these extremes: the common amphipod Paracalliope, Naididae worms, the hydroptilid caddisfly Oxyethira and the midge larva Tanytarsus. Of these, Tanytarsus was the only taxon to show a decline in relative abundance consistent with increasing stress exposure.

Table 3 Results of PERMANOVA testing for the effects of ranked stressor exposure, sampling habitat (littoral, benthic), and season (autumn, spring).

Source	d.f.	SS	MS	Pseudo-F	P(perm)
Organic
Exposure	3	8894	2964.7	2.59	0.055
Habitat	1	15791	15791	1.79	0.001
Season(Habitat)	2	18440	9218	8.31	0.001
Intensity × Habitat	2	2885	1442	1.26	0.289
Intensity × Season(Habitat)	5	5738	1148	1.03	0.408
Residual	46	51048	1110
Total	59	101000
Sediment
Exposure	3	3006	1002	1.10	0.343
Habitat	1	22124	22124	24.25	0.001
Season(Habitat)	2	8979	4490	4.92	0.001
Intensity × Habitat	3	3067	1022	1.12	0.334
Intensity × Season(Habitat)	3	2531	844	0.92	0.517
Residual	25	22812	912
Total	37	75018
Thermal (autumn)^1
Exposure	3	1140	379.91	0.50	0.86
Habitat	1	16199	16199	21.49	0.001
Intensity × Habitat	3	1706	568.78	0.75	0.577
Residual	17	12815	753.84
Total	24	49688
Thermal (littoral)
Exposure	2	1540	770	1.77	0.044
Season	1	3529	3529	8.10	0.001
Intensity × Season	2	1360	680	1.56	0.092
Residual	14	6100	436
Total	19	14309

Significant effects (P < 0.05) are indicated in bold.

¹Air-lift samples were not collected from thermal sites in spring so sampling habitat was tested only in autumn and seasonal effects were tested only for littoral samples.

Longitudinal discontinuities were evident for some of the seven structural metrics investigated (Fig. 3 richness and % abundance of Diptera not shown because no relationships were evident). Benthic taxa richness was very low around Huntly in autumn for transects above and below the thermal discharge, but then increased in some downstream samples to levels approaching those in the organic reach (Fig. 3). In contrast, the number of taxa in littoral habitats did not change noticeably along the section of river studied. %Crustacea declined markedly in littoral samples around Huntly in autumn, partly reflecting increases in %Trichoptera (mostly Hydroptilidae) and %Oligochaeta which also increased in relative abundance around the Waipa confluence. The increase in %Trichoptera for littoral samples around Huntly occurred for transects both upstream and downstream of the thermal discharge, although variability was high. %EPT-H abundances were low in littoral samples throughout the length of the river, but in benthic habitats were higher (> 20% of numbers) in some samples from the upper part of the river and comprised > 5% of numbers in some samples collected downstream of Huntly (Fig. 3).

Figure 3 Longitudinal changes in five macroinvertebrate compositional metrics with distance down river from T1 (see Fig. 1). Air-lift samples across transects were combined and samples with no taxa or only one taxon were included in this analysis. F-values indicate significance of stressor exposure (no parentheses) with distance down river as a covariate (parentheses). Arrow in upper left panel indicates location of the wetland tributary mentioned in the text.

For the structural metrics derived from T1–T16 samples, stressor exposure effects were statistically significant only for %Crustacea and %Oligochaeta in littoral habitats (ANOVA, P < 0.05; Fig. 3). Season was a statistically significant factor affecting %Crustacea in benthic habitats, total taxa richness in littoral habitats (both lower in spring), and %Trichoptera in benthic habitats (lower in autumn). There were no interactions between stressor exposure and season for any of the metrics tested in either habitat type. The covariate distance down river was significant for all metrics derived from benthic samples, including Diptera taxa richness (F = 37.04, P < 0.001) and percent abundance (F = 12.25, P < 0.01). For littoral habitats, distance down river was significant for %Crustacea, %Oligochaeta (Fig. 3) and Diptera taxa richness (F = 6.50, P < 0.05).

Discussion

Stressor responses

As with large rivers globally, the Waikato River is subject to multiple cumulative pressures, and disentangling the relative effects of individual stressors and their interaction with natural environmental gradients is difficult to achieve in the absence of comparable control rivers. Our approach to this issue was to differentiate stressor effects as much as possible by sampling across transects above and below large contrasting discharges, most of which were known to form plumes that would affect stressor exposure both longitudinally and laterally below discharges, thereby enabling us to sample near-field to far-field zones of influence. Clinal effects, such as decreasing riverbed sediment particle size, were expected to be more evident in benthic deepwater habitats than in littoral habitats, so congruence in responses between these habitat types would give some insights into the relative significance of stressor versus clinal influences. Sites not exposed to particular stressors occurred upstream of, or on opposite banks to, the discharges, although these ‘control’ sites would have experienced a mix of influences from upstream. However, other discharges to the river are small, with limited spatial and temporal influence relative to those studied (e.g. Hickey et al. 2001; NIWA 2010), and would therefore be expected to have considerably less effect on biota. Moreover, pest fish species presence was constant throughout the section of river studied, although koi carp are larger and more abundant below Ngaruawahia where they may have greater ecological effect (Hicks et al. 2010). Historical legacies can also influence the contemporary ecological potential of large rivers, and constrain biological responses to stressors. The Waikato River has experienced channel deepening in the upper part of the study section due to gravel dredging and bed degradation caused by retention of sediment by upstream dams, while in the lower section stop-banking and wetland drainage have all had enduring effects on natural ecological function (Howard-Williams et al. 2010). The results of this study must therefore be viewed within the context of these historical legacies and the prevailing mix of cumulative stressors.

Other studies have demonstrated that macroinvertebrates can respond directly to increases in water temperature (e.g. Quinn et al. 1994; Carolli et al. 2011), and suspended and deposited sediment (e.g. Ryan 1991; Larsen et al. 2011), while nutrient enrichment is more likely to have indirect effects through changes to aquatic vegetation cover and dissolved oxygen concentrations in off-channel areas (Houser & Richardson 2010). Work in experimental mesocosms to determine individual and combined effects of multiple stressors has shown that deposited sediment has a greater effect on invertebrates than nutrient enrichment or moderate temperature elevation (Piggott et al. 2012), although it is unclear how these findings would scale up to large rivers. Near-field effects of particular stressors were indicated in the present study for community composition but only in littoral habitats where changes in %Crustacea and %Oligochaeta were associated with stressor exposure, at least for sediment and thermal discharges. The muted response of other macroinvertebrate groups in the present study could partly reflect: (1) stressors not being sufficiently intense or well-resolved to detect responses; and (2) historical legacies, habitat constraints, clinal changes or geomorphic discontinuities limiting the range of colonists that could potentially be affected at particular sites. Angradi et al. (2009) reported contradictory reach-specific responses of some macroinvertebrate metrics on mid-continent US rivers, apparently reflecting variability in prevailing stressor gradients among rivers and highlighting the need to examine stressor-specific impacts to resolve biotic patterns in large rivers.

Spatial and temporal patterns

Macroinvertebrate communities in the section of Waikato River studied were characterised by a few dominant taxa, particularly in littoral habitats dominated by organic substrates (willow roots, wood and macrophytes) and below the confluence with the Waipa River. Upstream of this confluence, water quality is higher (Tulagi 2012) and air-lift sampling indicated bed substrates are more compacted and often dominated by gravels, providing habitat for a wider range of benthic taxa. Communities colonising littoral and benthic habitats were most similar in this upper part of the river compared to downstream suggesting that the pool of local colonists influenced what was found in different habitats within this upper reach. Taxonomic composition throughout the river was more similar for littoral than benthic habitats, as indicated by within group similarity, highlighting the consistent physical structure of littoral habitats and changes in benthic substrate particle size down river evident from air-lift sampling. In support of this interpretation, distance down river was a significant covariate for community metrics derived from benthic samples but not for most littoral sample metrics. Collier et al. (2013) collected similar air-lift samples from a range of gravel-bed rivers around the country and retrieved a median of 417 individuals per sample (K Collier, University of Waikato, unpubl. data). We therefore consider that the low abundances in some benthic samples from the Waikato River reflect habitat suitability rather than sampling efficiency.

Delayed effects of tributary inputs mediated by river morphology were suggested by downstream changes in macroinvertebrate community composition being most pronounced several kilometres below the confluence of the sediment-laden water from the Waipa River, where the main channel widens as it enters the lower floodplain (Pingram 2014). Slower-flowing water in this floodplain section potentially enabled suspended sediment to settle out onto benthic habitats where it can influence invertebrate community composition (Poulton et al. 2003). This process may account for the clear discontinuity in macroinvertebrate community composition that was evident around Huntly based on the sampling carried out, irrespective of whether transects were above or below the thermal discharge. In addition, benthic taxa richness and percent abundances of littoral crustaceans were reduced, and there were corresponding increases in the relative abundances of oligochaetes in littoral samples, suggesting deposition of fine sediment along river edges.

Tributary junctions are known to create spatial discontinuities that enhance physical heterogeneity and introduce novel carbon sources leading to localised increases in biotic diversity in rivers (Kiffney et al. 2006; Rice et al. 2006; Grant et al. 2007). Jones (2010) proposed that waterway segments below lake outlets are likely to be more productive than lakeless systems because of resource subsidies, more so in small streams than large rivers where effects are likely to be most pronounced along edges downstream of confluences. The onset of the floodplain section in the lower Waikato River is also where several nutrient-rich riverine lakes discharge into the river, and growth of algae fuelled by these nutrient inputs may partly account for increases in littoral %Trichoptera abundance, which was largely a reflection of greater numbers of the algal-piercing caddis Oxyethira. The potential importance of tributaries was also highlighted by the patches of habitat near confluences, most notably for a wetland tributary (arrow in Fig. 3 upper left panel), that supported quite diverse macroinvertebrate communities downstream of Huntly, despite an overall decline in benthic taxa richness below the Waipa confluence. The Waipa River and nearby inflows from the forested Hakarimata Ranges may also contribute novel taxa to the river around tributary junctions (Collier & Lill 2008), although sampling in the present study was not intended to capture local diversity hot spots.

Macroinvertebrate community composition can reflect flow conditions as well as water quality in large rivers (Applegate et al. 2007), and higher flow levels during spring may partly account for observed temporal differences in community composition in the present study. Nevertheless, interactions of season with habitat or stressor exposure were all non-significant, suggesting consistent responses irrespective of high or low flow seasons. Flow history over shorter timescales could potentially influence colonisation of littoral habitats, with recently inundated areas less likely to be colonised by aquatic taxa. This would be expected to be a particular issue in sections of river influenced by daily hydropeaking operations which can lead to water level variations of up to 1.86 m per day (Pingram 2014). However, the consistent taxa richness in littoral samples for both seasons and across the length of the river suggests samples were collected at sufficient depth for any effect on results to be minimal in the present study.

Implications for large river biocriteria

Several workers have attempted to develop multimetric macroinvertebrate indices for large rivers, and have reported varying results depending on prevailing environmental conditions and stressors (e.g. Wessell et al. 2008; Angradi et al. 2009; Wiegel & Dimick 2011). A meta-analysis of studies investigating macroinvertebrate responses to pressure for rivers of order ≥ 5 identified a range of metrics indicating responses to human pressure (Table 4). Unique combinations of metrics were identified in many of those studies, suggesting variable responses among rivers depending on the physiographic setting, pool of colonist species, combination of stressors and the sampling methods or substrates used. EPT richness was used in eight studies but was not used in the present study because the EPT fauna was limited (mean of less than one taxon per sample, range one to four taxa), and previous work had shown that this metric did not reflect land use pressure on the Waikato River in the same way as other non-wadeable rivers in the region (Collier & Hamer 2013). Some midge groups such as Tanypodinae have been proposed as indicators because of a propensity to increase in abundance following disturbance (Royer et al. 2001; Applegate et al. 2007). The chironomid Tanytarsus responded consistently to the gradient of stress across all transects in the present study, but single taxa are generally not reliable indicators for biomonitoring purposes. We did not attempt to define functional feeding groups, tolerance levels or trait categories for macroinvertebrate taxa because of limited information on the biology of several large river invertebrate species, but this would be a fruitful avenue for further research aimed at developing biomonitoring indices for large rivers (e.g. Statzner et al. 2005).

Table 4 Details of various studies on 5th–10th order rivers investigating macroinvertebrate community responses to human stressors that have applied structural or functional metrics (excludes trait, ratio and tolerance indices).

Country/region	Dominant stressor	Habitat	Substrate	Metrics	Reference
USA, Idaho	n.d.	Benthic mid-channel	n.d.	Total richness, EPT richness, %Elmidae, %Dominant taxon, %Predators	Royer et al. (2001)
USA, various	Water quality	Near-shore	Masonite plates	Total richness, EPT richness, No. Diptera taxa, Total no. ind., No. Ephemeroptera ind., %EPT, %Hydroptilidae, %Diptera, %Tanypodinae, %Amphipoda, %Oligochaeta	Applegate et al. (2007)
USA, Michigan	Multiple	Across transect	All available	Total richness, Diptera richness, Plecoptera richness, EPT richness, %Trichoptera, % Dominant taxon, FFG diversity	Wessell et al. (2008)
USA, various	Multiple	Benthic littoral	n.d.	Multimetric index from 10 metrics that varied among reaches	Angradi et al. (2009)
USA, various	Agriculture	Across transect < 1 m deep	All available	Diptera richness, EPT richness, Coleoptera richness, %Oligochaete and leech taxa, %Predator taxa richness, %Collector-filterer	Blocksom & Johnson (2009)
USA, Missouri	Multiple	Benthic, groynes	Rocks, gravel-sand-silt	Total richness, EPT richness, Chironomidae richness, %EPT, %Ephemeroptera, %Chironomidae, %Dominant taxon	Jackson et al. (2010)
Russia, Tver	Organic pollution, pesticides	Benthic	Mud, sand, cobble, organic matter	EPT richness	Schletterer et al. (2010)
USA, Wisconsin	Multiple	Rock, snag	Masonite plates	Insecta richness, EPT richness, %Insecta, %3-dominant taxa, %Gatherer, %Scraper	Wiegel & Dimick (2011)
New Zealand, various	Agriculture	Littoral sweep, deepwater benthic	Sand-gravel, wood, macrophyte, bank	Total richness, Trichoptera richness, Diptera richness, %Diptera	Collier et al. (2013)
New Zealand, Waikato	Agriculture	Littoral	Wood, macrophyte, bank	EPT richness (excluding Hydroptilidae)	Collier & Hamer (2013)

EPT, Ephemeroptera, Plecoptera, Trichoptera; FFG, functional feeding group; n.d., no data. % = percentage of total numbers unless otherwise stated.

An issue affecting the development of ecological indicators for large rivers is the potential for some taxa to swamp macroinvertebrate metrics based on broad taxonomic groupings (Seegert 2000). These ‘swamper’ taxa can include Hydropsychidae on rocks in fast-flowing sections of river (Jackson et al. 2010), Oligochaeta in sections dominated by fine sediments (Battle et al. 2007), and Crustacea such as the migratory shrimp Paratya and the amphipod Paracalliope which can be highly abundant in littoral habitats (Collier et al. 2013; this study). Jackson et al. (2010) suggested that excluding dominant taxa from metric calculations may make observed responses more comparable to wadeable stream assessments, and Seegert (2000) proposed removing interruptive species that comprised more than 50% of the catch when developing fish metrics. We were unable to test this approach in the present study because rare taxa were not enumerated in the fixed count protocol and could not therefore be included in revised % abundance calculations, highlighting the need for full counts to provide more flexibility in the development of metrics based on large river macroinvertebrate communities.

Conclusions

Collectively, results of the present study indicate that littoral faunas were more responsive than benthic faunas to stressor exposure and were less likely to display clinal patterns. This finding supports the sampling of river edges during stable flow conditions for large river biomonitoring (see also Flotemersch et al. 2006). Benthic sampling, however, may assist in identifying sedimentation effects, and the selection of sampling methods and metrics used therefore needs to be context dependent to achieve the most reliable outcome for biomonitoring objectives. Interpretation of biomonitoring data from large rivers also needs to take account of geomorphic discontinuities, such as tributary junctions and floodplain morphology, which appeared to influence macroinvertebrate community patterns in this large river. Analysis of full-count samples in the future should provide more flexible datasets for interpreting stressor effects by allowing the exclusion of abundant generalist taxa that may swamp sensitive taxa responses to stressors.

Acknowledgements

This paper is dedicated the the memory of Stephen Moore – taxonomist and photographer extraordinaire. We thank Callum Bourke, Kris Tapeti, Bryan Clements, Nathan Singleton, Kane Lynn, Mike Pingram and Richard Glass for assistance on various aspects of this study. Mike Pingram and two anonymous reviewers provided helpful comments on the manuscript.