Introduction

Birth and death processes are at the core of population ecology and demography. Obviously, most attention has been paid to the birth and death rates of individual organisms. Sometimes, however, it might be more convenient to work at a higher level of biological organization (Boughton and Malvadkar 2002). Students of colonial insects such as ants and termites often focus on the birth and death of colonies. The life time of an ant colony may coincide with that of the queen, in which case one might argue that work is still done at the level of the individual organism. But, as queen supersedure within a colony occurs in many species, such is not necessarily true. In plant ecology, Watt’s seminal paper on the dynamics of plant communities resulted in the idea of gap-phase dynamics, which still plays a central role in forest ecology (Watt 1947; Gratzer and others 2004; Hunter and others 2015). In forests, gaps appear as a result of, for example, storms or of the fall of a single canopy individual. The regeneration rate of the gap is then related to the size of the gap and characteristics of the surroundings (van der Maarel 1996). Similar thoughts are developed in coastal ecology, where hard substrata are discontinuous pieces of habitat surrounded by sand and mud. These hard substrata can almost entirely be covered by epifaunal invertebrates such as sponges, tunicates and mussels. Clearings are regularly made either by physical forces or by predators. The recolonization of such clearings will then, just as in the case of forest gaps, depend upon size and the characteristics of the surroundings (Connell and Keough 1985). Eventually these ideas evolved to a more general theory of patch dynamics, where the birth, growth or shrinkage, and death of patches rather than of individual organisms are the central issue (Levin and others 1993; Levin and Paine 1974; Picket and White 1985). Metapopulation ecology, introduced by Hanski, makes a further step by using entire populations, connected into a metapopulation, as the main unit of interest (Hanski and Gilpin 1991). That is, the birth and death of connected populations are studied.

It is interesting to note that the choice of Watt was not so much based on theoretical, but merely on practical considerations. He acknowledges that ‘the ultimate parts of the community are the individual plants,’ but adds ’but a description of it in terms of the characters of these units … is impracticable …’ (Watt 1947). In fact, the same argument is used by students of colonial insects, subtidal epifauna and island populations. It is much easier and it requires less data to describe the fate of patches, communities or populations than that of individual organisms.

Our interest is in the dynamics of intertidal bed-forming bivalve populations living on the soft sediments of the Dutch Wadden Sea, such as the blue mussel Mytilus edulis and the recently introduced Pacific oyster Crassostrea gigas. These intertidal bivalve beds are the habitat of many benthic and epibenthic species and attract numerous birds and fishes, for which these beds provide a rich food source (Goss-Custard and others 1982; Waser and others 2016). As such, they have a high conservation value, and conservationists were worried to see that in the late 1980s, in a period with low spatfall and ongoing fisheries on mussels on these beds, hardly any intertidal bed remained (Ens 2006; Herlyn and Millat 2000). Since then, these beds have partly recovered, but doubts remain as to what level compared to the pristine situation (Dankers and others 2001). It has even been claimed that due to sediment disturbance by mechanical fishing activities, the Wadden Sea has undergone drastic changes and collapsed into an alternative stable state, where recovery of bivalve beds is basically impossible without large-scale restoration projects (Eriksson and others 2010). Recently, quite some effort has been put into initiating such restoration projects, but without much success (van der Heide and others 2014; de Paoli and others 2015).

The aim of our study is to investigate if large-scale restoration projects are indeed needed for the recovery of bivalve beds. To this end, we estimate present recruitment and survival rates of intertidal bivalve beds, predict population size and fluctuations therein and compare that with available recent and historical data, using the bivalve bed as the study unit, instead of the individual organisms. It must be noted though that in one respect, the system cannot return to its original state, due to the invasion of the Pacific oyster (Troost 2010). Although the Pacific oyster was introduced in the Dutch Wadden Sea around 1978, it did not really spread until the late 1990 s, that is, well after the disappearance of the intertidal mussel beds (Reise 1998; Diederich 2005; Nehls and others 2006; Troost 2010). Whereas in the past bivalve beds consisted entirely of mussels, nowadays bivalve beds may consist entirely of mussels or oysters, or a mixture of these two species. Apart from studying present-day recruitment and survival of bivalve beds in general, we also estimate these parameters for different types of bivalve bed.

To summarize, we examine whether present recruitment and survival rates of mussel–oyster beds are sufficient to achieve and maintain historic population levels.

Materials and Methods

Definition of a Bivalve Bed

In 2002, Dutch, German and Danish scientists agreed upon a common definition of littoral mussel beds to facilitate trilateral comparisons (CWSS 2002; Herlyn 2005; Nehls and others 2009). A mussel bed was defined as a benthic community structured by blue mussels that may consist of an irregular collection of more or less protruding smaller patches, separated by open spaces. Patches should be larger than 1 m in diameter and within 25 m distance from each other. A collection of humps smaller than 1 m in diameter should have an areal coverage of more than 5% to be grouped within one patch (Figure 1). Following the expansion of the Pacific oyster into many mussel beds of the Wadden Sea, the original definition gained a more general character and is now used to define different types of bivalve beds with varying amounts of mussels and oysters. In the Netherlands, beds are defined as pure mussel beds if the cover of Pacific oysters is less than 5% and as pure oyster beds if the cover of mussels is below 5%. In all other cases, beds are denoted as mixed beds.

Figure 1
figure 1

Visual representation of the demarcation of mussel beds [from the Wadden Sea Quality Status Report 2009 (Nehls and others 2009)].

Full size image

We follow this definition, but because it does not say anything about the time component, we added the rule that patches from different years that overlap in space belong to the same bed. If two patches from the same year spatially overlap with the same bed in the previous year or in the following year, then these two patches are considered to belong to the same bed, even if they are further away from each other than 25 meter (Figure 2). Within this definition, beds can neither split nor merge. They can only be born, survive or die. The definition requires an iterative procedure to uniquely classify all patches in beds.

Figure 2
figure 2

Changes over time of mussel bed 1999.048. Each color represents a different year. Note that some green patches are only grouped within the same bed, as they overlap with the same blue bed, that in itself consists of patches that only form one bed because they overlap with the same red bed. The thick red line indicates a distance of 50 m (Color figure online).

Full size image

Field Data

Each spring, Wageningen Marine Research (formerly IMARES) and the private company MarinX, map all the littoral mussel and oyster beds in the Dutch Wadden Sea. New bed locations and the presence or absence of already known beds are determined by an inspection of the area from an airplane flying at an altitude of 500 m. Subsequently, researchers walk around the beds during low tide to demarcate the contours, and their tracks are recorded by GPS and imported into GIS. Characteristics of the beds, such as percentage cover of mussels and oysters, are recorded. Due to time constraints, not all beds can be visited each year on the ground. For unvisited beds, data from the previous year are used. Occasionally, it happened that a bed was only discovered in its second year (which can be derived from the size distribution of the individual mussels), in which case the data set of the previous year is updated to include that bed. Bed distribution data for the period 1999 – 2013 are used. Details of the sampling procedure are in Steenbergen and others (2006).

Survival Analysis

Mussels and oysters spawn in late spring and early summer and larval settlement occurs a few weeks after spawning. These beds become visible to the human eye in autumn. The survey is each spring, and the first appearance of a bed in the study (‘recruitment’) thus occurs almost a year after its ‘birth.’ We assume that the beds disappear (‘death’) shortly after their last appearance in the study. More precisely, we assume that beds have a lifetime in years equal to the number of times they occurred in the study. However, such assumption could not be made for beds that were still present in 2013, the last year of the study. We only know that they lived for at least a specific number of years. For example, when a bed first appears in the study in 2010 and is still present in 2013, we only know that its lifetime is at least 4 years. A survival analysis using Cox’s proportional hazard model was performed (Klein and Moeschberger 2003), as this analysis, which estimates an age-dependent survival curve, is capable of dealing with such so-called right-censored data. For beds that were already present in 1999, which is the first year of our study, the year of birth and thus the age are unknown and cannot easily be imputed either. These beds were therefore necessarily left out of the analysis. Each year, each bed was classified as either mussel, oyster or mixed. Bed type was therefore included in the model as a time-dependent qualitative covariate. The proportional hazard assumption of Cox’ model implies that the hazard, which is roughly speaking the rate of disappearing, of, for example, a mussel bed of a specific age, is proportional to the hazard of an oyster bed of the same age. The model further contained various quantitative and time-independent covariates (log bed size, orbital speed and duration of exposure of the midpoint of each bed), which means that these variables are assumed not to change over the study period. For orbital speed and duration of exposure, variables that do not change much over time, there was but a single measurement available. For longitude and bed size, annual measurements were available, but they hardly changed over time, as, for example, indicated by the high Pearson’s correlation coefficient between the initial value and the average over the entire period (r = 1.000 and r = 0.978, respectively). For both these variables, the initial value was used in the analysis. All quantitative covariates are standardized, which means that their mean becomes zero and their standard deviation one. This allows for an easier interpretation and comparison of the estimated regression parameters.

Depth data were derived from a bathymetric grid at a resolution of 20 by 20 m throughout the Dutch Wadden Sea. This bathymetry is based mostly on LIDAR soundings performed by Rijkswaterstaat between 2006 and 2012. For details on this so-called cycle 5 data set, see Elias and Wang (2013). Duration of exposure was obtained by combining this bathymetric grid with water level data of 15 tidal gauge stations encompassing the Dutch Wadden Sea. At all stations the sea level height is recorded every ten minutes. Any combination of three stations defines a triangular plane, which allows interpolation of the water level at all locations within the triangle. For details, see Rappoldt and others (2014). A location at a certain time is considered inundated if the estimated water level is above its height on the bathymetric grid. By calculating the inundation over the period 1999–2013 for all bathymetric grid points, we obtained a relative exposure time. Wave orbital speed was modeled using SWAN, a 2D horizontal wave model (Booij and others 1999), also based on the cycle 5 bathymetric grid. A total of 1480 different environmental scenarios based on 3 varying parameters (water level, wind direction and wind speed) were considered. For the period 1991–2013, wind (source KNMI) and water level data (source Rijkswaterstaat) at 1 h intervals were used to create time series of wind direction (8 bins: N, NE, E, SE, S, SW, W and NW), wind speed (5 bins: 2–6, 6–10, 10–14, 14–18 and 18–22 m/s) and water levels (37 bins: from − 0.9 to 2.7 m around mean tide level, in 0.1 m intervals) for each point of the bathymetric grid. Orbital velocities were modeled for all 1480 possible parameter combinations. Subsequently, trimmed mean orbital velocity was determined for each grid point. For details, see Donker (2015).

Simulation Model

We firstly developed a simple simulation model of bed dynamics, where in each year of the simulation one number is randomly selected from the observed 14-year series (2000–2013) of recruitment area data (Figure 3). Kendall’s rank correlation test showed that there was no significant trend over time in recruitment area (τ = −0.30, p = 0.16). Summing the area of all beds that were recruited in the sampled year thus provides the initial area of a simulated cohort. For each cohort and each subsequent year, the fraction of the initial area that survives is simply taken from the estimated survival curves for pure mussel beds and oyster/mixed beds. The population was simulated 1000 times, each time for a period of 30 years. Next we took a more complex approach, again randomly selecting recruitment years, but now also randomly selecting the survival time for each recruited bed separately. The survival time was selected by comparing the predicted survival curve with a pseudo-randomly generated uniformly distributed number between 0 and 1. The model including the three binary variables bed size (larger or smaller than median), exposure time (shorter or longer than medium) and bed type (mussel or oyster/mixed), was used to predict all survival curves. We either assumed that all beds are mussel beds or that all beds are oyster/mixed beds.

Figure 3
figure 3

Surface area of new beds in spring versus year. Blue indicates the area of pure mussel beds (Color figure online).

Full size image

As the survival function could only be estimated up to an age of 14 years, we extrapolated the curve using a constant annual increase in log hazard rate. We estimated this annual increase parameter, assuming a linear increase in log hazard rate over the last 6-year age period (Figure 4).

Figure 4
figure 4

Log cumulative hazard rate for pure mussel beds (upper panel) and oyster/mixed beds (lower panel) versus age. Red line shows the linear relationship fitted by using the last 6 observed values. Blue dots are extrapolated values (Color figure online).

Full size image

Software

All analyses were performed using the R platform (R Core Team 2016). We used the packages sp, maptools, rgeos, rgdal, raster and spdep for spatial analysis and the package survival for survival analysis. Scripts are available from the corresponding author.

Results

The classification of all observed or interpolated bed contours over the period 1999–2013 resulted in a total of 1436 beds, categorized as either mussel, oyster or mixed in each year of their existence (Table 1). A plot of the surface area of each bed in a specific year against the area of the same bed in the following year revealed no clear pattern of growth or shrinkage (Figure 5). Most beds only show relatively small changes in area, either a slight growth or a slight shrinkage.

Table 1 Original Data Can Be Summarized by a 1436 by 15 Matrix Where Each Element Contains the Type of Bed
Full size table
Figure 5
figure 5

Area of individual beds in two succeeding years. Points above the diagonal line point to growth and those below the line to shrinkage of beds. Points exactly on the diagonal are most likely the result of imputed data.

Full size image

The fate of all beds in terms of their year of birth and year of death is summarized in Table 2. The first striking result is that more than half of the beds, that is, 751 out of 1436, are on the main diagonal, which tells that these beds die after about one year. Another 232 and 121 die in the following two years. This result is confirmed by estimating the survival curves for each overall bed type separately (Figure 6). The survival curves, which show the probability of being alive versus age, drop down quickly in the first few years. Yet, after an age of about 4–5 years the survival functions go down very slowly, which implies that the hazard rate (the chance of disappearing) becomes very low. The baseline survival function of Cox’s proportional hazard model confirms this result. The model parameter estimates show the effect of the various standardized covariates (Table 3). For example, for those beds that are one standard deviation larger than the mean bed size, the hazard rate is only 61% of the average hazard rate, and beds that have a one standard deviation higher exposure time, have a 9% higher hazard rate (Table 3). The effect of orbital speed was weak and not significant, the strongest effect being due to bed type. Oyster and mixed beds have a much lower hazard rate and thus a much higher survival than pure mussel beds, the hazard rate for oyster beds is only 30% of that of pure mussel beds, and for mixed beds this figure is even slightly lower 27% (Table 3, Figure 6). The difference in hazard rate between oyster and mixed beds is thus very small, and the Cox model where oyster and mixed beds are pooled in one class was not significantly different.

Table 2 Year of ‘Birth’ (Rows) Versus Year of ‘Death’ (Columns)
Full size table
Figure 6
figure 6

Survival curves of the three different bed types separately. Average plus or minus SE.

Full size image
Table 3 Results of Cox Proportional Hazard Model
Full size table

The proportional hazard assumption was tested for each covariate of the fitted Cox model by correlating the scaled Schoenfeld residuals with time, after deleting the nonsignificant variables longitude and orbital speed and pooling the oyster and mixed classes. For all variables that were included in the model, the test for independence between residuals and time gave a nonsignificant result. The same was true for the global test for the model as a whole (χ2= 5.70, p = 0.13). The proportionality assumption was also visually inspected by plotting the log hazard rate for all eight combinations of Cox’s model, using binary variables for the two continuous variables bed size and exposure time. For both variables, the two classes are given by being smaller or larger than the median. The plots show that all lines are more or less parallel and thus confirm the proportionality assumption test (Figure 7).

Figure 7
figure 7

Log cumulative hazard rates versus age of the different bed types as classified by three binary variables, that is, bed size (larger or smaller than median), exposure time (shorter or longer than medium) and bed type (mussel, blue lines or oyster/mixed, red lines (Color figure online).

Full size image

The expected total bed area follows from multiplying the expected lifespan of a bed, which equals the integral of the (extrapolated) survival function, times the annual average of the sum of all ‘newborn’ bed areas (Figure 3). The expected lifespan is 3.43 years for pure mussel beds and 10.8 years for oyster/mixed beds, and the average total area of newborn beds is 3.07 km. So the expected total bed area is 10.5 km if all beds are mussel beds and 33.2 km if all beds are oyster or mixed beds. The simulations showed that the variability around this expectation is large and has a high serial correlation (Figure 8). There was hardly any difference in terms of mean predicted bed size or 90% prediction intervals between the simple and the more complex modeling approach. Longer runs of 1000 years yielded a 90% prediction interval for mussel beds from 5 to 23 km and for oyster/mixed beds from 21 to 53 km.

Figure 8
figure 8

Simulated total surface area of all bivalve beds for the first 30 years after the complete disappearance in 1990. Blue lines refer to simulations where all beds are assumed to be pure mussel beds, whereas for green they are all pure oyster or mixed beds. Solid lines are examples of a single realization, and dashed lines indicate the 90% prediction intervals. Red line shows the true total bed area over the study period 1999–2013 (Color figure online).

Full size image

The simulations also pointed to a high serial correlation, which is of course due to the large variability in recruitment and the average lifetime of a bed being much larger than one year. It might be helpful to notice that our simple simulation model resembles a moving average process of a large order (but with a nonnormal distribution), and that the simulated series could be well described by a first-order auto-regressive process, sometimes called the Markov process (Chatfield 1989), with a parameter value of α = 0.89, also pointing to a high serial correlation.

Discussion

Using Cox’s proportional hazard model with fixed-time covariates implies that some simplifying assumptions were made. First, the assumption that the hazard rates of two beds are proportional means that the relative risk of disappearing of one bed in comparison with that of another bed does not change with age or time. So the model is not able to accurately describe the possible case that, say, mussel beds will only have a higher risk compared to oyster beds when they are both young, but the same risk when they are both old. Further inspection of the separate survival (Figure 6) and log hazard curves (Figure 7) for the various groups showed, however, that the assumption of a proportional hazard rate is quite reasonable. For the quantitative covariates, the data were split into two groups (above and below the median value) to simplify the inspection. Second, applying fixed-time covariates ignored possible changes over time in the value of the covariates. For example, the initial bed area was used to describe bed size, but bed size may change. But changes over time were relatively small compared to differences among beds. For depth, orbital velocity and inundation, only data from a single survey were used. Annual data are not available. But most likely these variables will not change much over time, for example, shallow areas generally remained shallow over the entire study period (Elias and others 2012). For bed type, annual data were available, and since it appeared that many beds changed type during their lifetime, we therefore used bed type as a time-dependent qualitative covariate.

It is easy to understand that oyster and mussel beds differ in hazard rate, since they do have a totally different physical structure. Mussel beds consist of a dynamic meshwork of individual mussels attached to each other via temporary byssus threads, lacking a permanent anchorage in the substrate. Van de Koppel and co-workers emphasize that this flexibility induces spatial self-organization, thereby increasing resilience against disturbance (van de Koppel and others 2005). Nonetheless, mussel beds are vulnerable to storms and ice scouring (Nehls and Thiel 1993; Donker and others 2015). In contrast, Pacific oysters are permanently attached to each other via an organic-inorganic adhesive. Dead oysters remain attached to each other and oyster larvae preferentially settle on conspecifics, creating rigid and persistent structures (Walles and others 2015). Furthermore, there is also a difference in risk of predation. All size classes of mussels on an intertidal mussel bed are subject to predation by a suite of predators, most notably shorecrabs and shellfish-eating birds (Zwarts and Drent 1981; Smallegange and Van der Meer 2003; van de Kam and others 2004). In contrast, Pacific oysters are only subject to predation when small (Dare and others 1983; Mascaró and Seed 2001; Markert and others 2013; Walles and others 2015). Mixed beds had about the same survival rate as pure oyster beds, and we expect them to become the dominant type of bed in the future. In 2015, mixed beds already predominated in the Wadden Sea, comprising 1152 ha (62%) of a total of 1846 ha bivalve beds (van den Ende and others 2016).

We set out to study survival and recruitment of bivalve beds to investigate the claim that the ecosystem had collapsed after overfishing in the late 1980s and that beds could only recover with the help of large-scale restoration projects. To see whether bed areas are nowadays indeed much lower or whether bed dynamics are very different from the past, we must compare our results to historical information.

The simulation study showed that the large observed variability in annual recruitment of new beds produces individual trajectories that can differ strongly (Figure 8). Due to strong serial correlation, periods of several decades of bed area far above or far below average bed area are no exception. So the question is, whether historical data fall outside these wide prediction bands. Unfortunately, even in the well-studied Dutch Wadden Sea, data of total bed area before 1990 are uncertain and historic data spanning several decades are entirely lacking. In an extensive overview of all data reported in governmental reports, Dankers and others (2003) conclude that in the 1970s and 1980s between 10 and 56 km of bed area was present. A slightly different approach by the same authors arrived at a somewhat smaller, but still large interval of 17–48 km (Dankers and others 2003). The indicated uncertainty is not so much due to interannual variability, but merely to the use of different methods, such as aerial or ground sampling, and different ways of demarcating beds. Our approach yielded a 90% prediction interval of 21–53 km for oyster/mixed beds, which is entirely due to interannual variability. Other long-term studies also point to a huge variability and a very right-skewed distribution of mussel recruitment among years. For example, Beukema’s 1970–2008 series at Balgzand, a tidal flat in the westernmost Dutch Wadden Sea, only contained four years with good mussel recruitment (Beukema and others 2010). A recent study by Beukema confirms these results (Beukema and others 2015). This is a general phenomenon for most bivalve species in soft-sediment habitats. The same Beukema series revealed six outstanding recruitment years for the Baltic tellin Limecola balthica and also only four years for the cockle Cerastoderma edule. In the Wash, spatfall of cockles and mussels was poor in most years between 1990 and 1999. Significant recruitment that increased the shellfish stock occurred in only 14% and 19% of the study years for mussels and cockles, respectively (Dare and others 2004). Thus, as already indicated by the simulation study, full recovery of bivalve bed area after the disappearance in the early 1990 s to pre-1990 levels might take some time. If one also considers the uncertainty about the pre-1990 area and the expectation that future beds will most likely be mixed ones, it seems premature to invest effort in costly large-scale restoration programs.

Using beds as the unit of observation, and estimating their recruitment and survival, provided some insight in the underlying processes and relevant environmental factors that govern total bed area. Though it is of interest to managers and conservationists, it does not tell much about the fate, in terms of recruitment and survival, of individual mussels and oysters. Beds are continuously replenished with new cohorts (Figure 9), and we have not investigated the link between bed survival and recruitment and survival of individual animals. So we do not know, for example, whether most recruitment in terms of individuals occurs in existing beds or in newly formed beds. It also cannot be ruled out that stable beds with low hazard rates are linked with relatively low survival of individuals, but high recruitment of new cohorts. Perhaps beds that are constantly renewed are more stable than beds that lack regular replenishment. The link between these two levels of biological organization, the individual bivalve and the bed is an interesting topic for future studies.

Figure 9
figure 9

Size–frequency plots of mussels on a bed at Balgzand, showing that beds are constantly renewed by the disappearance of old cohorts and the renewal with recruits (unpublished data from Andreas Waser. See Waser and others (2016) for the precise location of the bed).

Full size image