Impacts of a bivalve mass mortality event on an estuarine food web and bivalve grazing pressure

ABSTRACT

Declines in bivalve populations have been quite common worldwide, often associated with coastal development, pollution and climate change. In addition to the impacts of these chronic stressors, occasional mass mortality events may have severe consequences on ecosystem services and biodiversity. In this study, we examined the impact of a mass mortality event of the clam (Austrovenus stutchburyi) on an estuarine food web and the grazing pressure exerted by the bivalve population. In February 2009, c. 60% of the clam population died in Whangateau Harbour, New Zealand. Population clearance rate calculations suggest that the clam population do not exert significant top-down control on phytoplankton biomass in the estuary, and thus the impact of the mortality event on bivalve grazing pressure was less severe than the reduction in abundance would suggest. A trophic model shows that phytoplankton play a limited role in the estuary food web, which is instead dominated by microphytobenthos and clams. This study highlights the importance of microphytobenthos in shallow estuaries, and the application of the trophic model is a useful tool that can identify key components of the ecosystem and could help inform monitoring programmes.

KEYWORDS:

• Clearance rate
• trophic model
• microphytobenthos
• Austrovenus stutchburyi
• Whangateau Harbour (New Zealand)

Introduction

Suspension feeding bivalves act as key species in estuarine ecosystems, being integral to benthic–pelagic coupling and altering rates of nutrient regeneration (e.g. Smaal & Prins 1993; Thrush et al. 2006). Due to their position in the food web, they can exert considerable top-down control on phytoplankton populations and provide a substantial food source for higher trophic levels (see review by Newell 2004). Previous research has often focussed on the pelagic grazing pressure exerted by suspension feeding bivalves (e.g. Hily 1991; Dame & Prins 1998; Lonsdale et al. 2007; Klimmerer & Thompson 2014) or their influence on nutrient cycling and primary productivity (e.g. Thrush et al. 2006; Jones et al. 2011a). Indeed, in some systems, the loss of suspension feeders has resulted in an increased incidence of phytoplankton blooms, for example, following reductions in eastern oyster (Crassostrea virginica) beds in Chesapeake Bay, USA (Jackson et al. 2001; Kemp et al. 2005) or clams (Mercenaria mercenaria) in Great South Bay, NY, USA (Lonsdale et al. 2007). However, the focus on control of phytoplankton populations by suspension feeding bivalves overlooks the importance of other food sources, such as microphytobenthos that may dominate production in shallow estuaries with high water clarity and extensive intertidal flats (e.g. MacIntyre et al. 1996). In these systems, infaunal suspension feeding bivalves often dominate macrofauna communities (Thrush et al. 2003; Newell 2004; Alfaro 2006), and resuspended microphytobenthos may provide an important food source for these bivalves (e.g. Kang et al. 2003; Safi 2003).

Populations of infaunal suspension feeding bivalves have declined globally in many estuaries, likely associated with increasing pressure from expanding coastal human populations (Piersma et al. 2001; Peterson 2002; Kraeuter et al. 2008). The infaunal suspension feeding bivalve, Austrovenus stutchburyi, is a dominant species on intertidal flats in New Zealand estuaries, often present in large, high-density (c. 1000 ind. m⁻²) beds (Dobbinson et al. 1989; Whitlatch et al. 1997). A. stutchburyi is a key species in these systems, capable of influencing microphytobenthos productivity and nutrient cycling, as well as providing a food source for shorebird and ray species (Thrush et al. 1991, 2006; Sandwell et al. 2009; Jones et al. 2011a; Pratt et al. 2014). A. stutchburyi populations are declining in some areas, however, likely due to adverse effects associated with sedimentation, pollution and over-harvesting (De Luca-Abbott 2001; Norkko et al. 2006; Cummings et al. 2007).

In addition to chronic stressors adversely affecting bivalve populations, mass mortality events may have severe consequences on ecosystem services or biodiversity (e.g. Jonsson & Andé 1992; Horiguchi 1995). In New Zealand, a mass mortality event significantly impacted bivalve beds in Whangateau Harbour (North Island) in February 2009, an estuary with an abundant A. stutchburyi population. It was estimated that c. 60% of the A. stutchburyi population was affected by the die-off, amounting to a loss of c. 250 million clams in an estuary 7.5 km² in area (Ministry of Fisheries 2009). This event was attributed to a combination of infection by a coccidian parasite and a mycobacterium, and high temperatures coinciding with mid-afternoon low tides (Ministry of Fisheries 2009). Similar mass mortality events have also been observed elsewhere, for example, in Germany, where the bivalve was Cerastroderma edule (Thieltges 2006), but the effects of these events on the wider food web are largely unknown, despite the fact that bivalve populations may provide an important food source for higher trophic levels, such as shorebirds, rays and other fish (e.g. Thrush et al. 1991; Le Port 2003). For example, a reduction in both the abundance and condition of cockles (C. edule) in the Dutch Wadden Sea due to shellfish dredging has been linked to the loss of red knots (Calidris canutus islandica) from the area (van Gils et al. 2006).

Our aim in this study was to evaluate the impact of a mass mortality event on an estuarine ecosystem, focusing on its effects on: (i) the grazing pressure exerted by A. stutchburyi, which is potentially an important ecosystem service regulating pelagic phytoplankton biomass; (ii) the estuarine food web. In order to do so, the population filtration rate was calculated for before and after the mass mortality event and compared to other key ecosystem parameters (water residence time and phytoplankton production time) to assess the effect of the mortality event on filtration rates at the ecosystem scale. Moreover, a mass-balance trophic model (Ecopath) was configured to quantify functional attributes (e.g. biomass, production and consumption) of the entire system prior to the mass mortality event, and also to allow the establishment of a possible A. stutchburyi biomass range to compare with the biomass following the mass mortality event, and to determine the relative importance of phytoplankton, microphytobenthos and detritus in sustaining estuarine consumers. A mass-balance trophic model for a New Zealand estuary has not previously been published, making this a key case study on the evaluation of estuarine food webs.

Methods

Study area

Whangateau Harbour is located on the north-east coast of New Zealand (36°19′S, 174°46′E), 60 km north of the city of Auckland (Figure 1). It is representative of a common estuary type in New Zealand (Hume et al. 2007), and has been the focus of a number of studies, which provide a substantial database (reviewed by Kelly 2009; Townsend et al. 2010) suitable for parameterisation of a trophic model and for ecosystem scale calculations. Whangateau is a small barrier-enclosed lagoon (area 7.5 km², mean depth 1.6 m), which has extensive intertidal sandflats (85% of total estuary area), is well-flushed (the tidal prism is 81% of total estuary spring-tide volume) and receives little freshwater input (Hume et al. 2007). Tides are semi-diurnal, with a range of c. 2 m. Substrate transitions from sand in the main estuary body to mud in the sheltered parts of the arms, and the resulting complex and varied habitat supports diverse benthic macrofauna and fish communities (e.g. Francis et al. 2005; Alfaro 2006).

Figure 1. Whangateau Harbour, New Zealand. Land is shaded dark grey, intertidal flats are white and subtidal areas are light grey.

Austrovenus stutchburyi grazing pressure in Whangateau Harbour

Estuary-wide A. stutchburyi bed clearance rates were estimated from the density-dependent empirically derived equations given in Jones et al. (2011b) and surveyed densities before and after the mortality event (Walshe et al. 2006; Pawley & Ford 2007; Pawley 2011). These surveys focused on the high-density beds; for the rest of the intertidal area, A. stutchburyi density was assumed to be 50 ind. m⁻² (Townsend et al. 2010). The bed clearance rates measured by Jones et al. (2011b) used clams of a size (20–25 mm shell height) similar to those in Whangateau prior to the mass mortality event, and so the equations are likely to be appropriate for this period. Following the mass mortality event, however, mean clam size was reduced in high-density beds from 23.4 to 20.5 mm (Pawley 2011). To account for the reduction in size on bed clearance rates, we used a previously determined allometric relationship between biomass and shell length for A. stutchburyi (R² = 0.93, n = 140; C. Pilditch unpub. data), and the relationship between clearance rate and biomass for the European cockle (C. edule) derived by Smaal et al. (1997):where M = ash-free dry weight (g) and L = shell length (cm), andwhere CR = clearance rate (m³d⁻¹ind⁻¹).

To calculate A. stutchburyi population clearance rate, the bed clearance rates were scaled by the estuary intertidal area and converted to the volume of water cleared per day. We assumed that A. stutchburyi feed for 8 h d⁻¹ because tides are semi-diurnal and beds are located in the mid-tide zone. This feeding period is consistent with previous research that established that A. stutchburyi feed for 4 h per tide (Beentjes and Williams 1986).

We scaled the A. stutchburyi population clearance rates pre- and post-mortality event to the estuary volume to calculate bivalve clearance time, and compared these to the residence time and rate of phytoplankton turnover to quantify the effect of the mortality event on bivalve grazing pressure. Bivalve clearance time was calculated as follows:where BC_T = clearance time (d), E_V = estuary volume (m³) and Q_B = volume of water cleared by A. stutchburyi in 1 day (m³d⁻¹). Water residence time (R_T; d) was estimated from E_v and the volume of water exchanged between the estuary and the coast in 1 day (Q_p; m³d⁻¹). Estuary volume and the tidal prism were derived from NIWA’s Coastal Explorer (retrieved from http://wrenz.niwa.co.nz/webmodel/coastal on 15 August 2011). Phytoplankton production time (PP_t; d) was estimated from areal phytoplankton biomass (g C m⁻²) and production (g C m⁻²d⁻¹) values for other North Island New Zealand estuaries (Vant & Budd 1993; Gallegos & Vant 1996; Wilcock & Kemp 2003).

Trophic model of Whangateau Harbour

Model description

Food web models are used to explore the energy flow and interactions between multiple trophic levels (Gaedke 1995; Fulton et al. 2003; Christensen et al. 2008). Ecopath is a well-established mass-balance (i.e. steady state) food web model, which has previously been applied to a wide range of coastal systems (Christensen & Pauly 1992), including estuaries (Milessi et al. 2010; Byron et al. 2011a). The model estimates the biomass and food consumption of species (or groups of species) in an ecosystem using a parameterisation based on an assumption that a mass-balance is achieved over an arbitrary period, typically 1 year (Christensen et al. 2008).

Parameterisation is based on two governing equations; the first equation describes production for each group i (which may be a single species or a group of species that share a common functional role), whereby mortality of a prey group (a reduction in that state variable) equates to consumption for a predator group (increase in the state variable), so that:where P_i is production for a group i, Y_i is total fishery catch, B_i is biomass, M2_i is the total predation rate, E_i is net migration (emigration–immigration), BA_i is biomass accumulation, and EE_i is ‘ecotrophic efficiency’, which is the proportion of production that is used within the system. The term given by P_i (1 – EE_i) is also sometimes described as ‘other mortality’.

Ecopath consists of a series of linear equations (of the form of the equation above) for each group in the system, which are then solved for one of the following four parameters: biomass (B_i), production/biomass ratio (P_i/B_i), consumption/biomass ratio (C_i/B_i) or ecotrophic efficiency (EE_i). The other three parameters, along with catch, net migration, biomass accumulation, assimilation and diet composition, must be assigned for each group. After ensuring mass-balance between groups, the energy balance within a group is obtained using the second governing equation:where C_i is consumption, P_i is production, R_i is respiration and UA_i is unassimilated consumption. A balanced model is defined as one in which EE < 1 for all model groups. If EE ≥ 1 for one or more groups, the model is deemed unbalanced (because consumption of a group cannot be greater than its production), and the input parameters must be adjusted to obtain a balanced solution (Christensen and Walters 2004).

Model input, balancing and outputs

The food web of Whangateau Harbour was defined with 16 groups, where species that shared similar functional roles were combined together in a single functional group. A full list of species included in each functional group is given in Table 1 and the groups are listed in Table 2. All available information on biomass, P/B and C/B ratios and diet composition for each of the model groups was assembled from a significant body of literature (including c. 50 scientific articles, c. 15 reports, c. 5 university theses and authors’ unpublished data) and converted into appropriate units for the model, that is, g carbon m⁻² for biomass and g carbon m⁻²yr⁻¹ for production and consumption. The origin of the input data is given in Table 1 and described in detail in Jones (2011). Quantitative data on diet composition for each of the species/groups were scarce. Where only qualitative estimates were available, the proportion of prey in the diet of a group was assumed to be proportional to the fraction that its biomass contributed to the total biomass of all the prey items for the group (Pinkerton et al. 2008). Unassimilated consumption (i.e. the fraction of food that is egested and directed to detritus) was assumed to be 0.3 for herbivores and 0.2 for carnivores (Christensen et al. 2008). Recreational A. stutchburyi harvest was also included as a fishery term in the model based on data in Kearney (1999).

Table 1. Input data sources for groups in Whangateau Harbour trophic model.

Model group	Species	References
(1) Shorebirds	New Zealand dotterel (Charadrius obscurus), banded dotterel (Charadrius bicintus), pied oystercatcher (Haematopus ostralegus), variable oystercatcher (Haematopus unicolor), bar-tailed godwit (Limosa lapponica baueri), pied stilt (Himantopus himantopus), black-backed gull (Larus dominicanus), red-billed gull (Larus scopulinus), pied shag (Phalacrocorax varius varius)	Battley et al. (2007), Bradford-Grieve et al. (2003), Dowding and Chamberlin (1991), Heather and Robertson (2000), Jiang and Gibbs (2005), Kelly (2009), Lundquist and Pinkerton (2008), Medway (2000), Parrish (2000, 2001, 2002, 2006), Pinkerton et al. (2008) and Thrush et al. (1994)
(2) Rays	Short-tailed stingray (Dasyatis brevicaudata), eagleray (Myliobatis tenuicaudatus)	Chubb et al. (1981), Curtis and Shima (2005), de Paula et al. (1993), Francis (1988), Francis et al. (2005), Froese and Pauly (2005), Grace (1972), Haedrich and Merrett (1992), Jiang and Gibbs (2005), Kelly (2009), Le Port (2003), Lundquist and Pinkerton (2008), Morrison and Carbines (2006), Morrison et al. (2002), NIWA (2011), Palomares and Pauly (1998), Pinkerton et al. (2008), Russell (1983), Taylor and Willis (1998), Thrush et al. (1991, 1994), Usmar (2010), Webb (1973) and Wolff et al. (1998)
(3) Piscivorous fish	Kahawai (Arripis trutta), spotted stargazer (Genyagnus monopterygius), kingfish (Seriola lalandi), barracouta (Thyrsites atun), John dory (Zeus faber)
(4) Invertebrate-feeding fish	Rock cod (Acanthoclinus fuscus), yellow eyed mullet (Aldrichetta forsteri), red moki (Cheilodactylus spectabilis), red gurnard (Chelidonichthys kumu), hiwihiwi (Chironemus marmoratus), triplefin spp. (Forsterygion spp.), spotty (Notolabrus celidotus), banded wrasse (Notolabrus fucicola), snapper (Pagrus auratus), leatherjacket (Parika scaber), small bastard red cod (Pseudophycis breviuscula), yellow belly flounder (Rhombosolea leporine), sand flounder (Rhombosolea plebeian) clingfish (Trachelochismus melobesia), goatfish (Upeneichthys lineatus)
(5) Planktivorous fish	Koheru (Decapterus koheru), trevally (Pseudocaranx dentex), blue maomao (Scorpis violacea)
(6) Herbivorous fish	Parore (Girella tricuspidata), piper (Hyporhamphus ihi)
(7) Carnivorous invertebrates	Whelks (Cominella spp.), oyster borer (Lepsiella scobina)	Ansell (2001), Brey (2001), Cammen et al. (1984), César and Armendáriz (2007), Cummings et al. (1997), Grant and Hay (2003), Hartill et al. (2000), Hewitt (2008), Hicks (1985), H. Jones unpublished data, Kearney (1999) Kelly (2009), Lin et al. (2007), Lundquist and Pinkerton (2008), Ministry of Fisheries (2009), Needham et al. (2010), Manickchand-Heileman et al. (1998), Morrisey et al. (1999), Pawley and Ford (2007), Pawley (2011), Rybarczyk and Elkaim (2003), Needham et al. (2014), Stevens et al. (2006), Stewart and Creese (2004), Townsend et al. (2010), Thrush et al. (2003) and Walshe et al. (2006)
(8) Austrovenus	Little-neck clam/cockle (A. stutchburyi)
(9) Suspension feeding invertebrates	Pipi (Paphies australis)
(10) Deposit feeding invertebrates	Mud crab (Austrohelice crasse), wedge shell (Macomona liliana), nut shell (Nucula hartvigiana), polychaete worms, amphipods, benthic copepods, echinoderms
(11) Zooplankton	Various	Bradford-Grieve et al. (1993), Brey (2001), Froneman (2001), Fulton (1984), Giles (2002), Lara-Lopez and Neira (2008), Leandro et al. (2007) and Uye and Liang (1998)
(12) Macroalgae	Flapjack (Carpophyllum flexuosum), sea lettuce (Ulva sp.), Neptune’s necklace (Horomosira banksii), coralline algae (Corallina spp.)	Byron et al. (2011b), Pinkerton et al. (2008), Taylor et al. (2008), Tomczak et al. (2009) and Ortiz and Wolff (2002)
(13) Microphytobenthos	Various	Cahoon and Safi (2002), Gallegos and Vant (1996), Halliday and Cummings (2009), H. Jones (unpublished data), MacIntyre et al. (1996), Needham et al. (2014) and Vant and Budd (1993)
(14) Phytoplankton	Various	Scarsbrook (2008), Gallegos and Vant (1996), Vant and Budd (1993), Vant and Safi (1996) and Wilcock and Kemp (2003)
(15) Bacteria	Various	Lundquist and Pinkerton (2008)
(16) Detritus		Kelly (2009), H. Jones (unpublished data) and Rybarczyk and Elkaim (2003)

Table 2. List of abbreviations for model parameters, model groups and summary statistics.

Abbreviation	Denotes	Units
Model parameters
B	Biomass	g C m^−2
P	Production	g C m^−2yr^−1
C	Consumption	g C m^−2yr^−1
R	Respiration	g C m^−2yr^−1
A	Assimilation	g C m^−2yr^−1
EE	Ecotrophic efficiency	–
UA	Unassimilated consumption	–
Model groups
Bir	Birds	–
Ray	Rays	–
FPi	Piscivorous fish	–
FIn	Invertebrate-feeding fish	–
FPl	Planktivorous fish	–
FHe	Herbivorous fish	–
BCa	Carnivorous benthic invertebrates	–
Aus	A. stutchburyi	–
BSu	Benthic suspension feeding invertebrates	–
BDe	Benthic deposit feeding invertebrates	–
Zoo	Zooplankton	–
Mac	Macroalgae	–
MPB	Microphytobenthos	–
Phy	Phytoplankton	–
Bac	Bacteria	–
Det	Detritus	–
Aus harvest	Recreational harvest of A. stutchburyi	–
Summary statistics
TPP	Total primary production	g C m^−2yr^−1
TR	Total respiration	g C m^−2yr^−1
TB	Total biomass	g C m^−2
TST	Total system throughput	g C m^−2yr^−1
TL	Trophic level	–
MTI	Mixed trophic impact	–

Input data were examined to ensure values conformed to simple ecological and physiological rules. For example, biomass, and production and consumption ratios should decrease with increasing trophic level (Link 2010). The input parameters for biomass, P/B and C/B required only minor adjustments to the initial estimates for C/B for birds and some benthic invertebrate groups to balance the model. Model output was examined to ensure that values for respiration and assimilation were also within realistic ranges. Ratios of respiration to assimilation were <1, and decreased with trophic level. Ratios of respiration to biomass, which are an expression of the activity level of a group, were also within the ranges reported in the literature (Bradford-Grieve et al. 2003; Jiang & Gibbs 2005; Christensen et al. 2008).

Ecopath outputs fractional trophic levels (TLs) for each model group to recognise that predators can consume prey across a range of TLs. Producers and detritus have, by definition, a TL of 1 and TL for consumers is calculated as 1 + the weighted average of the prey TLs (Christensen et al. 2008). Summary statistics provide overall measures of total system throughput (i.e. the sum of all flows through each group), consumption, respiration, export and flows to detritus, as well as total and primary production, and gross efficiency of any fisheries (i.e. catch/net primary productivity). The mixed trophic impact (MTI) routine in Ecopath evaluates the impact that a very small increase in biomass of one group will have on other groups in the system.

Post-mortality event bivalve biomass and possible range in A. stutchburyi biomass

We assume that the balanced Ecopath model represents the estuary food web before the mass mortality event. A possible range for A. stutchburyi biomass was established using the model to place the reduction in A. stutchburyi following the mass mortality event (estimated to be 1.58 g C m⁻²; Jones 2011; Pawley 2011) in the context of a theoretical minimum biomass. To establish the possible biomass range, we adapted a process that has previously been used to calculate ecological carrying capacity for harvested shellfish in other systems (Jiang & Gibbs 2005; Byron et al. 2011a). In these instances, shellfish biomass is increased in consecutive model runs until the model is no longer balanced (i.e. Ecotrophic efficiencies [EE] for one or more groups is greater than 1), and the point just prior to this is defined as the ecological carrying capacity. We first decreased A. stutchburyi biomass in consecutive model runs until the model was no longer balanced, to quantify the minimum A. stutchburyi biomass the system might be able to support, and then increased A. stutchburyi biomass in consecutive model runs to quantify the ecological carrying capacity.

Results

Austrovenus stutchburyi grazing pressure in Whangateau Harbour

Bivalve bed clearance rate was calculated to range from 94 to 276 L h⁻¹m⁻² prior to the mass mortality event, but decreased to between 66 and 174 L h⁻¹m⁻² after the mass mortality event (Table 3). Before the mass mortality event the total volume of water cleared by A. stutchburyi per tidal cycle (i.e. the population clearance rate) was estimated to be between 565,100 and 1,653,300 m³, which is much less than the tidal prism (9,491,105 m³). After the mass mortality event, the population clearance rate was reduced slightly to between 529,600 and 1,391,500 m³. The water cleared by A. stutchburyi before the mass mortality event was equivalent to c. 9% of total harbour volume, and after the mass mortality event it was reduced slightly to c. 8% of total harbour volume.

Table 3. Whangateau Harbour A. stutchburyi population clearance rates calculated for three flow speeds (2, 5 and 15 cm s⁻¹). Mean A. stutchburyi density and size are as reported in Pawley (2011).

	Before mass mortality event			After mass mortality event
	High-density beds	Low-density areas	Total harbour	High-density beds^a	Low-density areas	Total harbour
Density (ind m^−2)	495	50	NA	372	50	NA
Size (mm)	23.4	23.4	NA	20.5	23.4	NA
Bed clearance rate (L h^−1m^−2)^b
2 cm s^−1	94.4	36.9	131.2	65.5	36.9	102.3
5 cm s^−1	48.1	23.1	71.2	34.3	23.1	57.4
15 cm^−1	276	49.5	325.5	173.8	49.5	223.2
Population clearance rate per tidal cycle (m^3)^c
2 cm s^−1	242,000	704,800	946,800	168,000	704,800	872,800
5 cm s^−1	123,400	441,700	565,100	87,900	441,700	529,600
15 cm^−1	707,700	945,600	1,653,300	445,900	945,600	1,391,500
Average	357,700	697,400	1,055,100	233,900	697,400	931,300

^aBCR reduced by 22% for high-density beds after mass mortality event to account for reduction in mean size.

^bBed clearance rate (BCR) at each flow speed calculated from equations in Jones et al. (2011b).

^cPopulation clearance rate calculated by scaling BCR by area (641,500 m² for high-density beds and 4,780,000 m² for low-density areas) and assuming A. stutchburyi feed for 4 h per tidal cycle.

The bivalve clearance time (i.e. the estuary volume divided by the volume of water cleared by A. stutchburyi per day) was 5.5 and 6 days, before and after the mass mortality event, respectively (Table 4). This is much greater than both the water residence time for Whangateau Harbour (0.61 days) and typical phytoplankton production time (1.4 days; Vant & Budd 1993; Gallegos & Vant 1996; Wilcock & Kemp 2003).

Table 4. Ecosystem scale parameters for Whangateau Harbour.

Physical characteristics and turnover rates	Value
Estuary volume at spring tide^a (m^3)	11,663,589
Spring tidal prism^a (m^3)	9,491,105
A. stutchburyi population clearance rate per tidal cycle (m^3)
Before mass mortality event	1,055,000
After mass mortality event	961,200
Bivalve clearance time (days)
Before mass mortality event	5.54
After mass mortality event	6.07
Water residence time (days)	0.61
Phytoplankton production time (days)	1.40

^aSource: NIWA Coastal Explorer.

Trophic modelling

Pre-mortality event food web

The trophic model of Whangateau Harbour is shown in Figure 2, and input/output parameters and summary statistics are described in Tables 5 and 6, respectively. The benthic groups equated to 87% of total living biomass (A. stutchburyi accounted for 30%), and of the primary producers, microphytobenthos were by far the largest group, making up 90% of total primary producer biomass. Consumption by the benthic groups (92 g C m⁻²yr⁻¹) far outweighed that of the pelagic groups (4 g C m⁻²yr⁻¹) with A. stutchburyi accounting for 43% of total consumption and exceeded only by that of deposit feeding benthic invertebrates (50% of total).

Figure 2. Trophic model for Whangateau Harbour. Circle size is proportional to biomass. See Table 1 for group name abbreviations.

Table 5. (a) Balanced model input and output parameters, and (b) diet matrix, where numbers represent percentage prey (rows) taken by predators (columns).

	(a) Input and output parameters										(b) Diet Matrix
Group	TL	B	P/B	C/B	EE	P/C	UA	R/A	P/R	R/B	Bir	Ray	FPi	FIn	FPl	FHe	BCa	Aus	BSu	BDe	Zoo	Bac
Bir	3.234	0.005	0.20	90	0.000	0.002	0.2	0.997	0.003	71.8
Ray	3.133	0.128	0.11	3.90	0.000	0.028	0.2	0.965	0.037	3.01
FPi	3.930	0.017	0.27	3.50	0.000	0.077	0.2	0.904	0.107	2.53
FIn	3.211	0.132	0.56	6.67	0.659	0.084	0.2	0.895	0.117	4.78	1.9		68.7
FPl	3.180	0.010	0.40	6.99	0.774	0.057	0.2	0.928	0.077	5.19			5.2
FHe	2.058	0.033	0.34	13.48	0.854	0.025	0.3	0.964	0.037	9.10			16.1
BCa	3.195	0.110	1.60	9	0.134	0.178	0.2	0.778	0.286	5.60	1.7	0.5		1.6
Aus	2.100	3.380	2.18	12	0.214	0.182	0.3	0.740	0.350	6.22	49.5	85.0		49.0			50.1
BSu	2.100	0.170	2.18	12	0.136	0.182	0.3	0.740	0.350	6.22		0.7		2.5			2.5
BDe	2.300	3.200	3.00	15	0.121	0.200	0.3	0.714	0.400	7.50	46.9	13.8	8.0	46.8	20.0		47.4
Zoo	2.150	0.022	20	80	0.279	0.250	0.3	0.643	0.556	36.0			2.0	0.1	80.0	5.0		0.1	0.1
Mac	1.000	0.235	18		0.100											95.0
MPB	1.000	3.000	53.3		0.279													80.0	80.0	20.0	50.0
Phy	1.000	0.110	260		0.161													10.0	10.0		20.0
Bac	2.000	1.000	40	125	0.472	0.320	0.3	0.543	0.842	47.5								9.9	9.9	30.0	15.0
Det	1.000	151			0.609															50.0	15.0	100

Note: Figures in bold indicate those that are calculated by the model (i.e. model output).

TL, trophic level; B, biomass (g C m⁻²); P, production (g C m⁻²yr⁻¹); C, consumption (g C m⁻² yr⁻¹); EE, ecotrophic efficiency; UA , unassimilated consumption; R, respiration (g C m⁻² yr⁻¹); A = assimilation (g C m⁻²yr⁻¹). For group name abbreviations see Table 2. For derivation of input data see Table 1.

Table 6. Summary statistics for balanced models representing: original model for Whangateau Harbour (before the A. stutchburyi mass mortality event); the estuary at minimum A. stutchburyi biomass; and the estuary at ecological carrying capacity (i.e. maximum A. stutchburyi biomass).

Parameter	Original model	Min. A. stutchburyi biomass	Ecological carrying capacity
Sum of all consumption (g  C  m^−2)	221	189	363
Sum of all exports (g  C  m^−2)	96.0	112	22.4
Sum of all respiratory flows (g C  m^−2)	96.7	80.2	170
Sum of all flows into detritus (g  C  m^−2)	245	262	172
TST (g  C  m^−2)	659	643	727
Sum of all production (g  C  m^−2)	251	245	277
Mean trophic level of the catch	2.10	2.10	2.10
GE (catch/net primary production)	0.00011	0.00011	0.00051
TPP (g C  m^−2)	193	193	193
TPP/TR	1.99	2.40	1.13
Net system production (g C  m^−2)	96.0	112	22.4
TPP/TB (yr^−1)	16.68	21.65	8.24
TB/TST	0.018	0.014	0.032
TB (excluding detritus) (g  C  m^−2)	11.55	8.90	23.4
Total catches (g  C  m^−2)	0.022	0.022	0.099

TST, total system throughput; GE, gross efficiency; TPP, total primary production; TR, total respiration; TB, total biomass.

EE were low for some groups, such as the primary producers and benthic invertebrates (EE < 0.4), but were higher (>0.6) for other groups, such as invertebrate-feeding fish. Of the total flow into detritus, c. 60% was consumed by, and c. 40% exported from, the system (Table 5). Of the total system throughput (659 g C m⁻²yr⁻¹), 34% was consumed, 15% respired, 37% flowed to detritus and 15% was exported from the system (Table 6). Fractional TLs ranged from 1 for primary producers to 3.9 for piscivorous fish (Table 5). The greatest proportion (58%) of living biomass was in herbivores, whilst primary producers accounted for 29% of biomass, and only 3% of living biomass was at a trophic level of 3 or above.

The MTI analysis indicated that, of the top predators, birds and rays had little impact on any other group; piscivorous fish, however, had a large negative impact on their prey, and a positive impact on carnivorous benthic invertebrates, zooplankton and macroalgae (Figure 3). Negative impacts on prey species were obvious, for example, on zooplankton by planktivorous fish, and on phytoplankton and microphytobenthos by A. stutchburyi. Microphytobenthos had a positive impact on most groups, especially A. stutchburyi and other benthic suspension feeders. A. stutchburyi had relatively large impacts (either negative or positive) on most other groups, and the harvest of A. stutchburyi (i.e. recreational fishery) had little impact on A. stutchburyi or any other groups.

Figure 3. MTI analysis of Whangateau Harbour model. Direct and indirect impacts of a small increase in the biomass of one group (vertical axis) on all other groups (horizontal axis) are shown. Positive impacts are shown in black and negative in white. Impacts are relative and circle size is proportional to the size of the impact. See Table 2 for group name abbreviations.

Austrovenus stutchburyi biomass range and post-mortality event biomass

A. stutchburyi biomass was decreased in small increments from the original model value of 3.38  to 0.72 g C m⁻², at which point the model became unbalanced (i.e. EE > 1) (Table 7). A. stutchburyi biomass following the 2009 mass mortality event (1.58 g C m⁻²) was nearly three times greater than the minimum possible biomass as determined in this study (i.e. 0.73 g C m⁻²).

Table 7. Effect of changes in A. stutchburyi biomass on Whangateau Harbour trophic model under two scenarios: (1) decreasing biomass whilst maintaining a constant harvest at 0.022 g C m⁻²yr⁻¹, (2) increasing biomass (and harvest proportional to biomass) to calculate ecological carrying capacity.

Multiplier	Biomass (g C m^−2)	Harvest (g C m^−2yr^−1)	Mass-balance changes
Original model representing estuary before mass mortality event
1	3.38	0.022
Scenario 1: biomass decreased and harvest maintained at 0.022 g C m^− ^2yr^− ^1
0.8	2.70	0.022	Balanced
0.6	2.03	0.022	Balanced
0.47	1.58*	0.022	Balanced
0.25	0.845	0.022	Balanced
0.22	0.73	0.022	Balanced
0.21	0.72	0.022	Unbalanced: A. stutchburyi EE = 1.005
Scenario 2: biomass and harvest increased to calculate ecological carrying capacity
2.0	6.76	0.044	Balanced
4.0	13.52	0.088	Balanced
4.5	15.21	0.099	Balanced
5.0	16.90	0.110	Unbalanced: Microphytobenthos EE = 1.09

Note: Minimum and maximum A. stutchburyi biomass is highlighted in bold and A. stutchburyi biomass after 2009 mass mortality event is indicated by *.

Increasing A. stutchburyi biomass in small increments suggests that the biomass could increase by 4.5 times without exceeding the ecological carrying capacity of Whangateau Harbour (Table 7). When biomass was increased further than this the grazing pressure on microphytobenthos was too great (i.e. the model became unbalanced).

Total system throughput (i.e. the sum of all flows through each group) decreased slightly (by c. 2%) when A. stutchburyi biomass was reduced to the minimum possible, but increased by 10% when A. stutchburyi biomass was increased to ecological carrying capacity (Table 6). The sum of all exports increased by c. 15% at the minimum A. stutchburyi biomass, but decreased by c. 75% at ecological carrying capacity.

Discussion

Bivalve grazing pressure pre- and post-mortality event

Our calculations of A. stutchburyi population clearance rate in Whangateau Harbour reveal that the volume of water cleared by the bivalve population in each tidal cycle is much less than the tidal prism, and the ecosystem scale parameters reveal that the bivalve clearance time (c. 6 days) is much greater than both the water residence and phytoplankton production times (0.61 and 1.4 days, respectively). The impact of the 2009 mass mortality event on the grazing pressure exerted by the A. stutchburyi population was therefore probably less severe than the reduction in mean size and density of the population would suggest. Jones et al. (2011b) found that A. stutchburyi bed clearance rates do not scale proportionally with A. stutchburyi density, as individual clearance rates are suppressed at high bed densities. Thus the reduced clam density following the mass mortality event is partially offset by higher individual clearance rates. Regardless, the population clearance rate calculations for both before and after the mass mortality event suggest that A. stutchburyi do not exert a significant top-down control on phytoplankton in the harbour. In fact, the trophic model indicates that pelagic phytoplankton play a limited role in the estuary food web.

Pre-mortality event food web

The trophic model shows a system dominated by the benthic groups, particularly A. stutchburyi, deposit feeding benthic invertebrates and microphytobenthos. The latter group was the most important primary producer in terms of biomass and productivity, and A. stutchburyi and deposit feeding benthic invertebrates were responsible for the vast majority of total consumption. Furthermore, relatively large impacts on numerous other groups for both microphytobenthos and A. stutchburyi confirm that these groups are important components of the system, particularly in sustaining higher TLs.

In estuarine ecosystems where pelagic phytoplankton populations are low (as in this case) and there is limited top-down control, then bottom-up control on phytoplankton productivity may be more important (Prins et al. 1998). It is unlikely that light is limiting pelagic phytoplankton growth in such a shallow estuary with a highly abundant and productive microphytobenthos community (Cahoon and Safi 2002), but water column nutrient concentrations in estuaries similar to Whangateau tend to be low ( and concentrations are typically <0.01 mg/L; e.g. Tay et al. 2012). This highlights the functionally important link between benthic macrofauna, such as A. stutchburyi, and microphytobenthos. Microphytobenthos have been recognised as important primary producers in shallow estuarine systems worldwide (see reviews by MacIntyre et al. 1996; Miller et al. 1996), and in many of these systems there are abundant natural or cultivated suspension feeding bivalve populations. Experimental studies have shown that resuspended microphytobenthos may represent an important food source for bivalves in many estuarine systems, for example in Marennes-Oleron Bay in France (Kang et al. 1999), in Halifax, Nova Scotia (Grant et al. 1990) and in Kwangyang Bay in Korea (Kang et al. 2006). The bivalves, in turn, can enhance microphytobenthos productivity by increasing nutrient availability through excretion and bioturbation of nutrient-rich porewaters (Sandwell et al. 2009; Jones et al. 2011a). Our model results support these studies by highlighting the important role of microphytobenthos in our study system. Removal of microphytobenthos from the model revealed that phytoplankton and detrital production alone were not enough to sustain the consumer groups (i.e. a balanced solution was not obtained with benthic invertebrates and zooplankton feeding only on phytoplankton and detritus). It would seem that in this estuary (and likely in other similar systems) suspension feeding bivalves are not a ‘natural eutrophication control’ (sensu Officer et al. 1982); rather they appear to perform the opposite function by stimulating benthic primary production.

The trophic model also shows a system characterised by high flows to detritus for some groups and relatively high detritus export (c. 40%) from the system. Some other coastal lagoon systems show similar characteristics, for example, high microphytobenthos productivity and significant export of organic material to the open ocean (Haese & Pronk 2011), as well as high flows to detritus and dominance of benthic groups (Milessi et al. 2010). The ‘detritus export’ in our model was 96 g C m⁻²yr⁻¹, which includes detritus that is buried in the sediment and that which is physically exported out of system. Sedimentation rates in Whangateau Harbour are likely to be relatively low across most of the estuary due to the limited development of the catchment, and the well-flushed nature of the estuary (Kelly 2009). However, mangrove expansion indicates that net sedimentation is occurring in the more sheltered reaches of the estuary, and suggests that detritus burial is likely in these areas at least (Townsend et al. 2010).

Whangateau Harbour is ebb-dominated (Grace 1972) and has a large tidal prism (80% of estuary volume), and thus the system likely conforms to the ‘outwelling hypothesis’, where a substantial amount of detritus may be advected to the coastal environment (Odum et al. 1979; Childers et al. 2000). The flux of organic matter from Whangateau has not been measured directly, but comparison with studies from other estuaries reveals that an export of c. 100 g C m⁻²yr⁻¹ is not unreasonable, for example, 50–80 g C m⁻²yr⁻¹ for Manukau Harbour, New Zealand (Vant et al. 1998), 453 g C m⁻²yr⁻¹ for North Inlet, South Carolina (Dame et al. 1986), and c. 320 g C m⁻²yr⁻¹ for Swartkops Estuary, South Africa (Baird et al. 1987). This export of organic material to the coastal ocean is likely to provide an energy source for plankton and benthic consumers, which may significantly influence coastal productivity (Dame and Allen 1996; Schlacher et al. 2009; Savage et al. 2012).

Post-mortality event bivalve biomass in the context of possible A. stutchburyi biomass range

In this study, a possible range for A. stutchburyi biomass was established using the trophic model to place the reduction in A. stutchburyi following the mass mortality event in the context of a theoretical minimum biomass. We found that A. stutchburyi biomass was able to be substantially decreased, to 22% of its original value, without unbalancing the model. The minimum biomass (0.73 g C m⁻²) was less than the A. stutchburyi biomass estimated for the estuary after the mass mortality event (1.58 g C m⁻²), suggesting that the system is able to assimilate such a reduction in clam biomass. In other ecosystems (particularly those with commercially fished shellfish populations), low bivalve stocks are believed to have contributed to mass mortalities of shorebirds, such as oystercatchers (e.g. Atkinson et al. 2003). However, to the best of our knowledge, there were no adverse effects reported on any other components of the ecosystem following the mass mortality event in 2009. Furthermore, recent surveys (in 2013 and 2014) suggest that the A. stutchburyi population in Whangateau Harbour has recovered (Berkenbusch et al. 2015). We recognise that incorporating spatial and temporal dynamics into the modelling process might reveal that the system is more sensitive than the current results suggest, and this study does not address other potential effects of the mass mortality event. For example, there may have been spatially and temporally localised effects associated with organic matter decomposition, and changes in sediment composition associated with an increase in shell hash, that cannot be resolved by the ecosystem-level calculations and mass-balance modelling carried out in this study.

We estimated the maximum A. stutchburyi biomass (i.e. ecological carrying capacity) to be 4.5 times greater than the pre-mass mortality event level (fuelled by the high microphytobenthos biomass in the ecosystem), perhaps supporting anecdotal evidence of historically larger individual bivalves and bivalve populations (Marsden & Adkins 2010). A. stutchburyi populations in estuaries around New Zealand in recent decades may have been lower they were historically, due to habitat degradation, sedimentation and over-fishing (e.g. Cummings et al. 2007; Marsden & Adkins 2010). In other systems, detritus has been shown to make up a significant proportion of suspension feeding bivalves’ diet (Kang et al. 1999; Decottignies et al. 2007). If detritus were to be a more important component in A. stutchburyi diet (than is assumed in this model), then there is the potential for even greater ecological carrying capacity; in the current model, the detritus is not a large food source. However, as a substantial amount of detritus, phytoplankton and resuspended microphytobenthos is likely to be removed from the estuary by tidal flushing, ecological carrying-capacity calculations should be viewed with caution.

Conclusions

This study highlights the important role of microphytobenthos in sustaining the food web in Whangateau Harbour, and suggests that phytoplankton production, and the interactions between suspension feeding bivalves and phytoplankton are of reduced importance. Instead, the A. stutchburyi population plays a key role in facilitating the transfer of energy from microphytobenthos to higher TLs. Our results are likely to be applicable to other shallow, well-flushed estuaries that are characterised by extensive intertidal flats, and abundant microphytobenthos and suspension feeding bivalve populations.

The trophic model indicates that A. stutchburyi biomass may be reduced or increased considerably without exceeding the ecologically carrying capacity or resulting in food limitation for higher TLs. The reduction in A. stutchburyi biomass in Whangateau Harbour following a mass mortality event in 2009 was within the range identified in the model, suggesting it may have had a limited effect on trophic dynamics. However, this study does not address other potential effects (e.g. from organic matter decomposition or changes in sediment composition) of the mass mortality event.

The trophic model provided a framework to synthesise available knowledge from a wide range of sources, quantify the relative importance of ecosystem components, which is difficult to accomplish any other way when components are coupled in a complex and nonlinear way, and finally, to identify knowledge gaps. Although this study incorporated a comprehensive overview of available data, there were still limited data available for quantifying fish biomass in Whangateau Harbour (and in New Zealand estuaries in general). Further complicating parameter estimation, a number of fish species migrate to and from estuaries, or use estuarine habitats as nursery grounds before dispersing to open coast environments (Morrison et al. 2002; Francis et al. 2005). Finally, the model did not include temporal dynamics due to lack of time series data for validation. It is possible that dynamic features of the die-off event such as anoxia, sediment composition changes and rapid increase in detrital material might mean that the system is somewhat more sensitive to biomass changes than the current mass-balance model suggests.

This study highlights the importance of ensuring microphytobenthos productivity in shallow estuarine ecosystems is not compromised by catchment or estuarine management practices that increase turbidity and/or sedimentation (e.g. Pratt et al. 2013). Furthermore, the application of an ecosystem-level model that can identify key components of a system is a useful tool that can inform and prioritise monitoring programmes. National and international reporting programmes are often focused on environmental performance indicators, and results from ecosystem-level models may aid in the selection of indicator species to ensure that programmes include components of the ecosystem which have a significant impact on trophic dynamics.

Acknowledgements

We thank the following people for responding to requests for data and/or advice: Megan Carbines (Auckland Council), Richard Ford and Alan Moore (Ministry for Primary Industries), Mark Morrison, Malcolm Francis and Carolyn Lundquist (all National Institute of Water and Atmospheric Research). Matt Pinkerton (NIWA) provided helpful comments on an earlier draft of this manuscript. Associate editor: Associate Professor Andrea C Alfaro. HFEJ and CAP designed the study, HFEJ undertook the modelling, and wrote the manuscript. CAP, DPH, and KRB assisted in interpretation of results and editing the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Hannah F. E. Jones http://orcid.org/0000-0002-3231-4376

Additional information

Funding

This research was supported by a University of Waikato Doctoral Scholarship to HFEJ and CO1×1515 4.2.1 to CAP and KRB.