Assessing ecological community health in coastal estuarine systems impacted by multiple stressors

Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

Contents lists available at ScienceDirect

Journal of Experimental Marine Biology and Ecology journal homepage: www.elsevier.com/locate/jembe

Assessing ecological community health in coastal estuarine systems impacted by multiple stressors J.I. Ellis a,⁎, J.E. Hewitt b, D. Clark a, C. Taiapa c, M. Patterson d, J. Sinner a, D. Hardy c, S.F. Thrush e a

Cawthron Institute, Private Bag 2, Nelson, New Zealand National Institute of Water and Atmospheric Research Ltd., PO Box 11-115, Hillcrest, Hamilton, New Zealand Manaaki Te Awanui, 234a Waihi Rd, Tauranga, New Zealand d School of People Environment and Planning, Massey University, Private Bag 11-222, Palmerston North, New Zealand e Institute of Marine Sciences, Auckland University, Private Bag 92019, Auckland, New Zealand b c

a r t i c l e

i n f o

Article history: Received 23 December 2014 Received in revised form 11 August 2015 Accepted 3 September 2015 Available online 16 September 2015 Keywords: Canonical ordination Cumulative impacts Community health Multiple stressors Interactions

a b s t r a c t Increasing population pressure, urbanisation of the coastal zone and nutrient and sediment run-off from agriculture and forestry has increased the number of large-scale and chronic impacts affecting coastal and estuarine systems. The need to assess cumulative impacts is a major motivation for the current desire of managers and ecologists to define ecosystem “health” and “stress”. A number of univariate metrics have been proposed to monitor health, including indicator species, indicator ratios and diversity or contaminant metrics. Alternatively, multivariate methods can be used to test for changes in community structure due to stress. In this study we developed multivariate models using statistical ordination techniques to identify key stressors affecting the ‘health’ of estuarine macrofaunal communities. Macrofaunal and associated environmental samples were collected across 75 sites from within Tauranga Harbour, a large estuary located on New Zealand's North Island. The harbour receives discharges from urbanised, industrial, agricultural and horticultural catchments. Distance-based linear modelling identified sediments, nutrients and heavy metals as key ‘stressors’ affecting the ecology of the harbour. Therefore, three multivariate models were developed based on the variability in community composition using canonical analysis of principal coordinates (CAP). The multivariate models were found to be more sensitive to changing environmental health than simple univariate measures (abundance, species richness, evenness and Shannon–Wiener diversity) along an anthropogenic stress gradient. This multivariate approach can be used as a management or monitoring tool where sites are repeatedly sampled over time and tracked to determine whether the communities are moving towards a more healthy or unhealthy state. Ultimately, such statistical models provide a tool to forecast the distribution and abundance of species associated with habitat change and should enable long-term degradative change from multiple disturbances to be assessed. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The need to better understand the interactive and cumulative effects of multiple stressors is cited as one of the most pressing questions in ecology and conservation (e.g. Crain et al., 2008; Sala et al., 2000). With growing human population comes an increase in the diversity and intensity of anthropogenic stressors (Halpern et al., 2007) including habitat loss, over-exploitation of key species (Jackson et al., 2001) pollution (particularly excess nitrogen), invasive species and more recently climate change (Kappel, 2005; Sala et al., 2000; Venter et al., 2006; Wilcove et al., 1998). In New Zealand increasing population pressure, urbanisation of the coastal zone and nutrient and sediment run-off from agriculture and forestry has increased the number of large-scale and chronic impacts ⁎ Corresponding author. E-mail address: [email protected] (J.I. Ellis).

http://dx.doi.org/10.1016/j.jembe.2015.09.003 0022-0981/© 2015 Elsevier B.V. All rights reserved.

affecting coastal and estuarine systems. Potentially, this can lead to broader-scale changes to coastal and estuarine systems due to modification of habitats (Saiz-Salinas and Urkiaga-Alberdi, 1999; Smith and Kukert, 1996) and by influencing the health, abundance and distribution of benthic organisms. For example, elevated sediment loading to estuaries and coastal environments can lead to broad scale changes in ecology through modifying habitats (e.g., Saiz-Salinas and Urkiaga-Alberdi, 1999; Smith and Kukert, 1996) and, in particular, by influencing the health, abundance and distribution of benthic suspension feeders (Ellis et al., 2002). Benthic macroinvertebrates are frequently used as indicators of ecosystem health (Borja et al., 2000) because they provide and influence many critical ecosystem services including provisioning services such as food for birds, fish and humans and regulating and maintenance services such as primary and secondary production and the storage and cycling of nutrients (Thrush et al., 2013). Increased concentrations of silts and clay in suspension may significantly increase pseudofaeces production, decrease the amount

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

of algal food actually ingested, and may also damage bivalve gills (Bricelj and Malouf, 1984; Iglesias et al., 1996; Morse et al., 1982; Navarro and Widdows, 1997; Robinson et al., 1984; Stevens, 1987; Willows, 1992). Exposure to increased concentrations of suspended sediments for an extended time can, therefore, result in decreased amounts of energy available for growth and reproduction, and have deleterious effects on local populations. Ultimately this can result in functional changes including a loss of key suspension feeding species and a switch to deposit feeding communities. Low levels of nutrient enrichment in estuarine and coastal environments can have a positive effect on the benthos due to improved primary productivity, and therefore food availability. Beyond a critical point, however, excessive nutrient discharges can lead to accelerated eutrophication of coastal environments and adverse symptoms of over enrichment (Cloern, 2001; McGlathery et al., 2007). Metals can be essential for organisms as trace elements, however, at higher concentrations they can become toxic (ANZECC, 2000a). High exposure to heavy metals can cause physiological stress, reduced reproductive success, and outright mortality in associated invertebrates and fishes (Fleeger et al., 2003; Gagnaire et al., 2004; Nicholson, 1999; Peters et al., 1997; Radford et al., 2000). Estuaries and coastal ecosystems are particularly vulnerable to these stressors as they act as natural retention systems for sediments and heavy metal contaminants and are affected by nutrient run-off. The need to assess cumulative impacts is a major reason for the current desire of managers and ecologists to define ecosystem “health” and “stress”. Over the past decades, there has been a rapid development of indices of “ecosystem health” to assess the status of marine environments which are increasingly required under regulations like the US Clean Water Act and the EU Water and Marine strategy framework directives (de Juan et al., 2015). A number of univariate metrics have been proposed to monitor health including indicator species, indicator ratios and diversity or contaminant metrics (Borja et al., 2000; Gray and Mirza, 1979; Labrune et al., 2012; Margalef, 1958; Pielou, 1966; Rosenberg et al., 2004; Sanders, 1968; Shannon, 1948). In reviewing existing methods of defining and measuring ecological ‘health’ it was noted that many of the existing biological diversity indices do not differentiate amongst different types of taxa and are strongly affected by sample size (Dunn, 1994; Gappa et al., 1990). This limits their ability to detect changes in composition across different communities and habitats. Furthermore, it is not immediately apparent what differences or similarities in these indices actually mean to ecological functioning, as a similar diversity value can be obtained from communities with very different species (Clarke, 1993; Dufrene and Legendre, 1997). Alternatively, multivariate models that focus on community composition have been developed (see Anderson, 2008; Anderson and Willis, 2003; Hewitt et al., 2005). Community composition comprises both the number and type of taxa (or animals) that make up a biological community at a site, together with their relative abundances or biomasses. Defining community composition requires the same information needed to generate many univariate biological diversity indices; however, by preserving all the information on the abundance of specific taxa, a more sensitive, and more ecologically meaningful, response could be expected (Attayde and Bozelli, 1998; Gray, 2000; Hewitt et al., 2005, 2009; Pohle et al., 2001). The community composition found in areas largely unaffected by anthropogenic disturbances versus that found in more ‘impacted’ areas can be used as a benchmark against which to assess the relative health of community composition found at specific sites. Thus, relative ‘health’ can be defined in terms of the range of communities present in comparable locations that are not considered to be affected by anthropogenically-derived inputs and should serve to identify both acute effects and broader-scale chronic degradation. Community composition is generally determined using multivariate techniques, including ordination. Multivariate techniques have been applied successfully to indicate the effects of pollution (Ellis et al., 2000; Olsgard and Gray, 1995; Warwick et al., 1990) and subsequent studies have now shown that multivariate methods are better at determining

177

differences between communities with different degrees of contaminant disturbance than univariate measures (Hewitt et al., 2005). In the present study we examined benthic communities in a large harbour, which receives discharges from urbanised, industrial, agricultural and horticultural catchments. Therefore, we were interested in investigating whether multivariate methods would be useful to rank the health of intertidal sites based on predicted responses to multiple stressors including sedimentation, nutrients and contamination. The results of the ordination models were compared with some frequently used diversity indices. 2. Methods 2.1. Study site Tauranga Harbour is a large estuary (approximately 200 km2) located on the western edge of the Bay of Plenty on New Zealand's North Island. The harbour is predominantly shallow (b10 m deep), with intertidal flats comprising approximately 66% of the total area (Inglis et al., 2008). Catchment land use is predominantly pasture and indigenous forest with urbanisation concentrated in the south-east, near the city of Tauranga. Sampling was carried out from the December 2011 to February 2012. A total of 75 sites across the harbour were sampled for benthic macrofauna and associated sediment characteristics. Sites were chosen to reflect a range of habitats, including intertidal sand flats, shellfish beds, seagrass meadows and areas likely to be impacted by pesticides. Due to the large area sampled within a harbour, gradients in salinity and exposure are expected to be present with distance from river mouths. All sites were intertidal, therefore, depth and temperature were relatively constant across all study sites. Salinity was measured and did not vary much over the sites. Wave exposure was not measured directly, however, the coarse fraction of the sediment grain size was used as a surrogate in initial investigations. The results confirmed those of Hewitt et al. (2005) that wave exposure did not confound identification of effects of pollution on community composition. At each site, a 2 × 5 grid of ten plots (each 10 m × 10 m) was marked out and one replicate collected from each plot, yielding 750 samples overall. 2.2. Physico-chemical variables At each site, one 2 cm diameter core extending 2 cm deep into the sediment was collected from each of the 10 plots in the grid yielding 10 replicates for each site. The replicates were composited into a single sample and the sediment was analysed for grain size, organic matter (loss on ignition, LOI), nutrients (total nitrogen, TN; total phosphorous, TP), heavy metals (lead, Pb; zinc, Zn; copper, Cu) and chlorophyll-α (chl-α). Sediment grain size was determined by wet sieving and calculation of dry weight percentage fractions. Grain size fractions were gravel (N2 mm), very coarse sand (b2 mm and N1 mm), coarse sand (b1 mm and N500 μm), medium sand (b500 μm and N 250 μm), fine sand (b 250 μm and N125 μm), very fine sand (b 125 μm and N63 μm) and mud (b63 μm). Loss on ignition was measured by dry sediment loss after combustion at 550 °C (American Public Health Association 21st Edn 2540 D+ E (Mod); APHA, 2005). Total nitrogen was measured through the analysis and addition of total Kjeldahl nitrogen (TKN) and total oxidized nitrogen (TON). TKN was determined, after digestion with sulphuric acid using copper sulphate as a catalyst, using a Flow Injection Analyser (APHA 21st Edn 4500N C; APHA, 2005). TON was extracted from the sediment sample using potassium chloride and also analysed by a Flow Injection Analyser (APHA 21st Edn 4500N C; APHA, 2005). Sediment samples for total phosphorous and metals were dried at 30 °C then digested using a combination of nitric and hydrochloric acids, with heating to 95 °C for 30 min. Analysis for individual metal

178

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

and TP concentrations was by inductively coupled plasma mass spectrometry (USEPA 200.2 Digestion/ICP-MS; USEPA, 1994). Chlorophyll-α was extracted from a measured sediment sample weight in 90% acetone. The slurry was clarified using a centrifuge and read on a spectrophotometer at 665 and 750 nm, before and after acidification (Biggs and Kilroy, 2000). 2.3. Infauna Infaunal samples were collected to quantify benthic community structure at each site. One 13 cm diameter core extending 15 cm into the sediment was taken from each of the 10 plots in the grid yielding 10 replicates for each site. Core samples were sieved on a 0.5 mm mesh and macrofauna retained on the sieves preserved with ethanol. Only three of the ten replicates were sorted and identified to the lowest taxonomic resolution due to budgetary constraints. For models based on the relationship between community composition and pollution rank in Auckland, Anderson et al. (2002) found that models constructed from data using three cores at each site were less precise than models utilizing a greater number of replicates, however, they were not biased in any way. Macrofauna count data were used to calculate four univariate biological measures; total abundance (total number of individual organisms in a sample), number of taxa, Pielou's evenness (a measure of equitability, Pielou, 1966) and Shannon– Wiener diversity (a measure of community diversity, Shannon, 1948). 2.4. Statistical analysis 2.4.1. Step one: determine the environmental gradients Preliminary analysis using multivariate linear regression to select variables that explain the maximum variation in the community composition data cloud (McArdle and Anderson, 2001) was performed using Distance based Linear Modelling (DistLM Primer E; Clark and Gorley, 2006). DistLM was analysed using the square-root transformed Bray–Curtis similarities with a backward selection procedure (AIC selection criteria). This analysis identified percentage mud as the key grain size variable explaining benthic community composition which was subsequently used as the stress gradient for sedimentation. A number of New Zealand studies have demonstrated that mud context can be used as an indicator of stress related to terrigenous sedimentation (Thrush et al., 2003, 2005). However, for nutrients and contaminants a number of correlated variables had been measured and were determined by DistLM to be important. Principal component analyses (PCA; Jolliffe, 2002) are often used to reduce a set of correlated variables to orthogonal (uncorrelated) axes (e.g. Everitt, 2005). Here we used PCA to derive a single variable (the 1st principal component axis) that would characterize an overall gradient corresponding to increases in the concentrations of all nutrients, and this process was repeated for metals. Total phosphorus, total nitrogen and chl-α were chosen as key parameters for nutrients and copper, zinc and lead were identified as being important for contamination. PCAs were performed on the basis of square root transformed nutrient concentrations (TP, TN, chl-α) and log transformed metal concentrations (Cu, Zn, Pb) using the Primer v6 computer program (Clark and Gorley, 2006). The PC1 axis for nutrients (PC1nutrients) explained 91% of the variance in nutrient information and the PC1 axis for contaminants (PC1contamination) explained 85.5% of the variance in heavy metals (Table 1). In order to determine how these three environmental gradients might interact to affect benthic community composition, DistLM was used to partition the variance explained by the factors. DistLM was run seven times to obtain the percentage explained in benthic community composition (R2) by each of the three factors (sediment, nutrients and metals), then each pairwise combination and finally all three groups. DistLM allows the percentage explained by different factors to be partitioned out by retaining different variables in the models. These results can then be analysed by variance partition methods (Anderson

Table 1 Model details associated with the three canonical analyses of principal coordinates (CAP) models developed for Tauranga Harbour. Model

Sedimentation

Nutrients

Contamination

Key anthropogenic stressors (DistLM) Variables used in PCA Transformation Stress gradient (PC1 axis) Variation explained by PC1 axis

% mud content

TP, chl-α

Cu, Pb

– – % mud content –

TP, chl-α, TN Square root PC1nutrients 91%

Cu, Pb, Zn Log PC1contamination 85.5%

and Gribble, 1998; Borcard et al., 1992) to determine relative percentages explained by different components. 2.4.2. Step two: relate anthropogenic stress gradients to biotic assemblages In order to determine whether there was a significant relationship between the soft sediment faunal communities of the Tauranga region and each anthropogenic stressor (sedimentation, nutrients and contamination), the data was analysed using constrained ordination. A canonical analysis of principal coordinates, or CAP (Anderson and Robinson, 2003; Anderson and Willis, 2003) was applied, which allows a constrained ordination to be done on the basis of any dissimilarity or distance measure of choice (such as the Bray–Curtis measure; Bray and Curtis, 1957) and determines the axes that best discriminates an environmental gradient. A separate CAP analysis was done for each anthropogenic stressor. All CAP analyses were performed using Primer 6 (version 6.1.13) and PermAnova (version 1.0.3) based on Bray–Curtis dissimilarity. The model output was used to place benthic assemblages at each site along an anthropogenic stress gradient from healthy to impacted sites. A non-hierarchical clustering method (k-means; MacQueen, 1967) of the key variables for each model (percentage mud or PC1 axis) was used to identify possible groupings (ecological health categories; EHC) for each anthropogenic stress gradient (sedimentation, nutrients and contamination). Because k-means does not identify the best solution, the optimal number of sites occurring in each group was determined using the Calinski–Harabasz stopping criterion (Calinski and Harabasz, 1974). It is important to note that reference to healthy and impacted is relative to the sampling sites. For example, a site ranked as being highly impacted by sedimentation would have a high percentage of mud relative to other sites sampled in Tauranga Harbour, however, on a global scale it might not necessarily be considered to be highly impacted by sediments. 2.4.3. Step three: identify key species driving relationships between biotic assemblages and anthropogenic stress gradients DistLM modelling was used to determine key sensitive and pollution-tolerant species that may be driving the assemblage differences across the anthropogenic stress gradients for sedimentation, nutrients and contamination. All DistLM modelling was performed using Primer 6 (version 6.1.13) and PermAnova (version 1.0.3). 2.4.4. Step four: examine the effect of interactions between stressors We also wanted to explore the nature of important interactions between key anthropogenic stressors. We therefore modelled the dissimilarity in benthic community composition between sites with increasing nutrient levels for three distinct sediment levels (low, medium and high mud). This was based on ordination distance from nominated “healthy” sites where healthy sites were selected based on the CAP results for sediment, nutrient and contaminant CAP models. 3. Results 3.1. Physico-chemical variables Sediments within Tauranga Harbour were predominantly sandy (51–100% sand), with the exception of a site in the inner reaches of Te

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

Puna Estuary, which was primarily mud (76% mud; Fig. 1). Organic content of sediments in the harbour generally ranged from 0.9 to 4.5% LOI with the exception of one site. At 10% LOI, the organic content of sediments from the inner Te Puna Estuary site, the muddiest site sampled, was much higher than measured in the rest of the harbour. Inner harbour areas contained more mud and had higher organic content than outer harbour sites.

179

Nutrient concentrations also tended to decline with distance from the inner harbour region and associated rivers (Fig. 1). TN in sediments ranged from 140 to 1000 mg/kg and TP from 51 to 340 mg/kg. The inner Te Puna Estuary site showed comparatively high nutrient levels, with nitrogen and phosphorous concentrations of 1900 and 580 mg/kg, respectively. Sediment chl-α concentrations generally ranged from 1100 to 16,000 μg/kg, with particularly low concentrations (210 μg/kg) at

Fig. 1. Physicochemical variables sampled from sediments across 75 sites within Tauranga Harbour. a) Total nitrogen, b) total phosphorus, c) grain-size (as a percentage of gravel, sand and mud), d) zinc, e) lead, and f) copper. Major rivers and streams entering the harbour are shown with colours in (a) depicting modelled nitrogen loading (in ppb) for each segment (estimated from FENZ using CLUES; Woods et al., 2006; Leathwick et al., 2010). ANZECC (2000b) Low Interim Sediment Quality Guidelines (ISQG) for each metal are displayed in the legend.

180

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

Table 2 Relative percentage of variation in benthic community composition explained by different anthropogenic stressors using variance partitioning methods. Anthropogenic stressors

Relative % variation explained

Sedimentation Nutrients Contamination Sedimentation ∗ nutrients Sedimentation ∗ contamination Nutrients ∗ contamination Sedimentation ∗ nutrients ∗ contamination

4.9 2.7 7.5 6.1 0.7 0.8 1.7

Sedimentation (% mud), nutrients (total nitrogen, total phosphorous, chlorophyll-α), contamination (copper, lead, zinc).

an outer harbour site (Blue Gum Bay). There was no obvious relationship between chl-α and nutrient concentrations. Heavy metal concentrations in the harbour tended to be higher in inner areas compared with outer sites but all were below Australian and New Zealand Environment and Conservation Council (ANZECC, 2000b) Interim Sediment Quality Guidelines, which provide thresholds for possible biological effects (ISQG-Low; Cu 65, Pb 50, Zn 200 mg/kg; Fig. 1). However, concentrations were not below levels previously demonstrated to have effects on infauna (Hewitt et al., 2009). The inner Te Puna site had the highest Cu and Pb concentrations (6.1 and 13mg/kg, respectively), and the second highest Zn concentration (46 mg/kg) after the nearby outer Te Puna Estuary site (55 mg/kg).

3.2. Statistical results 3.2.1. Key anthropogenic stressors and interaction terms Variance partitioning (Anderson and Gribble, 1998; Borcard et al., 1992) identified high mud and contaminant concentration as the individual factors that explained the highest proportions of observed variation in community composition (Table 2). Some of the unexplained variation in benthic community composition is undoubtedly attributable to unmeasured abiotic factors, however, the amount of variation unexplained is not high for ordination analyses (Hewitt et al., 2002; ter Braak and Verdonschot, 1995). The sedimentation ∗ contaminant interaction could be a result of binding of heavy metals to fine silt/clays however the relative percentage explained by the interaction term was low (0.7%, Table 2). A relatively high percentage of variance was however explained by the sedimentation ∗ nutrients interaction. Comparing dissimilarity in benthic community composition between sites with increasing nutrient levels for low, medium and high mud revealed different responses to increasing nutrient levels (Fig. 2). At sites with high mud content, increasing nutrients resulted in high community variability and increasing dissimilarity to those sites that were expected to be less impacted. In low mud content sites, variability was low and there was a threshold effect of increasing nutrients. Initially there was no response but towards the end of the increasing nutrient gradient a 20% decrease in dissimilarity occurred indicating further homogenization of community composition.

Fig. 2. Dissimilarity in benthic community composition between sites based on ordination distance from nominated “healthy” sites with increasing nutrient values, at three distinct sediment levels. Light grey lines and shaded grey area represent the 95% confidence intervals.

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

181

3.2.2. Canonical analysis of principal coordinates models Results indicated that the sites identified as most impacted, for all three CAP models (sedimentation, nutrients and contamination), were located in the upper reaches of estuaries in some of the least exposed locations. In addition, the sensitivities of organisms characterizing sites that have different sediment textures, as well as contaminant and nutrient loadings, were found to vary. 3.2.3. Sedimentation canonical analysis of principal coordinates model A strong gradient of community change was observed in response to mud content of the sediment (R2 = 0.7683) suggesting that the multivariate model can be used to determine potential effects of changes in sediment mud content (Fig. 3a). The sedimentation gradient was based on the percentage mud content in the sediment, with the muddiest sites (48–49% mud) ranked as ‘8’ (impacted) and sandiest sites (b0.1–1.8% mud) as ‘1’ (healthy; Table 3). Over half (58%) of the sites were ranked in ecological health categories (EHC) less than ‘5’, with most (23%) in EHC ‘2’, suggesting fairly healthy communities with regard to sedimentation (Fig. 3). The organic content of the sediment increased with increasing ecological health category (mean 1.3% LOI in EHC ‘1’ increasing to 4.4% in EHC ‘8’), reflecting the tendency of organic material to accumulate in fine sediments. Sites in EHCs ‘7’ and ‘8’ (0.1% of sites), the most impacted ecological health category, were found in inner estuaries (Fig. 4) where deposition of sediments would be expected to be highest. Conversely, sites in EHCs ‘1’ and ‘2’ (32% of sites) tended to be in outer areas of the harbour. Infauna abundance was lowest at the most impacted sites (EHC ‘7’ and ‘8’, means of 85 and 56 per core, respectively) but no clear trend was observed amongst the other categories (means of 110–138 per core; Table 3). In general, species richness was slightly higher at healthier sites (means of 25–30 taxa per site) than more impacted sites (means of 20–24 taxa per site). Evenness showed no clear trend amongst categories and Shannon–Wiener diversity was only slightly lower at the most impacted sites (EHC ‘7’ and ‘8’). The key species differences at healthy versus impacted sites along an increasing gradient of siltation are provided in Table 4. Key species associated with high mud included the polychaete worms Nereididae, Scolecolepides benhami and Heteromastus filiformis and the deposit feeding bivalve Arthritica bifurca. Key species associated with low mud included the worms Scoloplos cylindrifer and Scolelepis sp. and the gastropod Halopyrgus pupoides. 3.2.4. Nutrient canonical analysis of principal coordinates model The nutrient CAP model was based on a constrained ordination of benthic community taxa in relation to the stress gradient (PC1nutrients) generated from the concentrations of TN, TP and chl-α at each site. While a significant community change was observed in response to this nutrient gradient, suggesting that the multivariate model can be used to determine potential effects of changes in nutrient concentrations, its correlation was lower (R2 = 0.5135) than the CAP models for sedimentation and contamination (Fig. 3b). Most of the sites (26%) were ranked in EHC ‘2’, with very few in the most healthy (‘1’) and impacted (‘6’) categories, suggesting a relatively healthy harbour, although slightly compromised, with respect to elevated nutrients (Fig. 3). The level of response to the PCA axis was most closely related to concentrations of nitrogen and phosphorous in the sediment, with increasing TN and TP concentrations as ecological health categories increased (i.e. became more impacted; Table 3). Organic content also increased along the nutrient gradient but chl-α was highest

Fig. 3. Canonical analysis of principal coordinates (CAP) for sedimentation, nutrient and contamination multivariate models based on data from 74 sites in Tauranga Harbour. Grey lines demarcate the ecological health categories for each model with lower numbers indicating a relatively ‘healthy’ community and higher numbers indicating a more ‘impacted’ community.

182

Table 3 Ecological health categories (EHC) for three multivariate models; sedimentation, nutrients and contaminants. n

Gravel

Sand

Mud

LOI

TN

TP

chl-α

Pb

Cu

Zn

N

S

J

H

1 2 3 4 5 6 7 8

7 17 10 9 14 8 7 2

b0.1–4.0 0.1–14.6 0.1–5.9 0.1–7.1 0.1–5.6 b0.1–1.6 0.2–10.2 0.3–1.5

94.2–100 82.8–96.7 62.7–94.0 83.4–91.0 51.6–87.5 73.9–81.3 60.0–68.4 50.7–51.0

b0.1–1.8 2.5–4.9 5.1–7.7 8.9–10.9 12.2–17.5 18.5–25.4 27.9–38.6 47.5–48.9

0.9–1.8 1.0–3.0 1.9–3.8 1.8–3.4 2.5–4.2 3.1–4.5 3.1–4.5 4.2–4.5

– – – – – – – –

– – – – – – – –

– – – – – – – –

– – – – – – – –

– – – – – – – –

– – – – – – – –

42–257 29–263 59–164 68–267 39–268 61–333 38–133 35–78

17–39 16–34 14–36 23–37 14–35 10–39 14–31 19–21

0.6–0.8 0.5–0.9 0.5–0.8 0.6–0.8 0.4–0.8 0.08–0.8 0.4–0.8 0.6–0.8

1.9–2.7 1.5–3.0 1.4–2.7 1.8–2.9 1.2–2.8 0.2–2.7 1.1–2.4 1.7–2.4

Nutrients

1 2 3 4 5 6

7 19 16 16 11 5

– – – – – –

– – – – – –

– – – – – –

0.9–1.9 1.4–3.4 1.2–3.5 2.4–4.3 2.4–4.2 4.0–4.5

140–220 250–460 280–520 390–690 540–700 140–220

51–110 91–160 120–180 130–220 180–260 51–110

210–4000 1200–11,000 1200–15,000 3600–16,000 2800–10,000 1100–9600

– – – – – –

– – – – – –

– – – – – –

49–170 29–257 46–268 61–333 35–239 78–133

17–39 17–33 14–33 10–39 14–34 13–31

0.6–0.8 0.5–0.9 0.5–0.8 0.08–0.8 0.4–0.8 0.4–0.7

1.9–2.7 1.5–3.0 1.4–2.7 0.2–2.9 1.1–2.4 1.1–2.4

Contamination

1 2 3 4 5

9 18 18 19 10

– – – – –

– – – – –

1.3–5.6 b0.1–13.0 2.6–34.2 3.6–48.9 12.5–47.5

– – – – –

– – – – –

– – – – –

– – – – –

b1.0–1.4 b1.0–1.9 1.2–3.1 1.4–5.1 3.0–5.6

b1.0 b1.0 b1.0–1.6 b1.0–2.2 1.3–3.0

b5.1–8.0 7.5–12.0 12.0–20.0 18.0–28.0 26.0–55.0

29–257 42–241 56–333 35–198 67–268

17–32 16–39 10–35 14–39 13–34

0.6–0.8 0.5–0.8 0.08–0.9 0.4–0.8 0.4–0.7

1.9–2.6 1.5–2.7 0.2–3.0 1.2–2.9 1.1–2.4

Notes: Ranges for key variables associated with each of the models are shown for each EHC. Inner Te Puna site excluded from CAP analysis because it was an outlier (outside the range of variation observed at the other sites). n = number of sites in each ecological health category; LOI = loss on ignition; TN = total nitrogen; TP = total phosphorous; chl-α = chlorophyll-α; Pb = lead, Cu = copper; Zn = zinc; N = total abundance per core; S = total number of taxa per site; J = Pielou's evenness; and H = Shannon–Wiener index. Gravel (N2 mm), sand (b2 mm and N63 μm), mud (b63 μm) and LOI are measured as a percentage; TN, TP, Pb, Cu, Zn are measured in mg/kg; chl-α is measured as μg/kg.

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

EHC Sedimentation

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

at the moderately impacted sites. Sites in EHCs ‘5’ and ‘6’, the most impacted categories, were generally found in estuaries along the inner coast of the harbour, whereas sites ranked lower tended to be situated in the outer harbour (Fig. 5). In general, the univariate biological measures were not as sensitive at detecting differences across the ecological health categories. No clear trend in macrofaunal abundance was apparent with infauna numbers highest at EHC ‘3’ and ‘4’ sites (means of 135 and 141 per core, respectively) and lowest at EHC ‘1’ sites (mean of 84 per core; Table 3). Species richness was similar in the first four categories (means of 24– 27 taxa per site) but slightly lower at EHC ‘5’ and ‘6’ sites (means of 22 and 19 taxa per site, respectively). Evenness showed no clear trend amongst categories and Shannon–Wiener diversity was only slightly lower at the most impacted sites (EHC ‘5’ and ‘6’). The polychaete Scolelepis sp. was associated with low nutrients while key species tolerant of elevated nutrient loadings included the polychaete worms, S. benhami and H. filiformis, and most of the amphipods (Table 4). 3.2.5. Contamination canonical analysis of principal coordinates model The contamination CAP model was based on a constrained ordination of benthic community taxa in relation to the stress gradient

183

(PC1contamination axis) generated from the concentration of heavy metals (Pb, Cu and Zn) at each site. A strong gradient of community change was observed in response to heavy metal concentrations in the sediment (R2 = 0.7075) suggesting that the multivariate model can be used to determine potential effects of changes in metal concentrations, despite the relatively low concentrations covered by the gradient (Fig. 3c). Most of the sites were ranked in EHCs ‘2’ (24%), ‘3’ (24%) and ‘4’ (26%; Fig. 3). All metal concentrations increased with increasing ecological health category (Table 3). The mud content of the sediment also tended to increase as ecological health became more impacted, reflecting the tendency of metals to bind with fine sediments. As with the other CAP models, EHC ‘5’ sites tended to be situated in inner harbour areas and EHC ‘1’ and ‘2’ sites further out (Fig. 6). As with nutrients, the univariate biological measures in general were not as sensitive as the multivariate ordinations at detecting differences across the anthropogenic stress gradient. Infauna numbers were lowest at EHC ‘4’ sites (mean of 91 per core) and highest at EHC ‘3’ and ‘5’ sites (means of 136 and 134 per core, respectively; Table 3). Species richness did not greatly differ across the first four categories (means of 25–28 taxa per site) and was reduced slightly only at the most polluted sites

Fig. 4. Canonical analysis of principal coordinates (CAP) for sedimentation in Tauranga Harbour for 74 sites. Colours indicate the ecological health categories where a green (low) ranking indicates a low effect of sedimentation (‘healthy’) and a red (high) ranking indicates a high effect (‘impacted’). Major rivers and streams entering the harbour are shown in blue.

184

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

Table 4 Key species identified along ecological health gradients for sedimentation, nutrient and contamination models as determined using Distance based Linear Models (DistLM). Model

Association

Species

Faunal group

Feeding mode

Sedimentation

Low mud

Scoloplos cylindrifer Scolelepis sp. Halopyrgus pupoides Scolecolepides benhami Heteromastus filiformis Arthritica bifurca Nereididae Scolelepis sp. Scolecolepides benhami Heteromastus filiformis Amphipoda indeterminata Orbinia papillosa Scolecolepides benhami Heteromastus filiformis Arthritica bifurca Amphipoda indeterminata

Orbinid polychaete Spionid polychaete Gastropod Spionid polychaete Capitellid polychaete Bivalve (deposit feeding) Nereid polychaete Spionid polychaete Spionid polychaete Capitellid polychaete Amphipod Orbinid polychaete Spionid polychaete Capitellid polychaete Bivalve (deposit feeding) Amphipod

D (surface/subsurface) D Microalgal and detrital grazer D (surface deposit) D (sub-surface deposit) D P D (surface) D (surface deposit) D (sub-surface deposit) D, P, G D D (surface deposit) D (sub-surface deposit) D D, P, G

High mud

Nutrients

Negative Positive

Contamination

Negative Positive

Species that respond negatively to increasing nutrients/contaminants are sensitive to elevated nutrients/contaminants, while species that respond positively to increasing nutrients/contaminants are more tolerant to that stressor and can be found at sites with high nutrient/contaminant loadings. Abbreviations for feeding mode: D = deposit feeder; P = predator/scavenger; S = suspension feeder; and G = grazer.

(EHC ‘5’; mean 21 taxa per site). Evenness showed no clear trend amongst categories and Shannon–Wiener diversity was only slightly lower at the most impacted sites (EHC ‘5’).

Key species associated with low contaminant loadings included the polychaete worm Orbinia papillosa while species tolerant of elevated contaminant loadings included the polychaete worms S. benhami,

Fig. 5. Canonical analysis of principal coordinates (CAP) for nutrients in Tauranga Harbour for 74 sites. Colours indicate the ecological health categories where a green (low) ranking indicates a low effect of nutrients (‘healthy’) and a red (high) ranking indicates a high effect (‘impacted’).

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

H. filiformis, amphipods and the deposit feeding bivalve A. bifurca (Table 4). 4. Discussion Distance based Linear Modelling identified sediments, nutrients and heavy metals as key stressors affecting the ‘health’ of macrofaunal communities in Tauranga Harbour. Three multivariate models were developed based on the variability in community composition using CAP analyses. The ecological assemblages generally reflected gradients of stress or pollution very well. However, the CAP models for sedimentation and contamination performed better than for nutrients. The multivariate models were more sensitive to changing environmental health than simple univariate measures (based on abundance, species richness, Pielou's evenness and Shannon–Wiener diversity). Species richness and Shannon–Wiener diversity only detected changes at the most impacted sites and no clear patterns in evenness were observed. No change in total macrofaunal abundance was evident across the ecological health gradients for nutrients and contamination and differences were only apparent at the most disturbed sites for sedimentation. This trend has also been reported in the literature where univariate measures found significant differences between the most and least disturbed sites, but not between smaller relative differences (Attayde and

185

Bozelli, 1998). The use of indicator species and diversity metrics are often based on contaminant toxicity tests and more generally on ecological theories of disturbance. Biological diversity indices, such as the Shannon–Wiener index, do not differentiate amongst different types of taxa which therefore limits their ability to detect changes in composition across different habitat types and communities (Gappa et al., 1990; Hewitt et al., 2005). Further a similar diversity value can be obtained from communities with very different species assemblages, therefore interpretation of what these indices mean to ecological functioning is also difficult (Clarke, 1993; Dufrene and Legendre, 1997). The increased sensitivity of the constrained ordination models based on community composition relative to univariate indicators means they are likely to be effective tools in detecting changes in ‘health’ and enable long-term degradative change from multiple disturbances to be assessed. This approach can be used as a management or monitoring tool where sites are repeatedly sampled over time and tracked to determine whether the communities are moving towards a more healthy or unhealthy state. Further the ordination models identify changes in benthic community composition along gradients of stress as well as categorising sensitive species. Changes in species composition can be linked with changes in functional components which provide information about how communities respond to environmental stress and in some contexts of their resilience (de Juan et al., 2015; Lavorel and

Fig. 6. Canonical analysis of principal coordinates (CAP) for contamination in Tauranga Harbour for 74 sites. Colours indicate the ecological health categories where a green (low) ranking indicates a low effect of contamination (‘healthy’) and a red (high) ranking indicates a high effect (‘impacted’).

186

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187

Garnier, 2002). The scientific challenge however is to scale up studies of functional components and “ecosystem health” to assess their relevance for monitoring diverse marine ecosystems as well as their dynamics and response to different sources of stress (de Juan et al., 2015). Hewitt et al. (2005) tested the utility of three multivariate constrained ordination techniques to assess changes in communities occurring along an anthropogenic disturbance gradient of stormwater pollution in two different habitats (estuaries and harbours) within a region of New Zealand. They concluded that the approach should be able to be utilized in other marine (or freshwater or terrestrial) systems where commensurable regional-scale multivariate databases exist. We have successfully applied this approach to another region and were able to develop models for a harbour system experiencing multiple sources of anthropogenic disturbances. Specifically, we extended this analysis by also developing models of macrofaunal species occurrence with respect to sediment mud content and nutrients as well as contaminant loadings from stormwater. This method may therefore provide a sensitive management tool for assessing ecological health to multiple stressors at regional scales. Our study also found that responses to two of the three stressors interacted. Specifically, variance partitioning methods indicated that there were intersecting responses to sedimentation and nutrients. The intersection term for these two stressors explained a higher percentage of the variance in benthic community composition than the other intersecting terms, which were generally small. Interestingly separating out benthic responses to nutrients into those in low, medium and high mud content sites we observed an interactive effect. We found variable but increasing dissimilarity with increasing nutrients for muddy sites and low variability with a threshold followed by a slight drop in dissimilarity for sandy sites. This disparity could be explained either by differences in nutrient dynamics or by community responses. For example differences in nutrient dynamics such as nutrient pore water transport processes in muddy versus sandy sediment dynamics may explain the observed responses. Furthermore, more species prefer sandy sites, allowing more different types of communities to form with high dissimilarity between them. A stressor could then decrease community similarities (at least over a relatively low stress range) such as the nutrient concentrations observed here. Our study provides evidence for interaction effects on estuarine macrobenthic community composition with multiple stressors, particularly between sedimentation and nutrients. Specifically we show shifts in the spatial heterogeneity of community composition, a phenomenon that would be very difficult to resolve with univariate analysis of community data. Interaction effects are likely to occur with nutrient enrichment of many coastal regions (Schindler, 2006). Further, these interactions are likely to increase as climate variability and extreme climatic events strongly influence the delivery of freshwater and associated nutrients as well as sediment loads to the coastal zone (Scavia et al., 2002), suggesting the need for further research, preferably linking infaunal responses to nutrient dynamics within the sediment and water interface. Within this study species displayed clear differences in abundance as a function of sediment type, nutrient loading or contaminant levels. Further work on the response of key species to stressors, detailed in Ellis et al. (under review), showed that most shellfish species populations were either sensitive to elevated mud, nutrient loading or contaminants, or sensitive to these stressors beyond a critical point. The response of polychaete worms to these stressors varied by species. 5. Conclusions Increasingly coastal and estuarine ecosystems are subject to multiple sources of anthropogenic stress. While many studies have documented the individual effects of these stressors on species and ecosystems, research on cumulative and interactive impacts of multiple stressors is less frequent (Crain et al., 2008). The rate of habitat change in estuarine and coastal ecosystems is likely to increase as a

consequence of global warming and increasing anthropogenic pressure on coastal areas. Within this study we successfully developed multivariate models to assess changes in communities occurring along anthropogenic disturbance gradients for sediments, nutrients and contaminants, building on models developed by Hewitt et al. (2005). Comparison with diversity indices suggested that the multivariate approach was more successful in defining changes and thus ecological health. The multivariate approach also allows monitoring of community composition to determine whether sites are improving or degrading over time for multiple stressors. Acknowledgements We thank the Ministry of Business, Innovation and Employment (MBIE) for funding this work (contract MAUX0907) and the Bay of Plenty Regional Council for providing co-funding for the data analysis. We greatly appreciate the help we received from various parties while undertaking the ecological survey; the University of Waikato, in particular Prof Chris Battershill, the Bay of Plenty Regional Council, Manaaki Te Awanui, Waka Digital and the many volunteers who assisted with field work. [SS] References Anderson, M.J., 2008. Animal–sediment relationships re-visited: characterising species' distributions along an environmental gradient using canonical analysis and quantile regression splines. J. Exp. Mar. Biol. Ecol. 366, 16–27. Anderson, M.J., Gribble, N.A., 1998. Partitioning the variation among spatial, temporal and environmental components in a multivariate dataset. Aust. J. Ecol. 23, 158–167. Anderson, M.J., Robinson, J., 2003. Generalized discriminant analysis based on distances. Aust. N. Z. J. Stat. 45, 301–318. Anderson, M.J., Willis, T.J., 2003. Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology. Ecology 84, 511–525. Anderson, M., Hewitt, J., Thrush, S., 2002. The development of criteria for assessing community health of urban intertidal flats. NIWA Report Prepared for Auckland Regional Council Technical Report No. 184, p. 54. ANZECC, 2000a. Australian and New Zealand Guidelines for Fresh and Marine Water Quality — Volume 2: Aquatic Ecosystems — Rationale and Background Information. Australian and New Zealand Environment and Conservation Council and Agriculture and Resource Management Council of Australia and New Zealand. ANZECC, 2000b. Australian and New Zealand Guidelines for Fresh and Marine Water Quality — Volume 1: The Guidelines. Australian and New Zealand Environment and Conservation Council and Agriculture and Resource Management Council of Australia and New Zealand. APHA, 2005. Standard Methods for the Examination of Water and Wastewater. 21st ed. American Public Health Association, Washington, DC. Attayde, J.L., Bozelli, R.L., 1998. Assessing the indicator properties of zooplankton assemblages to disturbance gradients by canonical correspondence analysis. Can. J. Fish. Aquat. Sci. 55, 1789–1797. Biggs, B., Kilroy, C., 2000. Periphyton Monitoring Manual. Published by NIWA for Ministry for the Environment 0-478-09099-4. Borcard, D.P., Legendre, P., Drapeau, P., 1992. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055. Borja, A., Franco, J., Perez, V., 2000. A marine biotic index to establish the ecological quality of soft-bottom benthos within European estuarine and coastal environments. Mar. Pollut. Bull. 40 (12), 1100–1114. Bray, J.R., Curtis, J.T., 1957. An ordination of the upland forest communities of southern Wisconsin. Ecol. Monogr. 27, 325–349. Bricelj, V.M., Malouf, R.E., 1984. Influence of algal and suspended sediment concentrations on the feeding physiology of the hard clam Mercenaria mercenaria. Mar. Biol. 84, 155–165. Calinski, T., Harabasz, J., 1974. A dendrite method for cluster analysis. Comput. Stat. 3, 1–27. Clark, K.R., Gorley, R.N., 2006. Primer v6: User Manual/Tutorial. PRIMER-E Ltd, Plymouth, UK. Clarke, K.R., 1993. Non-parametric multivariate analyses of changes in community structure. Aust. J. Ecol. 18, 117–143. Cloern, J.E., 2001. Our evolving conceptual model of the coastal eutrophication problem. Mar. Ecol. Prog. Ser. 210, 223–253. Crain, C.M., Kroeker, K., Halpern, B.S., 2008. Interactive and cumulative effects of multiple human stressors in marine systems. Ecol. Lett. 11 (12), 1304–1315. de Juan, S., Hewitt, J., Thrush, S., Freeman, D., 2015. Standardizing the assessment of functional integrity in benthic ecosystems. J. Sea Res. 98, 33–41. Dufrene, M., Legendre, P., 1997. Species assemblages and indicator species: the need for a flexible asymmetrical approach. Ecol. Monogr. 67, 345–366. Dunn, C., 1994. Gaps in GAP. Plant Sci. Bull. 40, 120–121. Ellis, J.I., Schneider, D.C., Thrush, S.F., 2000. Detecting anthropogenic disturbance in an environment with multiple gradients of physical disturbance, Manukau Harbour, New Zealand. Hydrobiologia 440, 379–391.

J.I. Ellis et al. / Journal of Experimental Marine Biology and Ecology 473 (2015) 176–187 Ellis, J.I., Cummings, V., Hewitt, J., Thrush, S., Norkko, A., 2002. Determining effects of suspended sediment on condition of a suspension feeding bivalve (Atrina zelandica): results of a survey, a laboratory experiment and a field transplant experiment. J. Exp. Mar. Biol. Ecol. 267 (2), 147–174. Ellis, J.I., Clark, D., Taiapa, C., Sinner, J., Patterson, M., Hewitt, J., 2015w. Species distribution models along multiple stressor gradients. J. Exp. Mar. Biol. Ecol. (under review). Everitt, B., 2005. An R and S-PLUS Companion to Multivariate Analysis. Springer, London. Fleeger, J.W., Carman, K.R., Nisbet, R.M., 2003. Indirect effects of contaminants in aquatic ecosystems. Sci. Total Environ. 317 (1-3), 207–233. Gagnaire, B., Thomas-Guyon, H., Renault, T., 2004. In vitro effects of cadmium and mercury on Pacific oyster, Crassostrea gigas (Thunberg), haemocytes. Fish Shellfish Immunol. 16 (4), 501–512. Gappa, J.J.L., Tablado, A., Magaldi, N.H., 1990. Influence of sewage pollution on a rocky intertidal community dominated by the mytilid Brachidontes rodriguezii. Mar. Ecol. Prog. Ser. 63, 163–175. Gray, J.S., 2000. The measurement of marine species diversity, with an application to the benthic fauna of the Norwegian continental shelf. J. Exp. Mar. Biol. Ecol. 250, 23–49. Gray, J.S., Mirza, F.B., 1979. A possible method for the detection of pollution-induced disturbance on marine benthic communities. Mar. Pollut. Bull. 10, 142–146. Halpern, B.S., Selkoe, K.A., Micheli, F., Kappel, C.V., 2007. Evaluating and ranking the vulnerability of global marine ecosystems to anthropogenic threats. Conserv. Biol. 21 (5), 1301–1315. Hewitt, J.E., Legendre, P., Thrush, S.F., Cummings, V.J., Norkko, A., 2002. Integrating results from different scales: a multi-resolution study of the relationships between Atrina zelandica and macrofauna along a physical gradient. Mar. Ecol. Prog. Ser. 239, 115–128. Hewitt, J.E., Anderson, M.J., Thrush, S.F., 2005. Assessing and monitoring ecological community health in marine systems. Ecol. Appl. 15 (3), 942–953. Hewitt, J.E., Anderson, M.J., Hickey, C.W., Kelly, S., Thrush, S.F., 2009. Enhancing the ecological significance of sediment contamination guidelines through integration with community analysis. Environ. Sci. Technol. 43, 2118–2123. Iglesias, J.I.P., Urrutia, M.B., Navarro, E., Alvarez-Jorna, P., Larretxea, X., Bougrier, S., Heral, M., 1996. Variability of feeding processes in the cockle Cerastoderma edule (L.) in response to changes in seston concentration and composition. J. Exp. Mar. Biol. Ecol. 197, 121–143. Inglis, G., Gust, N., Fitridge, I., Floerl, O., Woods, C., Kospartov, M., Hayden, B., Fenwick, G., 2008. Port of Tauranga: second baseline survey for non-indigenous marine species (research project ZBS2000/04). Prepared for Biosecurity New Zealand by NIWA. Biosecurity Technical Paper No: 2008/0, Wellington. Jackson, J.B.C., Kirby, M.X., Berger, W.H., Bjorndal, K.A., Botsford, L.W., Bourque, B.J., Bradbury, R.H., Cooke, R., Erlandson, J., Estes, J.A., et al., 2001. Historical overfishing and the recent collapse of coastal ecosystems. Science 293 (5530), 629–637. Jolliffe, I.T., 2002. Principal component analysis. 2nd ed. Springer Series in Statistics XXIX. Springer, NY 978-0-387-95442-4 (487 pp.). Kappel, C.V., 2005. Losing pieces of the puzzle: threats to marine, estuarine, and diadromous species. Front. Ecol. Environ. 3, 275–282. Labrune, C., Romero-Ramirez, A., Amouroux, J.M., Duchene, J.C., Desmalades, M., Escoubeyrou, K., Buscail, R., Gremare, A., 2012. Comparison of ecological quality indices based on benthic macrofauna and sediment profile images: a case study along an organic enrichment gradient off the Rhone River. Ecol. Indic. 12, 133–142. Lavorel, S., Garnier, E., 2002. Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail. Funct. Ecol. 16, 545–556. Leathwick, J.R., West, D., Gerbeauz, P., Kelly, D., Robertson, H., Brown, D., Chadderton, W.L., Ausseil, A.-G., 2010. Freshwater Ecosystems of New Zealand (FENZ) Geodatabase: Version one — August 2010 User Guide. Department of Conservation. MacQueen, J., 1967. Some methods for classification and analysis of multivariate observations. In: Le Cam, L.M., Neyman, J. (Eds.), Proceedings of the 6th Berkley Symposium of Mathematical Statistics and Probability vol. 1. University of California Press, Berkley, pp. 281–297. Margalef, R., 1958. Information theory in ecology. Gen. Syst. 3, 36–71. McArdle, B.H., Anderson, M.J., 2001. Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82, 290–297. McGlathery, K.J., Sundback, K., Anderson, I.C., 2007. Eutrophication in shallow coastal bays and lagoons: the role of plants in the coastal filter. Mar. Ecol. Prog. Ser. 348, 1–18. Morse, P.M., Robinson, W.E., Wehling, W.E., 1982. Effects of sublethal concentrations of the drilling mud components attapulgite and Q-broxin on the structure and function of the fill of the scallop Placopecten magellanicus (Gremlin). In: Vernberg, W.B., Calbrese, A., Thurnberg, F.P., Vernberg, F.J. (Eds.), Physiological Mechanisms of Marine Pollution Toxicity. Academic Press, pp. 235–259.

187

Navarro, J.M., Widdows, J., 1997. Feeding physiology of Cerastoderma edule in response to a wide range of seston concentrations. Mar. Ecol. Prog. Ser. 152 (1-3), 175–186. Nicholson, S., 1999. Cytological and physiological biomarker responses from green mussels, Perna viridis (L.) transplanted to contaminated sites in Hong Kong coastal waters. Mar. Pollut. Bull. 39 (1-12), 261–268. Olsgard, F., Gray, J.S., 1995. A comprehensive analysis of the effects of offshore oil and gas exploration and production on the benthic communities of the Norwegian continental shelf. Mar. Ecol. Prog. Ser. 122, 277–306. Peters, E.C., Gassman, N.J., Firman, J.C., Richmonds, R.H., Power, E.A., 1997. Ecotoxicology of tropical marine ecosystems. Environ. Toxicol. Chem. 16 (1), 12–40. Pielou, E.C., 1966. The measurement of diversity in different types of biological collections. J. Theor. Biol. 13, 131–144. Pohle, G., Frost, B., Findlay, R., 2001. Assessment of regional benthic impact of salmon mariculture within the Letang Inlet, Bay of Fundy. ICES J. Mar. Sci. 58, 417–426. Radford, J.L., Hutchinson, A.E., Burandt, M., Raftos, D.A., 2000. Effects of metal-based environmental pollutants on tunicate hemocytes. J. Invertebr. Pathol. 76 (4), 242–248. Robinson, W.E., Wehling, W.E., Morse, M.P., 1984. The effect of suspended clay on feeding and digestive efficiency of the surf clam, Spisula solidissima (Dillwyn). J. Exp. Mar. Biol. Ecol. 74, 1–12. Rosenberg, R., Blomqvist, M., Nilsson, H.C., Cederwall, H., Dimming, A., 2004. Marine quality assessment by use of benthic species-abundance distributions: a proposed new protocol within the European Union Water Framework Directive. Mar. Pollut. Bull. 49, 728–739. Saiz-Salinas, J.I., Urkiaga-Alberdi, J., 1999. Faunal responses to turbidity in a man-modified bay (Bilbao, Spain). Mar. Environ. Res. 47, 331–347. Sala, O.E., Chapin, F.S., Armesto, J.J., Berlow, E., Bloomfield, J., Dirzo, R., et al., 2000. Biodiversity — global biodiversity scenarios for the year 2100. Science 287, 1770–1774. Sanders, H.L., 1968. Marine benthic diversity: a comparative study. Am. Nat. 102, 243–282. Scavia, D., Field, J.C., Boesch, D.F., Buddemeier, R.W., Burkett, V., Cayan, D.R., Fogarty, M., Harwell, M.A., Howarth, R.W., Mason, C., Reed, D.J., Royer, T.C., Sallenger, A.H., Titus, J.G., 2002. Climate change impacts on U.S. coastal and marine ecosystems. Estuar. Coasts 25, 149–164. Schindler, D.W., 2006. Recent advances in the understanding and management of eutrophication. Limnol. Oceanogr. 51, 356–363. Shannon, C.E., 1948. A mathematical theory of communication. Bell Syst. Tech. J. 27, 379–423. Smith, C.R., Kukert, H., 1996. Macrobenthic community structure, secondary production, and rates of bioturbation and sedimentation at the Kane'ohe Bay lagoon floor. Pac. Sci. 50 (2), 211–229. Stevens, P.M., 1987. Response of excised gill tissue from the New Zealand scallop Pecten novaezelandiae to suspended silt. N. Z. J. Mar. Freshw. Res. 21, 605–614. ter Braak, C.J.F., Verdonschot, P.F.M., 1995. Canonical correspondence analysis and related multivariate methods in aquatic ecology. Aquat. Sci. 57, 153–187. Thrush, S.F., Hewitt, J.E., Norkko, A., Nicholls, P.E., Funnell, G.A., Ellis, J.I., 2003. Habitat change in estuaries: predicting broad-scale responses of intertidal macrofauna to sediment mud content. Mar. Ecol. Prog. Ser. 263, 101–112. Thrush, S.F., Lundquist, C.J., Hewitt, J.E., 2005. Spatial and temporal scales of disturbance to the seafloor: a generalized framework for active habitat management. In: Barnes, P.W., Thomas, J.P. (Eds.), Benthic Habitats and the Effects of Fishing. American Fisheries Society, Symposium Series, Bethesda, Maryland, pp. 639–649. Thrush, S., Townsend, M., Hewitt, J., Davie, K., Lohrer, D., Lundquist, C., Cartner, K., 2013. The many uses and values of estuarine ecosystems. In: Dymond, J.R. (Ed.), Ecosystem Services in New Zealand — Conditions and Trends. Manaaki Whenua Press, Lincoln, New Zealand. USEPA, 1994. Methods for the Determination of Metals in Environmental Samples: EPA 200 Series, Supplement I, EPA-600/R-94/111, Revision 2.8, May 1994. U.S. Environmental Protection Agency. Venter, O., Brodeur, N.N., Nemiroff, L., Belland, B., Dolinsek, I., Grant, W., 2006. Threats to endangered species in Canada. Bioscience 56, 903–910. Warwick, R.M., Platt, H.M., Clarke, K.R., Agard, J., Gobin, J., 1990. Analysis of macrobenthic and meiobenthic community structure in relation to pollution and disturbance in Hamilton Harbour, Bermuda. J. Mar. Biol. Ecol. 139, 119–142. Wilcove, D.S., Rothstein, D., Dubow, J., Phillips, A., Losos, E., 1998. Quantifying threats to imperiled species in the United States. Bioscience 48, 607–615. Willows, R.I., 1992. Optimal digestive investment: a model for filter feeders experiencing variable diets. Limnol. Oceanogr. 37, 829–847. Woods, R., Elliott, S., Shankar, U., Bidwell, V., Harris, S., Wheeler, D., Clothier, B., Green, S., Hewitt, A., Gibb, R., Parfitt, R., 2006. The CLUES project: predicting the effects of landuse on water quality — stage II. NIWA Client Report MAF05502.