Inconsistent Relationships of Primary Consumer N Stable Isotope Values to Gradients of Sheep/Beef Farming Intensity and Flow Reduction in Streams

Stable isotope values of primary consumers have been proposed as indicators of human impacts on nitrogen dynamics. Until now, these values have been related only to single-stressor gradients of land-use intensity in stream ecology, whereas potential interactive effects of multiple stressors are unknown. It also remains unknown whether stable isotope values of different primary consumers show similar relationships along gradients of stressor intensities. We sampled three common invertebrate grazers along gradients of sheep/beef farming intensity (0–95% intensively managed exotic pasture) and flow reduction (0–92% streamflow abstracted for irrigation). The δ15N values of the three primary consumers differed substantially along stressor gradients. Deleatidium δ15N values were positively related to farming intensity, showing a saturation curve, whereas Physella snail δ15N values were negatively related to farming intensity and Potamopyrgus snail δ15N values showed no relationship. In addition, Deleatidium stable isotope values responded positively to flow reduction intensity, a previously unstudied variable. An antagonistic multiple-stressor interaction was detected only for the mayfly Deleatidium, which occurred in streams experiencing up to 53% farming intensity. The lack of consistency in the relationships of the most important primary consumer grazers along the studied gradients may reduce their suitability as an indicator of anthropogenic N inputs.


Introduction
Most ecosystems are exposed to multiple environmental stressors acting simultaneously [1,2] and combined multiple-stressor effects can be difficult to predict because they can be larger or smaller than expected based on the effects of the individual stressors involved [3]. Stable isotope values have long been considered indicators of anthropogenic disturbance [4,5]. In aquatic environments, the use of δ 15 N (a measure of the abundance of the heavier isotope 15 N relative to the lighter isotope 14 N) as potential indicators of land-use intensification and nitrogen enrichment has been suggested for streams and rivers [6][7][8], ponds and lakes [9,10], estuaries [11] and coastal waters [12,13]. Agricultural land use can impose a variety of stressors on stream ecosystems, including nutrient enrichment, increased sediment load, higher light availability and augmented water temperatures due to removal of riparian vegetation [14,15]. Another important stressor is water diversion for irrigation, which commonly reduces stream discharge and flow velocity and changes sediment and temperature regimes [16]. Any of these stressors may alter stream nitrogen dynamics, for example by flow reduction intensifying in-stream nitrogen transformation processes such as denitrification by increasing water retention times, respiration rates and occurrence of anoxic conditions [17]. However, the combined effects of multiple agricultural stressors on stable isotope values in stream ecosystems still need to be investigated. Delta 15 N values differ between nitrogen sourced from precipitation, fertilizer and sewage [18,19], and also due to subsequent fractionation during assimilation by plants and other nitrogen transformation processes [20][21][22]. Fractionation is the enrichment of one isotope relative to another through biogeochemical processes that discriminate against the heavier isotope. Further, gaseous nitrogen losses through denitrification and ammonia volatilisation can lead to even more elevated δ 15 N values of the nitrogen initially derived from both inorganic fertilizer and manure [23]. Catchment land-use intensity has been positively related to δ 15 N values of stream water [19,24,25], sediment [26], aquatic plants [27], consumers [6,7,28] and fish [6].
Among stream consumers, invertebrate grazers represent an important link for nutrient and energy transfer from primary producers to higher trophic levels. Moreover, they have been used as bioindicators because they are widely distributed in running waters and generally reflect the isotopic values of their periphyton food sources (after taking into account enrichment through trophic fractionation) [29]. Compared to the high spatiotemporal variation of periphyton stable isotope values [30], grazers provide a more time-integrated and habitat-integrated measure of nitrogen sources and transformation processes [5,31,32].
Most of the existing stream studies have reported positive linear relationships between δ 15 N values and catchment land-use intensity [18,28], although Larson et al. [8] reported a saturation curve for δ 15 N values along a wide gradient of % agriculture. To our knowledge, no studies have investigated the combined effects of multiple, simultaneously operating agricultural stressors on consumer stable isotope ratios, even though non-additive effects such as antagonisms and synergisms have been documented from experiments assessing changes in diversity and ecosystem processes in freshwater ecosystems [33]. Nor have previous studies focused on multiple primary consumers: most have investigated the relationship between a single primary consumer taxon or calculated δ 15 N values from several taxa combined where the focal consumer taxon was not present at all sites [7,8]. Whether different primary consumer taxa show similar relationships along gradients of anthropogenic disturbance remains to be investigated.
Against this background, we studied consumer stable isotope values at 43 stream sites spanning wide gradients of land-use intensity (sheep/beef farming intensity and streamflow reduction) within an agricultural river catchment, to test the following hypotheses:

1.
δ 15 N values of common primary consumers (grazers) will show similar patterns along the gradients of catchment land-use intensity because they all feed on one resource (periphyton) and are therefore ingesting food with the same isotopic composition; 2.
δ 15 N values are positively related to catchment land-use intensities because higher farming intensity and greater flow reduction increase the input of heavy isotopes from fertilisation and increase the intensity of nitrogen transformation processes [17]; and 3. δ 15 N values follow antagonistic response patterns along the gradients of land-use intensity and flow reduction because both stressors have strong individual effects and their joint effects cannot exceed 100% [3].

Study Sites
The 43 sites were selected within a single catchment, the Manuherikia River in Central Otago, which is among the driest in New Zealand. All sites experienced similar climatic conditions. The sites included 3rd-5th order streams (subcatchment sizes ranged from 2.7 to 367.9 km 2 ) chosen to provide wide gradients of both % Farming Intensity and % Flow Reduction (Figure 1, further site details in [34]). The catchment contains sites with high flow reduction intensity in areas of low sheep/beef farming intensity and vice versa because water for irrigation is transported in water races from one sub-catchment to another. Our strategic selection of sites along the two land-use gradients ensured that % Farming Intensity and % Flow Reduction were uncorrelated (R 2 = 0.03, Figure 2).  Land cover in the Manuherikia River catchment is dominated by sheep and beef-cattle pastures (intensively managed 'high-producing exotic grassland', 25%), extensively managed 'low-producing grassland' (24%) and native tussock grassland (44%) [35]. Our measure of sheep/beef farming intensity was the percentage of each sub-catchment covered in 'high-producing exotic grassland', assuming that this land use would impose the strongest effects on nitrogen dynamics. Land cover data were available from the Landcover Database II [36] and a detailed delineation of stream reaches and their upstream catchments from the River Environmental Classification [37]. Overall, agricultural intensification in Central Otago has been relatively recent and land-use intensity is still increasing, therefore land-use legacy effects such as those reported elsewhere e.g., [38] are unlikely to be relevant in our case.  Land cover in the Manuherikia River catchment is dominated by sheep and beef-cattle pastures (intensively managed 'high-producing exotic grassland', 25%), extensively managed 'low-producing Flow reduction intensity was calculated from modelled data because the Manuherikia River catchment is subjected to five major water abstraction schemes, several dams and more than 238 individual water takes for which it is unknown exactly how much and when water is abstracted [39](Otago Regional Council, personal communication). Further, there is a complex network of water races transporting water within and between sub-catchments, making it almost impossible to estimate natural and current river flows for all reaches within the catchment. Moreover, it was not feasible to install flow gauges at all study sites. Kienzle and Schmidt [39] modelled streamflows for the Manuherikia catchment as runoff for hydrological response units based on rainfall, altitude, soil and vegetation properties while taking into account quickflow storage, groundwater storage and snow storage, using the ACRU model (Agricultural Catchments Research Unit; University of Natal, South Africa). They verified their model at four gauging stations with a high accuracy (ranging from 1.5 to 5.5%) for simulating the monthly totals. These four gauges were situated in headwater streams similar to our sites in size and geomorphology, but without any water abstraction. We therefore calculated flow reduction intensity as the percentage in streamflow reduction between the Dryland Scenario (current landcover but no water abstraction) and the Current Scenario by using the mean stream flow during five irrigation seasons from 1999/2000 to 2004/2005.
In the studied subcatchments, sheep/beef farming intensity ranged from 0 to 95% and flow reduction intensity from 0 to 92%. Invertebrates were sampled once from each site under base-flow conditions in Austral autumn, between 21 March and 4 April 2011.

Field Sampling and Sample Processing
Benthic stream invertebrates were sampled using a 500-µm mesh kick-net in pool and riffle habitats following standard methods for semi-quantitative collections [40]. We identified pools as deeper stream sections with slow, calm flow and riffles as shallower sections with relatively fast, turbulent flow. For each sample, we took invertebrates from 10 locations (depending on stream size from three to five different riffles or pools). Invertebrates were preserved in 70% ethanol and stored for two years prior to analysis. Although preservatives such as ethanol can affect δ 15 N values, other studies have shown that the effects of ethanol on benthic macroinvertebrate stable isotope composition was minor and non-significant when compared to non-preserved material [41,42].
We focused on the most abundant grazer taxa in our study systems because we expected these to play a pivotal role in the local food webs. The three most common and widespread grazers at our sites were the snails Potamopyrgus antipodarum (present at 39 of 43 sites, comprising 21.1% of the total number of invertebrates counted and identified [34]) and Physella acuta (37 sites, 1.8%), and larvae of the mayfly Deleatidium spp. (21 sites with max. 53% farming intensity, 9.9%). For stable isotope analysis, we selected ten similar-sized individuals per taxon and sample (mean body length 3.8 mm for Deleatidium, mean shell lengths 2.3 mm for Potamopyrgus and 2.9 mm for Physella; measured to the nearest 0.1 mm under a dissecting microscope; Olympus SZ51, Olympus, Tokyo, Japan) and stored them in 90% ethanol in 5-mL glass vials. Snail shells and detritus were removed from individuals by soaking them in 1mL 1M HCl overnight and rinsing with deionised water. Gut contents were not removed due to the small body sizes and because of strong relationships between stable isotope values of primary consumer tissue and gut contents [43].
Grazers were dried at 60 • C for 48 h and subsampled (0.8 mg) into tin capsules. Nitrogen isotopes were measured by combusting all material to N 2 gas in an elemental analyser (Carlo Erba Instruments model NC2500, Milan, Italy). Gases were separated on a packed molecular sieve GC column and sequentially sent to an isotope ratio mass spectrometer ('20/20 Hydra', Europa Scientific, UK). Isotope values are reported as δ values (parts per thousand deviations from atmospheric N 2 standards): where R is the heavy-to-light ratio of the isotope. A subset of samples was analyzed in duplicate and showed a mean standard deviation of 0.15% .

Data Analysis
We investigated the relationships of δ 15 N values of the three grazer taxa with sheep/beef farming intensity and flow reduction using general linear mixed models combined with information-theoretic model selection [44,45].
We examined 16 competing models for each grazer taxon including the global model (intercept plus five predictors: the first-order terms % Farming Intensity (FI) and % Flow Reduction (FR), habitat (riffle/pool), the second-order term FI × FI and the interaction FI × FR), simpler versions of the global model with one or more predictors removed, and the null model (intercept only). Sample sizes for each taxon were equivalent to the number of sites where the taxon occurred (39 sites for P. antipodarum, 37 for P. acuta, 21 for Deleatidium). All taxa were present at more than 16 sites, allowing us to examine all 16 models for each taxon. Site ID was included as a random effect to account for the spatial non-independence of the predictor habitat. The second-order term FI × FI was included to detect potential saturation-curve relationships such as reported by Larson et al. [8]. The second-order term FR × FR was not included because we had no hypothesis for this non-linear relationship. If the interaction or the second-order term were retained, the lower-order terms were also retained. The predictor 'habitat' was included because stable isotope values may differ between riffle and pool habitats due to different flow velocity conditions see e.g., [46]. All variables were centered by subtracting the sample mean (to improve interpretability of regression coefficients) and scaled with two standard deviations (to allow use of regression estimates as effect sizes) [47].
All models were ranked according to their AIC c values (small sample unbiased Akaike Information Criterion) [48]. The top model set was chosen by selecting all models within ∆AIC c ≤ 2 of the best model. Table 1 also shows the relative Akaike weights and marginal R 2 values (variance explained by the fixed effects) for the chosen top models. All analyses were computed in R (version 2.15) [49] using the packages lme4 [50] and MuMIn [51]. Table 1. Standardized partial regression estimates (effect sizes, ES), 95% confidence intervals CIs, delta AIC (Akaike Information Criterion) differences to the model ranked most highly (∆AICc), Akaike weights (re-scaled), marginal R 2 values (variation explained by the fixed effects) of the final models for the relationships of δ 15 N with sheep/beef farming intensity and flow reduction. Effect size categories after Nakagawa & Cuthill [52]: weak > 0.10, moderate > 0.30, strong > 0.50. Estimates in bold indicate where predictors had a significant effect on the response variable (95% CIs not including zero). FI = % Farming Intensity, FR = % Flow Reduction, FI × FI = second-order polynomial terms and FI × FR = interaction terms.   Deleatidium δ 15 N values were best modelled by the single-stressor model including the first-order and second-order terms for farming intensity, explaining 75% of the variation (curvilinear relationship, best described as a saturation curve; Table 1, Figure 4a). The second-best model for Deleatidium δ 15 N values was a complex multiple-stressor model, in which the relationship with one stressor gradient depended on the intensity of the second stressor. This model explained 76% of the variation in the data, with strong positive effect sizes for farming intensity and the interaction term, plus a weak positive effect size for flow reduction (Figure 4b). The best models for Potamopyrgus and Physella δ 15 N values indicated no relationship with farming intensity or flow reduction (null model; Table 1). The second-best models for Potamopyrgus (Figure 4d) and Physella δ 15 N values (Figure 4e) both included a positive but non-significant relationship with flow reduction (both models only explained 6% of the variation).

Primary
The interaction term in the second-best model for Deleatidium δ 15 N indicated an antagonistic interaction where the relationship with farming intensity was weaker at high flow reduction and the relationship with flow reduction was also weaker at high farming intensity (Figure 4b).
The predictor 'habitat type' (riffle versus pool) was not retained in any of the final models. Deleatidium δ 15 N values were best modelled by the single-stressor model including the first-order and second-order terms for farming intensity, explaining 75% of the variation (curvilinear relationship, best described as a saturation curve; Table 1, Figure 4a). The second-best model for Deleatidium δ 15 N values was a complex multiple-stressor model, in which the relationship with one stressor gradient depended on the intensity of the second stressor. This model explained 76% of the variation in the data, with strong positive effect sizes for farming intensity and the interaction term, plus a weak positive effect size for flow reduction (Figure 4b). The best models for Potamopyrgus and Physella δ 15 N values indicated no relationship with farming intensity or flow reduction (null model; Table 1). The second-best models for Potamopyrgus (Figure 4d) and Physella δ 15 N values (Figure 4e) both included a positive but non-significant relationship with flow reduction (both models only explained 6% of the variation).
The interaction term in the second-best model for Deleatidium δ 15 N indicated an antagonistic interaction where the relationship with farming intensity was weaker at high flow reduction and the relationship with flow reduction was also weaker at high farming intensity (Figure 4b).
The predictor 'habitat type' (riffle versus pool) was not retained in any of the final models.  Table 1 for details). The plots encompass the range of stressor values across the study sites and riffle/pool habitats (note the plotted range of farming intensities is smaller for Deleatidium). The interaction plot for Deleatidium (b) shows that the relationship with % Farming Intensity is more positive at lower % Flow Reduction (for ease of interpretation, % Flow Reduction was labelled as a factor with 2 levels; grey < 50% Flow Reduction, black > 50% Flow Reduction).

Do Primary Consumers Show Similar Relationships along Stressor Gradients
Ecologists have assumed that primary consumers mirror the stable isotope values of their periphyton food sources [32,53], and stable isotope values of different primary consumer taxa have been used interchangeably to assess anthropogenic disturbances in some previous freshwater studies [7,8]. In contrast, our three consumer taxa followed different response patterns in terms of direction and strength of observed relationships along our focal gradients of stressor intensity (see discussion below). Thus, we reject our hypothesis 1 that primary consumers show similar stable isotope patterns along stressor gradients.  Table 1 for details). The plots encompass the range of stressor values across the study sites and riffle/pool habitats (note the plotted range of farming intensities is smaller for Deleatidium). The interaction plot for Deleatidium (b) shows that the relationship with % Farming Intensity is more positive at lower % Flow Reduction (for ease of interpretation, % Flow Reduction was labelled as a factor with 2 levels; grey < 50% Flow Reduction, black > 50% Flow Reduction).

Do Primary Consumers Show Similar Relationships along Stressor Gradients
Ecologists have assumed that primary consumers mirror the stable isotope values of their periphyton food sources [32,53], and stable isotope values of different primary consumer taxa have been used interchangeably to assess anthropogenic disturbances in some previous freshwater studies [7,8].
In contrast, our three consumer taxa followed different response patterns in terms of direction and strength of observed relationships along our focal gradients of stressor intensity (see discussion below). Thus, we reject our hypothesis 1 that primary consumers show similar stable isotope patterns along stressor gradients.
Deleatidium δ 15 N values showed a strong relationship with sheep/beef farming intensity, whereas the best models for Potamopyrgus and Physella δ 15 N values indicated that this predictor was not able to explain a sufficient amount of variation in the data. Streamflow reduction only featured in the second-best models for Deleatidium and the two snails. The tighter relationship for Deleatidium may partly arise because of its absence from sites with more than 53% farming intensity that may be subject to high nitrogen concentrations and sediment levels [34]. Such physicochemical conditions may provide unsuitable habitat for this mayfly larva, which cannot feed on thick layers of filamentous green algae and biofilm attached to fine sediment [54,55], but they may suit snails such as Potamopyrgus that can exploit filamentous greens as well as epipsammic biofilms [54,56]. Thus, Deleatidium with its sweep-like mouthparts feeds preferentially on erect, tall-growing benthic diatoms (epilithon) and long filamentous algae are rarely ingested [57,58]. By contrast, Potamopyrgus possesses a radula that can scrape off tightly attached biofilms, enabling them to remove biofilm from surface and subsurface sediments [54] and even from fine or coarse particulate organic matter [59]. It is likely that Potamopyrgus can access a larger range of food resources than Deleatidium, helping to explain the snail's success in so many different habitats and as an invader [57]. The snail Physella also has a radula, therefore its primary feeding mode is likely to be similar to Potamopyrgus. In support of this theory, the relationship between the δ 15 N values of the two snails from the same sites was tighter than their respective relationships with Deleatidium ( Figure 3). Other possible explanations for differences in δ 15 N values between Deleatidium and the two snails might relate to their respective physiology, mobility and life history.
A limitation of our study is the lack of data on consumer diet composition and the δ 15 N values of all potential diet sources. However, given that we collected invertebrate consumers from the same habitats, differences in consumer δ 15 N values were probably related to ingestion of different components of the periphytic community. Potential food sources such as diatoms, filamentous green algae and leaf material differed in their isotopic values in a pasture stream in Waikato, another region of New Zealand dominated by farmland [56].

Multiple-Stressor Patterns of Primary Consumer Stable Isotope Values
Although survey-based studies such as ours are less suitable for explaining cause-and-effect relationships than manipulative experiments [60], the observed relationships between stable isotope values and our catchment-scale variables can be regarded as robust for one of our three response variables, Deleatidium, for which we were able to explain up to 76% of the variation.
At the catchment scale, Deleatidium δ 15 N values were positively related to rising sheep/beef farming intensity and also to increasing flow reduction, but δ 15 N values of the snails were either not or only weakly related to rising farming intensity (only partially supporting hypothesis 2). For δ 15 N values of Deleatidium, which showed by far the strongest patterns of the three consumers (but only occurred at 21 of the 43 sites and only in streams with up to 53.4% farming intensity), the best model was an antagonistic one in which the response pattern along one stressor gradient depended on the intensity of the second stressor (in support of hypothesis 3).

Differences in Consumers' Relationships with Sheep/Beef Farming Intensity
Four processes may have contributed to the increase in δ 15 N values in Deleatidium with rising sheep/beef farming intensity: inputs of industrial fertilizer (initial δ 15 N values between +3% and −3% ; mean around 0% ), inputs of animal waste products (δ 15 N values from +35% to −15% ; mean around +10% ), nitrogen transformation processes such as denitrification and ammonia volatilization in agricultural soils and streams [18,61], and atmospheric nitrogen deposition. The latter probably plays a lesser role in New Zealand than in many other developed countries, but all three other processes likely contributed to the observed δ 15 N patterns. High-producing exotic grasslands used for sheep/beef farming in New Zealand receive not just animal manure, but also regular inputs of industrial fertilizers, via a process called "topdressing" (often done using small light airplanes). Further, field observations of the streambed substrata during invertebrate sampling, combined with a strong sulfur smell at some study sites (K. Lange, pers. comm.), suggested that anoxic conditions in the sediment, which increase denitrification rates in a zone of very low oxygen conditions [17], may have occurred at some of the 43 sites. Because our study was observational rather than manipulative and we had no information on manure or fertilizer inputs from the stream catchments, we cannot determine the relative importance of these three processes in determining the patterns in our δ 15 N data.
The curvilinear relationship of Deleatidium δ 15 N values with farming intensity in our study contrasts with the positive linear relationships reported in previous studies in New Zealand [7] and elsewhere [6]. In these other studies, Anderson and Cabana [6] focused on nitrogen load from manure while Clapcott et al. [7] studied agriculture, forestry and urban areas, contrasting with our investigation of the percentage of high-intensity farming in an agricultural catchment. Anderson and Cabana [6] considered agricultural prevalence as a land-use intensity gradient, but their gradient only extended to a maximum of 52% of the catchment area. The possibility that manure from different animals may also differ in isotopic composition could also add to the observed variation between studies. Perhaps we were able to detect a positive relationship approximating a saturation curve because our sites covered a wider gradient of land-use intensity than Anderson and Cabana [6]. Larson et al. [8] also reported a non-linear relationship between δ 15 N and catchment land-use intensity (% agriculture) that was best explained by a saturation curve, and suggested the non-linearity may have been caused by biotic nitrogen accumulation itself being non-linear. Under such circumstances, phosphorus might have been the limiting factor for primary productivity so that further nitrogen inputs could not be utilized and turned into biomass. It is also possible that the bioavailability of nitrogen differed among sources such as inorganic fertilizer, manure and humic substances.
Potamopyrgus showed no relationship with farming intensity and we speculate that these snails may have switched to other resources where tightly-attached epilithic algae were unavailable or where other resources were more common, such as filamentous green algae, macrophytes, detritus or heterotrophic biofilms, all of which likely had δ 15 N values differing from that of tightly-attached epilithic algae [61].
The weak negative relationship of Physella δ 15 N values with farming intensity was unexpected. One possible explanation for the contrasting relationships of Physella and Deleatidium δ 15 N values is exploitation of different microhabitats. The snail tended to be associated with streams having intermediate farming intensities, whereas Deleatidium abundances were associated with lower farming intensities [34]. Therefore, it was unsurprising that the taxa were sympatric in only a few sites (Figure 3). These sites may have provided a variety of habitats including rocky surfaces for Deleatidium and macrophytes for Physella. Future studies that would determine stable isotope composition of all basal resources at a given site would permit gaining a better understanding of the consumer's resource use, as well as assessment of microhabitat use.

Positive Consumer Relationships with Flow Reduction Intensity
The positive relationships of δ 15 N values with flow reduction intensity in our study could be due to streams subjected to high flow reduction experiencing higher denitrification rates because of reduced flow velocities. Denitrification in streams mainly occurs in oxygen-depleted zones of bed sediments, and denitrification rates can be increased by high rates of respiration and anoxic conditions following excess rates of primary production related to increased water retention times [17]. Alternatively, weaker dilution of the concentrations of nitrogen from anthropogenic sources due to reduced stream flows could also result in increased δ 15 N values in the tissues of grazers.
For Deleatidium, the intensity of sheep/beef farming was more important in determining stable isotope values than hydrological alteration. In a companion study of biological invertebrate traits in the same river catchment [35], we also found that farming intensity had stronger effects on habitat availability and flows of matter and energy than flow reduction. It is likely that agriculture imposes more direct stressors on stream environments (e.g., sediment and nutrient inputs, reduced shading) than hydrological alteration (e.g., increased retention time and reduced peak flows). On the other hand, our choice of study sites spanning 3rd to 5th order streams may have included some differences in flow regimes between individual sites, and these may have contributed to the overall variation in our non-manipulative study, weakening the observed flow reduction effects.

Antagonistic Stressor Interaction
The interaction term between sheep/beef farming intensity and flow reduction was retained in one of the final models for Deleatidium δ 15 N values. This interaction was classified as an antagonism because the combined multiple-stressor effects on stable isotope values were weaker than one would have expected based on their respective individual effects. In general, it is difficult to find evidence for synergisms if both stressors have large individual effects and the sum of the individual effects exceeds 100% [3]. In our study, this was the case for Deleatidium where individual stressors effects already caused a strong increase in δ 15 N values. Moreover, antagonisms can also pose significant management challenges since the intensities of all interacting stressors may need to be reduced to achieve substantial recovery [62].
The interaction between farming intensity and flow reduction for Deleatidium could be interpreted as follows: Overall, δ 15 N values of Deleatidium increased with farming intensity (the strongest effect in the model) due to reasons discussed earlier on. Low levels of flow reduction appear to have little effect on this positive relationship. At high levels of flow reduction, however, nitrogen inputs from adjoining pastures via the flowing water might be decreased to such an extent that the positive effect of increased farming intensity on δ 15 N values is weakened.

Conclusions
Based on our findings, we believe that the mayfly Deleatidium may be well-suited as a bioindicator in stable isotope studies on agricultural impacts because its δ 15 N values were strongly related to sheep/beef farming intensity. Moreover, the mayfly provides a time-integrated measure (1-2 generations per year) [63,64] and occurs throughout New Zealand [57,65]. Further, sampling of this mayfly in the field is easy because it often reaches high densities, and the techniques for stable isotope analysis are well established. Interestingly, in terms of explaining variation in our stream survey data, δ 15 N values of Deleatidium (R 2 = 0.76) performed better than structural measures of invertebrate community composition (e.g., taxonomic richness of pollution-sensitive mayflies, stoneflies and caddis flies, R 2 = 0.42) [34] when detecting impacts at the catchment scale.
In stream surveys or experiments focusing on land-use effects, the same primary consumers are unlikely to be present at all study sites or in all experimental units, especially if these sites/units span broad gradients of stressor intensities. We caution against extrapolating 'missing' values for one taxon with data from other taxa (as done in [7] and [8]). Instead we recommend establishing baseline values and closely investigating resource use and multiple-stressor relationships for several primary consumers, ideally focusing on taxa with strict dietary preferences, such as Deleatidium.
Finally, we observed one complex interaction between paired stressors (out of three possible cases) in our analysis. Such non-additive response patterns have also been found for other functional or structural metrics in freshwater [60,66], terrestrial [67,68] and marine environments [62,69]. Thus, our findings suggest that non-linear and non-additive responses of consumer stable isotope values to multiple stressors may be fairly common and should therefore be considered in future studies.