Variability in annual recruitment success as a determinant of long-term and large-scale variation in annual production of intertidal Wadden Sea mussels (Mytilus edulis)

Abstract

To understand the background of the strong variation and recent decline of stocks and production of mussels (Mytilus edulis) on tidal flats of the Wadden Sea, we analysed long-term (twice-annual for 26 years) and multi-station (15 sites) estimates of numbers, mean individual weights, biomass, and annual production on Balgzand, a 50-km² tidal-flat area in the westernmost part of the Wadden Sea (The Netherlands). Somatic production was estimated from summed growth increments of soft tissues per half-year period and expressed in ash-free dry mass (AFDM). In adults, positive values in spring/summer regularly alternated with negative values in autumn/winter, when up to ∼25% (mean: 14%) of individual weight gains in the preceding season were lost. No weight losses were observed during the first winter of the life of mussels. The 26-year mean of net somatic tissue production P amounted to 5.5 g AFDM m⁻² a⁻¹ at a mean biomass B of 3.2 g AFDM m⁻²; the ratio P/B varied strongly with age composition of the mussel population and ranged between 0.5 and 3.0 a⁻¹ (mean: 1.7). Within the restricted areas of mussel beds, mean biomass and annual production values were two orders of magnitude higher. In the Wadden Sea, mussel beds cover a typical 1% of extensive tidal flat areas. Numerical densities of recruits showed straight-line relationships with subsequent life-time year-class production. Once recruits had reached an age of ∼10 months, their numbers predicted subsequent production within narrow limits. Production per recruit averaged 0.21 g AFDM for 10-mo recruits and was not related to recruit density. Local variation in annual production varied strongly, with maximal values between mid-tide and low-tide level, where recruitment was also maximal. Production per recruit was higher at low than at high intertidal levels. Frequently failing recruitment is indicated as the main cause of declining mussel stocks in the Wadden Sea. As in other bivalve species, a declining frequency of the occurrence of cold winters appears to govern declining recruitment success and consequently declining production.

Introduction

Mussels (Mytilus edulis L.) play an important role in the ecosystem of coastal waters and estuaries such as the Wadden Sea (Dankers and Zuidema 1995). Though the abundance of mussels in the Dutch Wadden Sea varies strongly from year to year (Dankers and Koelemaij 1989), their biomass can account for no less than about a quarter of the total benthic biomass of the tidal flats (Beukema 1976) and even more in the subtidal zone (Dekker 1989).

Its populations play an important role in coastal ecosystems by promotion of sedimentation (Verwey 1952), and by effects on exchange processes at the sediment-water interface (Van Duren et al. 2006), on concentrations of nutrients and phytoplankton (Dame and Dankers 1988; Asmus and Asmus 1991; Dankers 1993), and on other intertidal animals by creation of a specific mussel-bed faunal assemblage (Dittmann 1990; Saier 2002). They are an important food source for such specialised bird species as oystercatchers (Hulscher 1996) and eider ducks (Swennen 1976). In some parts of the Wadden Sea, birds can consume substantial proportions of annual mussel production (Faldborg et al. 1994; Nehls et al. 1997). Mussel fishery dredges for seed and half-grown mussels to stock culture plots. In some years, this fishery severely affected natural mussel beds (Smit et al. 1998; Dankers et al. 2004). Because re-appearance of new intertidal beds proved to be a slow process (Dankers et al. 2001, 2004), fishery for shellfish is by now restricted to subtidal areas. Particularly in the western half of the Dutch Wadden Sea, mussel stocks on tidal flats remained low for the last 15 years (Steenbergen et al. 2005) and nowadays mussel beds are few there (see map 1 in Goudswaard et al. (2006)).

The outstanding role of mussels in the Wadden Sea was convincingly demonstrated by an unintended large-scale “natural experiment” (Beukema et al. 1998) in 1990, when intensive mussel fishery removed nearly all intertidal mussel beds from the Dutch Wadden Sea at a time when subtidal mussels were also scarce (Dankers 1993; Smit et al. 1998). The consequences were far-reaching and included such diverse effects in 1990/1991 as changing sediment character over vast areas from silty to sandy (Zwarts 2004), enhancing phytoplankton concentrations by lower grazing pressure, enhancing growth rates of other bivalves by elevated food levels, and increasing mortality rates in a number of other bivalve species serving as alternative food for bivalve consumers (Beukema 1993; Beukema and Cadée 1996). Some consequences were even disastrous, as food shortage caused high mortality in some specialised bird species (Camphuysen et al. 1996, 2002).

In the present paper, we emphasise one aspect of the role of mussels on the tidal-flat ecosystem of the Wadden Sea: their contribution to the benthic secondary production. Because abundance of mussels (like many other species) is highly variable on tidal flats both from place to place and from year to year (Beukema et al. 1983, 1993; Brinkman et al. 2002; Dare et al. 2004), we base our production estimates on a long-term (∼30 years) data series of numbers and weights of macrofaunal benthos at 15 sites on Balgzand, a 50-km² tidal-flat area in the westernmost part of the Wadden Sea. So far, published studies on mussel production were only short-term (1 or at most 2 years) and mostly referred to a small number of sampling sites.

The primary aim of the present study is to study the year-to-year variability in annual mussel production and to understand its underlying causes, such as recruit numbers and life-time production per individual recruit. Ultimately, we intend to study the relationship between total primary and secondary production in the area. Long-term data on phytoplankton concentration and primary production are already available for the westernmost inlet of the Wadden Sea (Cadée and Hegeman 2002). There is evidence that phytoplankton abundance is the prime-limiting factor for total zoobenthic biomass in estuarine areas (Herman et al. 1999), including Balgzand (Beukema et al. 2002). In some major species, however, recruitment success rather than food supply appears to be decisive for production and resulting biomass (cockle Cerastoderma edule: Beukema and Dekker 2006; Macoma balthica: Van der Meer et al. 2001; Dekker and Beukema 2007). Unfortunately, sufficiently extensive data sets to evaluate the relationships between primary and secondary production (and biomass) in marine benthos appear to be scarce. The only way to obtain relevant data on total zoobenthic production is by estimating secondary production for a number of years and at several sites, separately in all important individual species and put together these species data. We started this effort with studies of secondary production in cockles C. edule (Beukema and Dekker 2006) and M. balthica (Dekker and Beukema 2007). Papers to come will deal with the bivalves Mya arenaria and Ensis americanus.

Methods

Since the 1970s, the macrozoobenthos community has been monitored at 15 fixed sampling sites (12 transects of 1 km each and 3 squares of 30 × 30 m, for locations see Fig. 1) at a 50-km² tidal-flat area called “Balgzand”, in the westernmost part of the Dutch Wadden Sea (at about 53°N and 5°E). The monitoring program comprises twice-annual (late winter: mostly February/March, and late summer: mostly August/September) estimates of numbers and biomass of all macrobenthic (those retained on 1-mm sieves) animal species. Details on sampling sites and methods can be found in earlier publications (e.g., Beukema and Cadée 1997). The sampling sites were scattered (almost) randomly over the tidal flats of the Balgzand area and covered the entire intertidal depth range (compare Fig. 1b of Beukema and Cadée 1997), which runs on average from about +6 to −8 dm, compared with mean-tide level (MTL). The level of near-shore sites is generally above MTL and these sites are immersed for less than half of the time. Levels of more offshore sites are increasingly lower and the two farthest sites even remain permanently immersed for more than half of the low tides. The sampling sites cover a wide range of sediment composition. Mud (material < 60 μm) content of the top 5 cm of the sediment ranges from >10% near the SW coast to <1% in the far offshore parts of Balgzand.

Fig. 1

Map of Balgzand, showing location of the 15 sampling sites: 12 transects of 1 km each (marked 1–12) and 3 squares of 900 m² each (A–C). Thin lines are the borders of the intertidal area. Mussel beds were present along parts of transect-4 (in all years 1970 till summer 1990), transect-5 (in all years 1970–2005, except 1990/1991), transect-8 (infrequent and short periods), and transects-9 and -10 (frequent short periods)

Mussels were particularly numerous in the central part of Balgzand, where transects 4–11 (Fig. 1) in several years cut across extensive mussel beds. A map of Balgzand with the exact location of the mussel beds in early 2006 can be found in Goudswaard et al. (2006). The other half of the sampling sites were located in areas where mussels were encountered only occasionally, in the form of scattered individuals or small clumps. Most of these mussel-poor sites were located either at high levels (above or around MTL: transects 1 + 2 and A + B, respectively) or close to the mean low-water level (12 + C).

Along each of the 12 transects, 50 cores were taken at equal intervals (20 m), covering a total of 0.45 m² (summer) or 0.9 m² (winter) per transect. The 18–32 samples taken per square covered 1–1.8 m². In winter, larger samples were taken than in summer, because numerical densities were generally lower in winter than in summer and we needed a certain minimal number of individuals to estimate their mean weight with sufficient precision. In the present paper, we use data on M. edulis of the homogeneous 28-years series 1978–2005 (before 1978, regular annual sampling was largely restricted to late winter). All data on abundance are expressed in n m⁻², on biomass in g AFDM m⁻² , and on production in g AFDM m⁻² a⁻¹. All M. edulis were sorted from the sieved samples, measured and opened by short immersion in boiling water to remove the soft parts. These were dried (two or three days at 60°C in a ventilated stove), weighed (individually or per 1-mm shell-length class), and incinerated. Ash-free dry mass (AFDM) is dry weight minus ash weight.

Values for annual mean biomass B were calculated as means of the estimates of the two annual samplings (in March and August) in the year. However, for the 0-group (the spat of the year, which were not yet present at the time of the March sampling) the annual mean was calculated as the average of the biomass at the August and the subsequent March sampling. These estimates of B were used to calculate annual P/B ratio values (in a⁻¹) for separate age classes and for the total population.

Values of (somatic) production P were estimated following the weight-increment summation method as described by Van der Meer et al. (2005). We made separate calculations for half-year periods (from time t to t + 1) and for individual sampling sites. As annual estimates we used the sum of two successive seasonal estimates, the one for the March-August season plus the one for the subsequent August–March season. For each of the 15 sites and 52 half-year periods between March 1979 and March 2005, we defined production P as the product of mean numerical density (0.5 N_t + 0.5 N _t+1) and weight change (w _t+1 − w _t ), where N _t and N _t+1 are the numerical densities (n m⁻²) at the start and end of the production period, respectively, and w _t+1 and w _t are the mean soft-parts AFDM per individual at these times. Estimates of N and w were obtained separately per year class or (in old animals) groups of cohorts (cohort or year class: all individuals born in the same year).

In mussels older than ∼2 years, it proved to be difficult to separate all individuals into groups of the same age because their shells usually lacked clear annual rings beyond the first one. Like other investigators, we overcame this problem in assigning older individuals to specific cohorts by means of length-frequency curves. As an example, Fig. 4 shows that strong cohorts can easily be traced in the total population for about 5 successive years. Fortunately, during most of the period of our study, new strong year classes arose with long intervals of 4 years, compare Fig. 1c of Beukema et al. (2001): in 1975, 1979, 1983, 1987, and 1991. As a consequence, the length–frequency distributions of the infrequent strong cohorts are bordered by broad length ranges occupied by only few individuals, i.e. inevitable errors in the allotment of individual mussels to specific cohorts result in only minor errors in the mean values characterising strong cohorts.

We did not succeed in estimating production during the first growing season. The first estimates on N and w available were of spat-sized animals in August, a few months after their settlement. At that time a substantial part of the growth in the first growing season was already completed. A minimal estimate of spat production before August would equal the spat biomass found in August, but the true production (including their biomass at settlement and any production by spat that had died between settlement in spring and sampling in August) must have been higher, maybe twice the biomass observed in August. We arbitrarily adopted a value of 2 times spat biomass in August (i.e., the product of their numerical density times their mean individual weight in August) as a rough estimate of production of spat-sized animals in their first few months of life.

During the 1969–1974 period, a preliminary study was executed at a single mussel bed near transects-9 and -10 (Fig. 1). Annually in September or October, the size (in m²) of the major mussel patches within a 1,100-m² area (i.e., the maximal size of the bed, as observed in summer 1970) was estimated and six samples of 0.1 m² each were taken within this area. By allotting each of these samples to the patches with a density looking similar to the sampled spot (including all patches in this procedure), the total number of mussels within the bed area was roughly estimated with minimal effort. All sampled mussels were measured and their AFDM was assessed per mm-group. All mussels were allotted to year classes as described above.

Estimation of total-Balgzand mussel production by this method of sampling mussel beds one-by-one proved to be too time consuming for a long-term routine program. For long-term estimates of mussel production, we had to choose for using the data obtained by the twice-annual long-transect sampling program, though this method was not tailored to mussel sampling.

Results

Seasonal patterns in growth and production

Like other bivalves in the Wadden Sea, mussels are growing faster in the spring/summer season than in the remainder of the year. Nevertheless, we usually observed significant length increases between late-summer and late-winter, too (Fig. 2a: t tests: length increases amounted to several times their standard errors in all half-year periods). The detailed observations in the intensively-studied musselbed area around transects-9 and -10 pointed to continued fast growth in early autumn (between August and October) and reduced growth between October and March. The continuing length growth between August and March, however, did not generally result in positive contributions to production in the autumn/winter season. At the same length, mussels lost substantial parts of the weight of their soft parts between late-summer/early autumn and late-winter, body-mass indices declining usually from ∼6 to ∼4 mg cm⁻³. As a consequence, mean AFDM per individual only increased in autumn/winter if growth was sufficiently fast (i.e., the length ratio March/August should be >∼1.15). Such high relative increases in length within half a year were generally restricted to mussels of spat size. Only in the first winter of their life, mussels showed positive changes in their mean AFDM (Fig. 2b: the points around the 12th month). At older ages (the points around the 24th and 36th month), mean weights declined between August and March. Consequently, mean production estimates for the entire Balgzand were always positive in the first autumn/winter (in 14 out of 14 observations, P < 0.01, sign test), but mostly negative in later autumn/winter periods of their life (in 17 out of 19 cases for the second winter and in 17 out of 21 cases for the third and later winters, P < 0.01 and <0.05, respectively). The numbers of observations used above were lower than the 26 winters of the total period of observation, because (to include only meaningful estimates) data were excluded if they referred to extremely low numbers of mussels, i.e., if the 15-station mean biomass at the start of the winter half-year was <0.2 g AFDM m⁻².

Fig. 2

Mytilus edulis. Half-year estimates of growth in a shell length (in mm with 1 SE) and b weight of soft parts (in g AFDM) during the first years of the life of mussels on Balgzand. Horizontal axis is a time scale in months, with 0 = January 1st of year of birth (the actual start of benthic life is around June 1st: the 0-points indicated at 5 mo). Examples of three cohorts are shown: those born in 1969 (solid squares), in 1973 (open circles), and in 1987 (solid triangles). The 1969 and 1973 cohorts were sampled on the mussel bed near transects-9 and -10 (see Fig. 1 for location and Fig. 4 for length distributions) and the one born in 1987 was sampled on transect 5 (i.e., the transect with the highest numbers of mussels). Numbers of observation per point >30

Seasonal weight losses of soft parts in bivalves (including mussels) are related to winter temperatures: the colder the winter, the lower the proportion of weight which gets lost (Zwarts 1991; Honkoop and Beukema 1997). As expected, this relationship was found to be also true for the proportions of March-to-August mussel production lost in the subsequent autumn/winter season (Fig. 3). In half-year periods including a mild winter, relative losses in mussel production (in the second winter) were more than double (∼20 vs. <10%) the losses observed in years characterised by a cold winter (r = 0.59, n = 19, P < 0.01).

Fig. 3

Mytilus edulis. Relationship between mean winter (Dec, Jan, Feb) water temperature (T, in °C) and proportion (pr, in %) of second-growing-season (spring/summer) production lost in subsequent part of the year (till March). Linear fit: pr = −2.8 + 3.7 T (n = 19, r = +0.59, P < 0.01). Years (19 out of 26) were included only if the spring/summer production had been >0.1 g AFDM m⁻²

The regular alternation of seasonal periods of weight gain and weight loss complicates an unbiased estimation of annual net production, as was more fully explained in Beukema and Dekker (2006). In short: production values based exclusively on late-summer samples would grossly overestimate net annual production, whereas values based exclusively on late-winter samples would result in underestimates. Therefore, annual net production values were calculated as sums of two successive seasonal estimates: an invariably positive spring/summer one (March–August periods) and a frequently negative autumn/winter one (August–March periods).

Production of an intertidal mussel bed

The studied mussel bed came into being after a heavy settlement in the summer of 1969. By the end of this summer, some ten million spat-sized mussels were present in an area of ∼1,100 m². In the subsequent 5 years, densities of new recruits were lower by one or two orders of magnitude compared to the strong year class born in 1969 (Table 1). The strong 1969-cohort numerically dominated the bed for several years (Fig. 4). Only by 1974, the cohort born in 1973 became more numerous. Mean individual weights increased rapidly in the first three growing seasons, but at older ages growth rates slowed down (Table 1).

Table 1 Mytilus edulis. Annual production (in kg AFDM m⁻² a⁻¹) by six separate cohorts (yearclasses 1969–1974, incl.) in a 1,100-m² mussel bed: before August 1969, between August 1969 and September 1970, between September 1970 and September 1971, etc., as estimated from mean individual weights (g AFDM) and numerical densities (10³ m^–2) recorded annually in late summer or early autumn (A August, S September, and O October)

Fig. 4

Mytilus edulis. Size frequency distributions (Percentage of total number of measured individuals) of mussels in mixed samples taken at several places within a ∼1,100 m² bed (located near transect-9 and -10, Fig. 1), annually in late summer or autumn. a August 1969 (n = 197, estimated number present ∼10⁷), b September 1970 (n = 489, ∼7 × 10⁶), c September 1971 (n = 118, ∼3 × 10⁶), d October 1972 (n = 245, ∼10⁶), e October 1973 (n = 158, ∼2 × 10⁶), and f September 1974 (n = 372, 10⁶); black columns the strong cohorts of 1969 and 1973

The data presented in Table 1 were used to estimate annual production values for the entire mussel bed. The annual production estimates (Table 1) are expressed in kg AFDM m⁻² a⁻¹ by dividing the annual total-bed production figures by 1,100 m², being the maximal size of the bed. As expected from their numerical dominance, the members of the 1969 cohort dominated the production for several years. The total 6-year (1969–1974, incl.) production amounted to 3.43 kg AFDM m⁻² , thus mean annual production in the bed area was 0.57 kg AFDM m⁻² a⁻¹. The share of the strong cohort born in 1969 over this 6-year period was 93%. It should be noted that the above estimate of 0.57 kg is an overestimate (by some 10 or 20%) of the true annual production as only late summer/early autumn sampling data were used to calculate this value.

Production in a large intertidal area

In the 1979–2005 period, the annual 15-site means of M. edulis production varied strongly from <1 to ∼17 g AFDM m⁻² a⁻¹ with a 26-year average of 5.5 ± 0.8 (SE) g AFDM m⁻² a⁻¹ (Fig. 5a). It makes little sense to add standard errors to the individual annual values, because these figures would reflect the huge place-to-place variability in mussel production on the heterogeneous Balgzand area (including many 0-values). The mean annual values of net production shown in Fig. 5a do not show a statistically significant long-term trend.

Fig. 5

Mytilus edulis. Long-term series of successive estimates for the mussel population on Balgzand (means of all 15 sampling sites): a net annual somatic production of soft parts P expressed in g AFDM m⁻² a⁻¹; b mean annual biomass B expressed in g AFDM m⁻² ; c annual P/B ratio expressed in a⁻¹. Year of growing season (Y) is indicated on horizontal axis. Only the long-term trend in biomass was significantly different from a horizontal line: B = 217 − 0.108 Y, r = −0.45, n = 26, P < 0.05

Mussel biomass on Balgzand was highly variable, mean annual values ranging from <1 to >8 g AFDM m⁻² (Fig. 5b), with a 26-years average of 3.2 ± 0.4 g AFDM m⁻². The long-term trend was significantly negative, with a mean annual rate of change of −0.11 g a⁻¹ r = −0.45, n = 26, P = 0.02). Though the decline in annual production was non-significant (Fig. 5a: r = −0.30, P = 0.13), the lowest values of annual production were found in the second half of the period of observation.

Mussels were most productive in the second spring/summer half-year of their life (Fig. 6a). In their first summer, individual weight gain was too low and in the third and later spring/summer seasons their numbers had become too low for high production values. In the autumn/winter half-years, growth and production were positive only in the first winter of life, whereas in the second and later winters mussels generally lost weight and as a consequence contributions to production of these age groups together were usually (in 22 of the 26 winters) negative. On average, the negative contributions of >1-year-old mussels amounted to 14% of the foregoing spring/summer production (without a clear relationship with age). In Fig. 6a, the size of the one positive white column is similar to the sum of the sizes of the subsequent negative white columns. This means that, on average, the positive autumn/winter contributions by juveniles roughly compensates the negative ones by adults. As a consequence, the net annual production values summed for all age groups together were close to the summed spring/summer production values (5.5 and 5.6 g AFDM m⁻² , respectively).

Fig. 6

Mytilus edulis. Relationship between age (year-1 starts in spring of year of birth) and parameters of production and biomass. a Mean half-year production (in g AFDM m⁻², with 1 SE, n = 26) in the 1979–2005 period on Balgzand (means of all sampling sites). Black columns refer to production completed within the main growing season (March–August), white ones to the remaining part of the year. Production during the first growing season (left-most black column) was not actually measured, but arbitrarily set at two times the August spat biomass. b Mean biomass (in g AFDW m⁻², with 1 SE, n = 26). Black columns refer to late-summer samples and white ones to late-winter samples. c Mean P/B ratio values (in a⁻¹, with 1 SE, numbers of observations indicated in brackets). No value included for the first year of life, because of the arbitrary nature of the first-growing season estimate of P. Double asterisks denotes statistical significance (P < 0.01)

Biomass values of individual cohorts were maximal at the end of their second summer (at an age of almost one-and-a-half year) and subsequently declined (Fig. 6b): apparently, losses by mortality became larger than production by growth. Late-winter biomass values were always lower than preceding summer values as a consequence of both individual weight loss and population losses by mortality. Summer biomass values were on average larger than preceding late-winter values (Fig. 6b) due to spring/summer growth exceeding losses by mortality.

With increasing individual weights at aging, annual weight gains comprise a declining proportion of mean individual weight and therefore the ratio P/B is bound to decline with age. In mussels, these values declined from >2 a⁻¹ in the second year of life to about 0.2 a⁻¹ in the fifth and later years of life (Fig. 6c). Only the declines from the 2nd to each of the later years proved to be statistically significant (Wilcoxon tests, P < 0.01). Annual P/B values for the entire population are thus bound to depend on age composition of the population. These P/B values varied strongly from ∼0.5 to ∼3 a⁻¹ (Fig. 5c), with a long-term average of 1.7 ± 0.2 a⁻¹ (n = 26). Particularly high values were found when the population was composed almost entirely of young animals, such as in 1979 + 1980, 1988, and 1991 due to exceptionally strong new cohorts born in 1979, 1987, and 1991, respectively. The extremely high P/B value of 1990 had a different cause: intensive fishing just before August had removed nearly all mussels from a then large stock on Balgzand, resulting in a high value of P but a low one of B (in August B was close to 0 instead of >10 g m⁻² if there had been no fishing). The low P/B values around 2001 were caused by failing recruitment during the entire 1998–2002 period, leading to populations of almost exclusively old individuals in the years 2000–2002, incl.

Production in any year was usually made up of contributions from several cohorts. Cohort production (over their subsequent life time) proved to be strongly related to their numerical density early in life (Fig. 7). The first estimates available of recruit densities were obtained in August, at an age of ∼3 months, and the numbers at that time only loosely predicted subsequent net production (Fig. 7a), though the correlation was significantly positive, with a mean production per recruit of 0.04 g AFDM. This amount is much less than maximal weights of soft parts of adults (1–2 g). Early death prevented substantial production in by far most of the individual recruits. The second estimate of recruit density was made in March, at an age of ∼10 months, when the recruits were still of spat size, but had survived their first winter. The by then much lower numbers of recruits (an average of only ∼15% of the August spat numbers were still present on Balgzand in March) predicted subsequent production much more precisely with an average amount per individual of 0.21 g AFDM (Fig. 7b). Recruit densities (as assessed in March) were closely correlated with subsequent life-time production of the cohort and explained >90% of its variability (regression analysis).

Fig. 7

Mytilus edulis. Dependence of subsequent life-time net cohort production of soft parts (P, g AFDM m⁻²) on initial recruit numbers (D, n m⁻²), as assessed either in August at an age of ∼3 mo or in March at an age of ∼10 mo. Data of 24 cohorts (born in 1979–2002, incl.). a August recruits. Linear fit: P = 1.84 + 0.03 D (n = 24, r = 0.50, P < 0.05) or P = 0.04 D if forced through the origin. b March recruits. Linear fit: P = 0.1 + 0.21 D (n = 24, r = 0.97, P < 0.0001) or P = 0.211 D if forced through the origin. c March recruits. Polynomal fit: P = −0.04 + 0.226 D − 0.00012 D² (n = 24, r = 0.95, P < 0.0001)

A curvilinear (polynomal) best fit (Fig. 7c) of the same data did not yield a closer relationship (r = 0.95 for the curvilinear vs. 0.97 for the straight line relation). In fact, the two lines were hardly different (compare Figs. 7b, c). The close and approximately linear relationships between recruit densities and subsequent life-time cohort production values suggest a rather low and non-consistent variability in mean production-per-recruit values of the cohorts. Production per March recruit did not significantly change with recruit density (r = −0.15, n = 23, P = 0.5), nor did production per August recruit (r = −0.30, P = 0.2). Thus, production per recruit was not significantly lower in strong than in weak year classes, and cohort production showed a close to linear increase with early cohort strength.

Local variation in production

Long-term mean annual production values varied significantly (1-way ANOVA, P < 0.0001) between the individual Balgzand sampling sites. Sites with high and low production were not randomly scattered over the area. Low production values were found particularly at the higher intertidal levels, i.e., at the seven sampling sites with a height higher than ∼4 dm below MTL (i.e., about half-way between the MTL and the low-water level). Production was also low at the two sites close to the low-water level (Fig. 8a, compare Fig. 1). Thus, only the central part of Balgzand (between ∼2 and ∼8 km from the SW coast with heights between ∼4 and ∼7 dm below MTL) showed high mussel production.

Fig. 8

Mytilus edulis. Relationships between intertidal height at each of the 15 sampling sites in dm below or above mean-tide level (MTL) and parameters of production and recruitment. a Long-term mean somatic production of soft parts in g AFDM m⁻² a⁻¹ with 1 SE (n = 26, the production in the years 1979–2004); numbers or letters indicate the different sampling sites (see Fig. 1). b Long-term mean density of recruits of the 26 cohorts born in 1979 to 2004 as observed in August at an age of ∼3 months (open symbols) and in March at an age of ∼10 months (solid symbols), in n m⁻² with 1 SE (n = 26). c Mean life-time production per March-recruit in g AFDM ind⁻¹ with 1 SE, as realised after the first (month of) March of their life (i.e., excluding the production at spat size). Points for transects 1, 2 and 3 (see Fig. 1) omitted, because of extremely low (on average < 0.1 m⁻²) recruit densities. Means of all cohorts born in the 1979–2002 period, thus n = 24 or less (because of cohorts with 0-recruit numbers at the site in March); numbers indicate numbers of observations

Mean annual recruitment (26-years mean numbers m⁻² of spat-size mussels in August or March) showed similar relationships with intertidal level (Fig. 8b). Spatial variation in mean local recruit densities explained most of the observed spatial variation in mean production. For August recruitment RAug and total mean production P of the mussels of all size groups, the relationship was P = −0.2 + 0.074 RAug, with r = 0.87 (n = 15, P < 0.001, best linear fit) and for recruit densities in March, RMa, production P after March (thus of adult-size mussels only) it was P = 0.75 + 0.27 RMa, with r = 0.98 (n = 15, P < 0.0001).

Though place-to-place differences in mean annual recruitment could explain almost the complete variation in mean annual production between the sampling sites, mean individual life-time production per (March) recruit also varied significantly (1-way ANOVA, P < 0.001) between sampling sites (Fig. 8c). In lower intertidal areas higher mean values in production per individual recruit were observed than in high areas. The negative relationship of production per recruit (P/R) with intertidal height (H) was statistically highly significant: P/R = 0.04 − 0.029 H (r = −0.86, n = 12, P < 0.001, best linear fit).

Discussion

Methodological comments

Reliable and precise estimates of annual production can be obtained only if the basic data on numerical densities and individual weights are sufficiently reliable and precise. Mussels are notoriously difficult to sample on a large scale as a consequence of the way they are distributed. The spatial distribution pattern of by far most mussels on tidal flats is not more or less randomly over large areas but highly clustered in beds of limited size, covering usually only one or a few percent of the total area. As a consequence, coefficients of variation for local densities are extremely high (Beukema et al. 1983). To obtain precise estimates of the size of mussel beds and the numbers present within the beds, special sampling methods are desirable, e.g. a combination of aerial photographs and ground sampling (Nehls 1999; Brinkman et al. 2003; Herlyn 2005). We experienced that such time-consuming procedures were not feasible within the framework of a large-scale and long-term sampling program of the entire benthic fauna. Therefore, we had to rely on the existing data gathered in a standardized sampling program that was aimed at a study of the total macrozoobenthos and was certainly not optimal for sampling mussels. Fortunately, the fixed transects of this program crossed some extensive and more or less permanent mussel beds, enabling long-term comparisons of production estimates without too much sampling error.

The question whether or not the hardly adequate method nevertheless resulted in representative estimates of mussel abundance can be answered by a comparison with the results of two much more extensive sampling programs on Balgzand executed by the Netherlands Institute for Fishery Research (RIVO). The data shown in Fig. 9 were kindly made available by P. Kamermans (personal communication, 2006) and the RIVO methods are described in Steenbergen et al. (2005). On average, the results of the 8-years RIVO data set (1998–2005) of a program specifically directed to an assessment of mussel stocks (involving measurements of the size of the beds) closely corresponded with our (NIOZ) estimates: RIVO = 1.07 NIOZ (r = 0.69, n = 8, P < 0.05). The comparison with the longer, but not mussel-directed, RIVO data set (1990–2005) from a sampling program designed for the assessment of cockle stocks showed a rather loose, though highly significant, correspondence: RIVO = 2.08 NIOZ (r = 0.92, P < 0.001). The RIVO and NIOZ data thus showed a similar year-to-year variability. The on average higher RIVO values arose from a few substantially higher annual values, but the number of years with higher RIVO values (9 out of 16) was not significantly different from a 50% expected by chance.

Fig. 9

Mytilus edulis. A comparison between results of independent mussel stock assessments on Balgzand by NIOZ and RIVO. Horizontal axis: 15-station means by the authors (NIOZ) in g AFDM m⁻², samples taken in March of each year. Vertical axis: (solid symbols) means of ∼70 stations sampled by RIVO in spring of 16 years (1990–2005), and (open symbols) results of specific mussel bed assessments by RIVO in spring of 8 years (1998–2005). The RIVO data were given in whole-mussel weights and converted to AFDM by multiplication by 0.07 (assuming conversion factors of 0.28 from total to flesh weight, and 0.25 from fresh flesh to ash-free dry weight). The line indicates exactly equal estimates

As a consequence of the observed pattern of seasonal changes in weight, including a season of weight loss in adult mussels (Fig. 2b), once-annual sampling could easily result in seriously biased production estimates (compare Fig. 3b of Beukema and Dekker 2006). Unbiased estimates can be obtained only by a proper choice of the annual sampling period. As explained by Beukema and Dekker (2006), the best procedure in such cases is to sample twice-annually, specifically within the two annual periods of relatively slow weight change, i.e., around the times of minimal and maximal individual weights. In Wadden Sea mussels, the available data point to optimal sampling periods in late winter and in late summer: Pieters et al. (1979) found minimal values for dry tissue weights in March and maximal ones in October in a mussel bed in the Wadden Sea, Dankers et al. (1989) and Dankers and Zuidema (1995) followed mean weights of soft parts of mussels at several places in the western Wadden Sea and usually found minimal weights in winter (mostly February) and maximal ones in late summer (August or September). Thus, our February/March and August/September samplings appear to be fairly adequately timed for a twice-annual sampling program of (adult) mussels. In some years, weight gain may have continued for one or two months after August, but model calculations (not shown) revealed that the resulting underestimates by a too-early August sampling was <10%. Craeymeersch et al. (1986) arrived at a similar conclusion.

The above pattern of alternating half-year periods of weight gain and weight loss appear to be general for (adult) bivalves in temperate areas. They have also been observed in mussels in areas outside the Wadden Sea, e.g., in British waters (Dare and Edwards 1975; Dare 1976; Bayne and Worrall 1980). Other bivalves show similar patterns; examples from the Wadden Sea include cockle C. edule (Zwarts 1991; Beukema and Dekker 2006), sand gaper M. arenaria (Zwarts 1991), Pacific oyster Crassostrea gigas (Cardoso et al. 2007) and the tellinids Scrobicularia plana (Zwarts 1991) and M. balthica (Beukema et al. 1985; Zwarts 1991; Dekker and Beukema 2007). Honkoop and Beukema (1997) experimentally showed that food availability and temperature determine whether individual weights increase or decrease. The balance between food intake and energy need becomes negative in periods when food supply is low and temperatures are high. This is particularly the case in mild winters, when weight losses were found to be relatively high (as compared to cold winters) in all bivalve species studied by Honkoop and Beukema (1997): M. edulis, C. edule, and M. balthica. As a result, higher proportions of production will get lost in mild than in cold winters, as shown in Fig. 3. In mussels, autumn/winter losses appear to be low (Figs 3, 6) as compared to cockles C. edule (Beukema and Dekker 2006) and M. balthica (Dekker and Beukema 2007). Negative net annual production values were not observed on Balgzand, but can occur elsewhere (Cusson and Bourget 2005).

Long-term variability in production

Studies of a considerable (>10 years) length on production of zoobenthos, enabling some insight into temporal variability and its underlying processes, appear to be rare. So far, they were apparently not available in M. edulis for other areas than Balgzand (Table 2). The studies in these other areas lasted in all cases between <1 and 3 years. In view of the substantial between-year variability in mussel recruitment (Beukema et al. 2001) and production (Fig. 5a), short-term studies are probably inadequate to obtain representative values of this species’ annual production or P/B ratio in any area.

Table 2 Mytilus edulis. Published estimates of mean annual somatic soft-parts production (P, in kg AFDM m⁻² a⁻¹) with estimates of annual mean biomass (B, in kg AFDM m⁻²) and P/B (a⁻¹). (a) Within natural intertidal mussel beds of stated size (in 10⁴ m²), and (b) for more extensive coastal areas including mussel beds and/or scattered mussels. Length of observation period indicated (years). To facilitate comparisons, some conversion factors were applied (1 g DM = 0.9 g AFDM, 1 g C = 2 g AFDM)

The long-term mean Balgzand somatic-production estimate of 5.5 g AFDM m⁻² a⁻¹ is well within the range of estimates obtained in other extensive areas (Table 2b: a range of 3–42 with a 6-area average of 15 g AFDM m⁻² a⁻¹). Production per unit area observed within the restricted areas of mussel beds are two orders of magnitude higher: the 6-area average of the estimates summarized in Table 2a is ∼750 g AFDM m⁻² a⁻¹ and the 6-years mean for the studied Balgzand bed amounted to a similar ∼570 g AFDM m⁻² a⁻¹ _. This ∼100-fold difference between large-area and restricted-bed averages agrees with the observation that mussel beds usually cover in the order of 1% of extensive tidal flat areas. In the 1995–2005 period, Steenbergen et al. (2005) found total mussel bed areas to range between 0.01 and 0.9 (mean 0.4) km² on the ∼50 km² of Balgzand tidal flats, i.e., on average 0.8% of the total intertidal area (detailed data were kindly made available by P. Kamermans of the RIVO). For the entire Dutch Wadden Sea, Ens et al. (2004: their Fig. 57b) compiled recent (1990–2004) data on the surface areas of mussel beds and found a range of <1 to ∼25 km², with a 15-years (1990–2004) mean of ∼8 km². This equals ∼0.7% of the total area of intertidal flats. For the ∼5 km² Königshafen area (on the Island of Sylt, German Wadden Sea), Nehls et al. (1997) state 2.5%. For the total of four political parts of the international Wadden Sea (Netherlands, Lower-Saxony, Schleswig-Holstein, and Denmark), each covering around 10³ km² of tidal flats, Dankers et al. (2001) list total areas of mussel beds. Before 1980, mussel beds usually covered 1–3% of the tidal flats, but these proportions declined to ∼0.5% (or even less in some years) in later years. The declines of intertidal mussel stocks were particularly strong (to near-extinction) in the western half of the Dutch Wadden Sea (see below), but were also substantial in extensive parts of the German Wadden Sea (Zens et al. 1997; Herlyn and Millat 2000; Nehls et al. 2006).

In particular intensive fishing and failing recruitment are held responsible for the recent declines of Wadden Sea mussel stocks (Ens et al. 2004). Before 1988, recruitment tended to be above-normal once per four years (viz. eight times between 1955 and 1988), but such successful recruitment years never occurred between 1988 and 2001 (data by Van Stralen, cited in fig. 56 of Ens et al. 2004). On Balgzand, mussel recruitment showed a similar pattern (see Fig. 1c of Beukema et al. 2001). Despite recent (after 1991) restrictions on mussel fishery in all intertidal parts of the Wadden Sea, recovery of mussel beds appears to be a slow process. Particularly in the western half of the Dutch Wadden Sea, few of the beds have returned up to now after their removal by intensive fishery in 1990; compare recent maps in Steenbergen et al. (2005) and Goudswaard et al. (2006) with older maps in Dijkema (1989). A lack of suitable hard substrates (i.e., existing mussel beds) may have contributed to the slowness of the recovery (compare Dare et al. 2004).

The long-term mean Balgzand production estimate of 5.5 g AFDM m⁻² a⁻¹ does not fully account for the total annual production by M. edulis, because it does not include gamete production nor the organic contents of the shells. Mussel shells contain high proportions of organic material and the amounts of shell organic material may almost equal those of the soft body (Gardner and Thomas 1987). A good guess for shell organic-matter production might be 75% of the spring/summer production of the soft parts, i.e., on average ∼4 g m⁻² a⁻¹. Gamete production in mussels is highly variable, depending a.o. on their condition and size/age (Seed and Suchanek 1992). A fair guess might be 25% of the mean biomass of adult mussels (Sprung 1983), i.e., ∼0.7 g AFDM m⁻² a⁻¹. Thus, total annual net production by mussels on Balgzand might have amounted to a long-term average of around 10 g organic material m⁻² a⁻¹. This value equals the similarly obtained one in the cockle C. edule (Beukema and Dekker 2006) and might account for almost 20% of the total macrozoobenthic production on Balgzand.

Local variability of production

The estimates of mean annual M. edulis production at the 15 Balgzand sites varied strongly, ranging from close to 0 to tens of g AFDM m⁻² a⁻¹ (Fig. 8a). They reflected more or less the local variation in annual recruitment, with maximal values half-way between the low-water level and MTL, and very low values at more extreme levels (Fig. 8b). The near-absence of mussels at very low and very high levels is in accordance with the results of a habitat-suitability study for mussel beds by Brinkman et al. (2002). Spatial variation in mean local recruit densities explained a high proportion of the observed spatial variation in mean production. Nevertheless, production per recruit was not the same at all intertidal levels, it tended to be higher at low than at high intertidal levels (Fig. 8c). Similar relationships with intertidal height were observed in other bivalves (Cerastoderma: Beukema and Dekker 2006; Macoma: Dekker and Beukema 2007). As in other bivalve species, mussel growth rates appear to be limited by short daily immersion times at the higher intertidal levels (Wanink and Zwarts 1993; Faldborg et al. 1994; Buschbaum and Saier 2001; Thieltges 2006). At the lowest levels, low annual recruitment in particular appeared to limit productivity (Fig. 8b). This may be caused by high mortality in spat in these areas, caused by epibenthic predators such as shrimps, as has been observed in some other bivalve species on Balgzand (Beukema and Dekker 2005).

Factors underlying production variability

By far most of the variability of life-time production of M. edulis cohorts on Balgzand could be explained by their initial numerical abundance. In particular recruit numbers assessed at an age of ∼10 months (Fig. 7b) explained >90% of the variability of the subsequent life-time somatic production. The same was true for the spatial variability in production (compare Fig. 8a with solid points in Fig. 8b).

As annual recruitment success was highly variable (Beukema et al. 2001), strength of mussel cohorts strongly differed and consequently their contribution to production (Fig. 7). The three most productive cohorts together accounted for more than half of the summed production of all 24 cohorts (1979–2002, incl.) that had (nearly) completed their life-time production in 2006. On the other hand, the 12 cohorts with the lowest production contributed together only ∼5% of the 24-cohort sum. This high between-cohort variability resulted in a high coefficient of variation (100 standard deviation of individual values divided by their mean) of cohort lifetime production, being no less than 145%. This is much higher than the coefficient of variation for annual production estimates of the entire population (all cohorts present together): 65%. The explanation for this lower variability will be that in most years several cohorts were simultaneously present (compare Fig. 4) and weak contributions by some cohorts were usually counterbalanced by a substantial one by a strong cohort. Such strong cohorts can dominate production for several years (compare in Table 1 the contribution of the 1969-cohort with the total annual production).

Initial numbers (recruit densities in March) predicted precisely subsequent (lifetime) production of a cohort (Fig. 7b: r = 0.97). Though the prediction of annual production of the entire population by initial (total: spat + adult) numbers in March was less precise (r = 0.64), the correlation found was still highly significant (P < 0.01, n = 26). Thus again variability of numbers appeared to be a decisive factor for production. On the other hand, production per recruit did not show significant relationships with numerical density. Apparently, there was no appreciable reduction in the production per recruit at high densities by food limitation: within the range of densities occurring on Balgzand we found no evidence for density-dependent growth rates resulting in density-dependent production per recruit. Only spatial variation in production-per-recruit values varied significantly (Fig. 8c), showing a significant relationship with intertidal height. Probably, feeding time was limiting growth and thus production per recruit at high intertidal levels.

In conclusion: the year-to-year variability can only be understood by insight into recruitment variability. Variability in annual mussel recruitment success proved to be highly synchronized with that in other bivalve species in the Wadden Sea (Beukema et al. 2001). In all major bivalve species of the Waddden Sea, temperature of the foregoing winter was significantly negatively correlated with recruit abundance in summer, as also shown by Strasser et al. (2003). Maybe, a climatic factor such as winter character does not directly affect bivalve recruitment, but exerts its influence by the action of such biotic factor as predation pressure on early life stages by epibenthic predators such as shrimps and crabs. Evidence for such mechanism has been presented in mussels in the Wash by Young et al. (1996) and in the bivalves M. balthica, C. edule and M. arenaria by Strasser (2002), Philippart et al. (2003), and Beukema and Dekker (2005). The slow recovery of mussel stocks observed in various parts of the international Wadden Sea in recent years (roughly after the overfishing in 1990) may have to do with the declining frequency of cold winters (Strasser et al. 2003; Nehls et al. 2006), though a lack of suitable hard substrate after the almost complete removal of mussel beds in 1990 cannot be excluded as an additional cause (Dare et al. 2004). If climate factors were the main factor governing year class strength in Wadden Sea bivalves, the declining trend of bivalve stocks (including those of the mussel: Fig. 5b) may be expected to continue as long as the expected global warming continues.