Drivers of macroinvertebrate community structure in unmodified streams

Study sites

The Tongariro National Park in the central North Island of New Zealand comprises a central volcanic massif of three mountains of predominantly andesitic geology and the Tihia-kakaramea volcanic massif to the north. For sampling, 47 first- to sixth-order streams and rivers were selected from within the Tongariro National Park that were subjected to minimal human interference by excluding sites with greater than 10% catchment pastoral land use. All sites had a minimum of 90% volcanic hard sedimentary geology. Ten of the 47 sites had flows modified by hydroelectric dams, but as these were run-of-river type dams, they were not able to hold back floods and thus kept relatively natural flow regimes. The flow regimes of the remainder of sites were unmodified and ranged from runoff-fed streams to stable spring-fed streams. Eight sites were situated within Pinus radiata plantation forestry, but these sites were limited to mature forest, which has been found to have similar macroinvertebrate communities to native forest (Quinn et al., 1997; Quinn, Boothroyd & Smith, 2004). These sites belong to eight fifth order drainages (Fig. 1). For a more detailed description of the study sites, see Tonkin, Death & Collier (2013).

Biological collections

Physicochemical sampling

Physicochemical assessment was performed concurrently with biological sampling and sampled variables can be found in Table 1. Bed stability was assessed using the substrate component of the Pfankuch stability index (Pfankuch, 1975), which consists of rock angularity, brightness, packing, percent stable materials, scouring and amount of clinging vegetation.

Substrate size composition was assessed using the ‘Wolman Walk’ method, by selecting and measuring (beta axis) 100 stones at 1-m intervals 45° to the stream bank in a zigzag manner (Wolman, 1954) and grouping into Wentworth scale size classes (Wentworth, 1922; Cummins, 1962). The percentage of these classes was then converted into a single substrate size index (SI) by summing their midpoint values weighted by their proportion (I assigned bedrock as 400 mm for use in the calculation) (Quinn & Hickey, 1990). Stream-slope was assessed as the height drop over 10–100-m sections depending on the size of the stream. Substrate heterogeneity was calculated using the Shannon diversity index, and embededness graded on a three-point scale from loose to tight.

Conductivity, temperature and pH were spot-measured using a Eutech ECScan and pHtestr2 (Eutech Instruments, Singapore) respectively. Depth and current velocity were measured using a Marsh McBirney flomate current meter (Marsh McBirney, Frederick, Maryland) at five equidistant points along the thalweg, and width at three points, of each study reach.

The percentage of coarse particulate organic matter (CPOM) was visually assessed in the substrate along five transects throughout each reach. The percentage of debris jams and undercut banks were visually assessed along the entire study reach, as well as the percentage of overhead cover shading the stream. The percentage of riparian vegetation was visually assessed along each study reach broken into the following categories: native forest, native scrub, planted forest, pasture and bare ground.

Finally, to assess the link between overhead cover and catchment bare ground and tussock cover, catchment land-use for each site was extracted from the River Environment Classification (REC) (Snelder & Biggs, 2002), and % bare ground and % tussock cover was combined for analyses.

Data analysis

All analyses were performed using R version 2.15.2 (R Core Team, 2013). To assess the strength of the link between stream width and % overhead cover, and whether this link was confounded by elevation, Pearson’s correlation was performed using the cor.test() function. Further, to test if this link was reflected by the catchment cover of tussock and bare ground, overhead cover was correlated with the combined tussock and bare ground percentage.

Spatial autocorrelation was examined using a Mantel test based on Pearson’s product-moment correlation and 999 permutations, with the mantel() function in the package Vegan (Oksanen et al., 2011). This was performed by comparing geographic with environmental and macroinvertebrate dissimilarity matrices calculated using the function vegdist() in Vegan. Geographic and environmental matrices were calculated using Euclidian distances and macroinvertebrate matrices were calculated using Bray-Curtis distances on log (x + 1) transformed macroinvertebrate data. Geographic data was simply the NZTM (New Zealand Transverse Mercator) easting and northing coordinates.

To explore the best set of environmental drivers of community metrics, partial least squares regression (PLSR) was performed using the plsr() function in the pls package (Mevik & Wehrens, 2007). First, the environmental variables were standardised using the scale() function. PLSR was then performed using the full set of 24 standardised environmental variables and models were validated using leave-one-out cross validation. This was carried out for each of the eight macroinvertebrate metrics separately to get the best possible prediction for each individual metric and the number of components was limited to two. For the cross-validated models, the root mean square error of the estimate was calculated to evaluate the model.

For all multivariate analyses, the raw invertebrate data was log (x + 1) transformed to reduce heteroscedasticity. To select the best subset of environmental variables explaining variation in the multivariate macroinvertebrate community structure, BIO-ENV (Clarke & Ainsworth, 1993) was performed with bioenv() in the Vegan package. This function selects the best subset of environmental variables by maximising the correlation between environmental (using Euclidean distances) and community distance matrices. Spearman rank correlations were used to correlate the matrices and Bray-Curtis distances for the community dataset.

To visually assess the multivariate structure of the macroinvertebrate community, non-metric multidimensional scaling (nMDS) ordination was performed using the metaMDS() function in the Vegan package. Again, Bray-Curtis distances were used and the number of axes was limited to two. To examine the gradient of effect of the variables identified in BIO-ENV, smooth surface thinplate splines were fitted using the ordisurf() function in Vegan. This procedure uses generalized additive models (GAMs) to overlay a smoothed response surface, which allows a more detailed interpretation than a simple linear vector.

Environmental variables

Conductivity ranged from 40 to 298 µS cm^-1, averaging 112.8 µS cm^-1 (Table 1). Spot temperature ranged from 6.6 to 17.6 °C, with a mean of 10.79 °C. Velocity ranged from 0.16 to 1.46 m s^-1, with a mean of 0.78 m s^-1 (Table 1). Mean depth was 27.19 cm and varied between 5.7 and 52.2 cm. Substrate size varied between 43.67 and 254.64 mm, with a mean of 143.97 mm.

Mean chlorophyll a was 1.87 µg cm^-2 and ranged from 0.03 to 5.02 µg cm^-2 (Table 1). Thin films were the dominant growth form of periphyton with a mean of 46.98% cover, followed by mats with 14.26% cover (Table 1). Filamentous algae was rarely present with a mean of 2.51% cover. Bryophytes were highly variable ranging from 0 to 90% cover, with a mean of 13.57% cover. Environmental data were not spatially autocorrelated (r = 0.03, p = 0.2), but there was a spatial association for the macroinvertebrate data (r = 0.08, p = 0.02).

Univariate metrics

A total of 97 taxa was collected from the 47 sites in this study. Insects dominated all the sites and were the most taxonomically rich, with 35 caddisfly (Trichoptera) taxa, 22 Diptera, 14 mayflies (Ephemeroptera), and eight stoneflies (Plecoptera).

The pooled number of invertebrates in the benthos ranged from 15 to 5978 individuals 0.5 m^-2 and the number of taxa collected at each site ranged from 4 to 45 taxa 0.5 m^-2 (Table 1). MCI ranged from 90 to 134.5, averaging 113.5, and QMCI ranged between 3.07 and 7.75, averaging 5.95. Percent EPT taxa and individuals averaged approximately the same, between 55 and 59%, but the range was much greater for %EPT individuals (7.85–90.85) than %EPT taxa (35.29–73.08; Table 1). There were, on average, 15 EPT taxa per site, ranging between 2 and 27 taxa.

PLSR

Key variables

The variables with the consistently highest (or most negative) loadings on either component across all eight macroinvertebrate metrics were % overhead cover, % bryophyte cover, % debris jam and depth (Fig. 2).

Percent overhead cover was one of the strongest predictors of the macroinvertebrate metrics (Fig. 2). Overhead cover had a loading greater than 0.4 on all of the first PLS components for all metrics except for density (richness = 0.49, MCI = 0.48, QMCI = 0.48, %EPT taxa = 0.52, %EPT individuals = 0.45, EPT richness = 0.5, Margalef’s = 0.5; Fig. 2) and % CPOM cover contributed slightly lower positive loadings to the same component (richness = 0.36, MCI = 0.37, QMCI = 0.37, %EPT taxa = 0.38, %EPT individuals = 0.35, EPT richness = 0.35, Margalef’s = 0.36; Fig. 2). On the same components that % CPOM and overhead cover were positively loaded, width exhibited weaker negative loadings between −0.29 and −0.36, as width and overhead cover were negatively correlated (r = −0.45, p = 0.001). However, overhead cover also declined with elevation (r = −0.33, p = 0.025), but there was no link between elevation and width (r = −0.08, p = 0.58). Overhead cover was highly negatively correlated with the percentage of the catchment with bare ground and tussock cover (r = −0.73, p < 0.0001).

Chlorophyll a had the strongest loading for either component predicting macroinvertebrate density with a loading of 0.57 on component 1 and on the same component % planted forest contributed −0.42 (Fig. 2). Percent bryophyte cover (−0.43), % debris jam (−0.59) and % CPOM cover (−0.43) all had strong negative loadings on the second component predicting the number of individuals (Fig. 2). As well as a strong negative loading on component 2 of the density model, debris jam contributed strong negative loadings on all of the remaining metrics (richness = −0.42, MCI = −0.32, QMCI = −0.39, %EPT taxa = −0.33, %EPT individuals = −0.37, EPT richness = −0.39, Margalef’s = −0.4; Fig. 2). Moreover, in addition to a positive contribution by % overhead cover on component 1, % periphyton mat cover contributed strongly negatively to predict QMCI (−0.43) and %EPT individuals (−0.42; Fig. 2).

Percent bryophyte cover also contributed strongly negatively to component 2 of the MCI (−0.51), QMCI (−0.42), %EPT taxa (−0.43) and %EPT individuals (−0.4) models (Fig. 2). Percent periphyton film cover only contributed strongly to the second component of MCI (0.45) and Margalef’s index (0.41). The Pfankuch stability index contributed negatively to component 1 of taxonomic richness (−0.32) and density (−0.33), but contributed little to any of the other models (Fig. 2).

Multivariate

The best model using BIO-ENV selected eight of the 24 environmental variables (% overhead cover, Pfankuch, chlorophyll a, depth, temperature, % debris jam, % planted forest, % bryophyte cover) and had a correlation (rho) of 0.562 with the multivariate community structure (Table 3). However, there was little improvement in the link with community structure from the model with four variables (rho = 0.525). This model included % overhead cover, Pfankuch, % planted forest and % bryophyte cover.

The nMDS ordination on log (x + 1) transformed macroinvertebrate communities produced a reasonable ordination with a stress of 0.159 (Fig. 3). Overlaying this ordination with individual GAM fitted smooth surfaces of each of the eight variables selected by BIO-ENV using thinplate splines indicates the main influence of these important variables was along nMDS axis 1 (Fig. 3). Sites positively loaded on axis 1, tended to be more resource rich with greater levels of chlorophyll a, bryophytes, overhead cover and more stable (i.e., Pfankuch).

One site (site 32), which had very low abundance (four taxa and 15 individuals collected), loaded much more negatively on axis 1 than any other site (Fig. 3). This site was situated in a plantation forestry area, which may explain the spline loadings of % planted forest declining along this axis. At the other end of the scale, site 45, which exhibited the highest loading on nMDS 1, was a small spring-brook with a high percentage of bryophyte cover and overhead canopy cover, with a macroinvertebrate community that reflected these physical factors.

Removing this site from the nMDS increased the two-dimensional stress from 0.159 to 0.182. Furthermore, running the BIO-ENV with this site excluded reduced the link between environmental variables and the invertebrate community from rho = 0.562 to rho = 0.533, and the set of explanatory variables remained similar, with the notable exclusion of % planted forest. This BIO-ENV selected nine variables: % overhead cover, Pfankuch, chlorophyll a, depth, temperature, debris jam, substrate size, % bryophyte, and % macrophyte. Finally, excluding site 32 from the PLSR on abundance (which was the only invertebrate metric % planted forest influenced strongly) removed the influence of planted forest on density, but also reduced the amount of variation explained by the first two components from 64% to 55%.

Splitting the invertebrate community into EPT and non-EPT taxa indicated similar multivariate structure (Mantel r = 0.62, p = 0.001), and BIO-ENV revealed similar drivers, with 10 and 11 variables selected for each of EPT and non-EPT data respectively (Table 3).

Of the eight variables explaining multivariate structure, temperature was the only one that contributed little to univariate metric prediction. The biggest loading for temperature on any of the PLS components, was −0.21 on component 1 predicting macroinvertebrate density. The remainder of environmental variables contributed to both univariate and multivariate prediction.

Metric use

Contrary to my hypothesis, local-scale environmental variables were able to explain substantial variation in taxonomic richness of benthic macroinvertebrate communities in these pristine mountain streams with two components of a PLSR (70%). This is a relatively strong explanatory power for a set of streams with a highly constrained range of landscape-scale factors such as geology, climate, land-use and altitude. Nonetheless, there was considerable natural variation in environmental variables measured at these sites, which may explain the strong link with the macroinvertebrate community. In accordance with my primary hypothesis, the environmental drivers of univariate macroinvertebrate metrics and those of multivariate community structure were essentially similar. Where regional differences are present in environmental conditions, the regional species pool is an important driver of local species richness, with the ability to override local scale influences in streams (Heino, Muotka & Paavola, 2003; Tonkin et al., 2014). However, without these regional differences, the importance shifts back to local-scale drivers, especially if large-scale factors such as geology are limited as they were in the present study.

I used a limited range of community metrics, which could potentially mask important information, given metrics both respond differently to different stressors (Yuan & Norton, 2003) and only represent a small part of community dynamics. A multimetric approach may overcome this reduction of complex relationships by simultaneously explaining several aspects of community structure and increasing the likelihood of incorporating a wider range of responses to different stressors (Karr, 1999). However, this approach is dependent on metrics displaying different trajectories from each other, and the metrics used here all exhibited similar responses to each other (http://dx.doi.org/10.6084/m9.figshare.1053134). Thus, a much wider range of metrics than used in this study would benefit a multimetric approach, such as functional diversity or taxonomic distinctness (Clarke & Warwick, 1998).

Most of the variation explained in this study was by resource supply metrics, which likely reflects the fact that no anthropogenic stressors were present at these sites. The metrics used in this study are commonly applied in biomonitoring situations and have been demonstrated to respond to various anthropogenic stressors including land-use degradation, flow reduction and forest harvesting (Lenat, 1988; James & Suren, 2009; Reid, Quinn & Wright-Stow, 2010; Shearer & Young, 2011). Not surprisingly, given the lack of human influences that would increase nutrient or organic loads, the metrics with the weakest link with environmental drivers were those designed to respond to organic enrichments, namely MCI and QMCI.

Overhead cover and resource supply metrics

Several environmental variables strongly influenced both macroinvertebrate metrics and multivariate structure. Of these, the most important were percent overhead canopy cover, percent bryophyte cover, stream-bed stability, debris jam cover, periphyton cover and biomass, and stream size. In a previous study on these streams, Tonkin, Death & Collier (2013) showed that disturbance and productivity were important drivers in open- but not closed-canopy streams, despite canopy cover not affecting periphyton biomass. Due to a high density of cover in many New Zealand forests, canopy cover often strongly influences stream communities due to lowered periphyton standing crops (Winterbourn, 1990; Death & Zimmermann, 2005). However, this forest presence does not necessarily shift the functional composition due to a comparatively low level of allochthonous material entering these streams (Winterbourn, Rounick & Cowie, 1981). Canopy clearance can, however, have far-reaching consequences and can follow the subsidy-stress pattern proposed by Odum, Finn & Franz (1979); where the initial response is positive before a critical threshold is reached that leads to negative effects on biodiversity of stream biota (Quinn, 2000; Allan, 2004). This phenomenon may be related to a greater periphyton food supply for macroinvertebrates, in the absence of shading (Death & Zimmermann, 2005).

While we are dealing with pristine streams in the present study, lack of canopy is not usually an independent stressor, rather it is associated with complete shifts in land-use, which can come with much more intense pressures than simply clearing riparian vegetation (Sponseller, Benfield & Valett, 2001). Nonetheless, Tonkin & Death (2012) found similar relationships between productivity, disturbance and diversity when comparing streams from the region assessed in this study and those of a region dominated by intensive agriculture. As expected, there was a negative link between canopy cover percentage and stream width in the present study, but this link was not as strong as one would expect. I sampled a range of stream sizes, from first- to sixth-order, but the only first order stream was spring-fed and the majority of second order streams had large spring inputs.

The weakness of this canopy-width link results from several of the smaller streams lacking riparian vegetation, or if present, the vegetation often consists of small shrubs and tussocks, as was evident with the strong correlation between catchment bare ground and tussock and overhead cover. The importance of canopy cover on these communities is, therefore, not a function of an underlying stressor such as land-use change, but a natural vegetation phenomenon. Riparian vegetation can structure stream communities in several ways, including provision of habitat, physical structure and resource supply to stream dwellers, as well as filtering land-use impacts (Naiman & Henri, 1997). The importance of canopy presence in this study, therefore, likely represents several other important habitat factors.

Along with overhead cover, bryophyte cover was an important driver of invertebrate communities in both the PLSR and BIO-ENV approaches, but it did not correlate with any invertebrate metric. Moss was clearly linked with stream stability and likely also reflects the flow regime of these sites. Nineteen of the sites had no moss and a further 13 had 5% cover or less, therefore, it is likely a non-linear effect. Many environmental factors drive non-linear threshold responses in ecological communities (Dodds et al., 2010), and while the PLSR was able to extract the most important drivers of community metrics, I cannot make any inference regarding non-linear effects because PLSR assumes linear relationships between the independent and dependent variables.

Plantation forestry

The only variable to influence multivariate structure but not any of the metrics, other than abundance, was the percentage of plantation forestry in the riparian zone. Yet, despite being a local-scale measure, all plantation forestry sites were situated within one defined area consisting of mainly mature Pinus radiata. These sites could be simply comprised of a different suite of taxa, but not ultimately affect the invertebrate metrics such as taxonomic richness. However, previous studies have found similar macroinvertebrate communities between sites in mature Pinus radiata forestry and native forest in New Zealand streams (Quinn et al., 1997; Quinn, Boothroyd & Smith, 2004).

While a clear negative link with nMDS axis 1 was apparent for plantation forestry, this pattern was driven largely by one site with very low abundance, and thus forestry sites did not group separately from the remainder of sites in ordination space. Moreover, plantation forest percentage was not related strongly to any of the other environmental variables. In fact, removing this site from the analyses removed the influence of planted forest on both community structure and abundance. Thus it may be that this single site was inflating the influence of planted forest on these communities although, when excluded, the explanatory success of the models was reduced.