Introduction

Over the last decades it has become well established that some organisms can have disproportionally strong effects on their abiotic environment, indirectly affecting other species. Such species, often called ‘ecosystem engineers’ by Jones and others (1994), typically promote their own preferred conditions at the local (‘patch’) scale (Bertness and Leonard 1997; Rietkerk and others 2004 and references therein). However, ecosystem engineering is often not only important locally, but may also have strong impacts at landscape scales (Wright and others 2002; Kefi and others 2007; Scanlon and others 2007). Apart from altering the spatial structure of the environment, ecosystem engineers may affect the spatial distribution and abundance of their resources (for example, nutrients, water, and light). This also alters resource availability for other species (Gutierrez and others 2003; van de Koppel and others 2006), which should in turn affect the spatial distribution of their consumers (for example, Hassell and May 1974; Folmer and others 2010; Piersma 2012). Although effects of prey-patchiness and ecosystem engineering on the distribution of species have been documented separately (for example, Hassell and May 1974; Wright and others 2002), assessments of the spatially extended effects of ecosystem engineers on resources, and their consumers have remained largely theoretical (Bagdassarian and others 2007; Olff and others 2009).

Reef builders like blue mussels (Mytilus edulis) and Pacific oysters (Crassostrea gigas) are striking examples of ecosystem engineers that impact their environment through habitat modification (Kröncke 1996; Gutierrez and others 2003; Kochmann and others 2008). At a local scale, mussels and oysters create hard substrate and increase habitat complexity, reduce hydrodynamics, and modify the sediment by depositing large amounts of pseudo-feces and other fine particles (Kröncke 1996; Hild and Günther 1999; Gutierrez and others 2003). However, in soft-bottom systems, their effects on sediment conditions typically extend beyond the direct surroundings of the reefs and may be detectable up to several hundreds of meters (Kröncke 1996; Bergfeld 1999). Many studies have demonstrated that reef builders have an important effect on the local benthic community (Dittmann 1990; Norling and Kautsky 2008; Markert and others 2009) and that the reefs themselves are important foraging grounds for avian consumers (for example, Nehls and others 1997; Caldow and others 2003). However, the spatially extended effects of such reef builders on this community remain largely unstudied.

Furthermore, possible implications of such spatially extended habitat modification on the community may also be important from a management perspective. In many intertidal soft-sediment systems, like the Wadden Sea, ecosystem engineers have disappeared due to multiple anthropogenic disturbances and many associated species disappeared with them (Piersma and others 2001; Lotze and others 2005; Kraan and others 2007; Eriksson and others 2010). For instance, in the Wadden Sea, 150 km2 of seagrasses disappeared in the 1930s (van der Heide and others 2007) and mussel beds were almost completely removed in the beginning of the 1990s and have only partly recovered thus far (Beukema and Cadee 1996). If spatial effects of ecosystem engineers are not recognized, such dramatic changes might result in unexpectedly strong losses in these ecosystems.

In this article, we investigate the effects of spatial habitat modification by mixed blue mussel and Pacific oyster reefs on the distribution of benthic prey and their consumers (shorebirds) at a sandy intertidal flat. We collected spatially explicit data on important abiotic variables and the biota in and around two reefs in the Dutch Wadden Sea. We used Structural Equation Modeling (SEM) to infer the relative importance of ecosystem engineering on the spatial distribution of recourses and consumers. Based on sediment and benthos data of 119 sampling stations at varying distances from the reefs and the spatial mapping of shorebirds, we constructed default models for four of the most commonly observed bird species that included all possible interactions between the birds and their environment. Next, we determined the relative importance of each interaction, using an approach with stepwise exclusion of variables.

Materials and Methods

Study Area

Our study area covered about 44 ha of intertidal mudflats, south of the island of Schiermonnikoog in the eastern Dutch Wadden Sea (53°28′15.75″N, 6°13′20.06″E). These intertidal flats contain a variety of macrobenthic invertebrate species (Beukema 1976) that are accessible to shorebirds twice a day (van de Kam and others 2004; van Gils and others 2006). The area contained two mixed reefs of blue mussels and Pacific oysters, established in 2002 (Goudswaard and others 2007 and unpublished data of our research group). The main cohort of bivalves was 7 years old, with several younger cohorts. Before the establishment of the two reefs, our study area consisted of a sandy intertidal flat without patches of hard substrata (van de Pol 2006 and unpublished data of our research group). The spatial relationships of the reef builders with the local and surrounding benthic community and associated shorebirds were examined at two adjacent study areas of 22 ha each (see Figure 2).

Benthic Sampling

Sediment, pore water, and benthic samples were collected in August 2009 on a predetermined 100 m grid with 46 additional random points. In total, 119 station points were sampled across the two study sites. All stations were identified during low tide using a handheld GPS. At each sampling station, we sampled and pooled three 5 cm deep sediment cores with a PVC corer with an area of 7.1 cm2. Sediment organic matter content in dried sediment (24 h, 70°C) was estimated as weight Loss On Ignition (LOI; 5 h, 550°C). Silt content (% sediment fraction < 63 μm) was determined by a particle size analyzer (Malvern). Redox potential was measured immediately after sampling with a multi-probe meter (556 MPS, YSI) in pore water that was extracted from the sediment with a ceramic cup into a vacuumized 50 ml syringe. Benthic samples were taken with a stainless steel core with area of 179 cm2 down to a depth of 20–25 cm. Samples were sieved over a 1 mm mesh and all fauna fixed in 4% formalin. In the laboratory, samples were stained with Rose Bengal, and fauna was identified to species level. Ash free dry mass (AFDM) of each species was determined by LOI (5 h, 550°C) after drying for 48 h in a stove at 60°C.

Bird Mapping

A 3.2 m high observation platform was constructed 100 m away from each of the two study sites in such a way that the platforms covered the respective sampling grids, that is, a reef and the associated gradient towards a sandy area, all within a radius of 500 m. The spatial distribution of shorebirds was determined during four tidal cycles between 18 August and 8 September 2009. Positions of individual birds were determined using the newly developed Telescope-Mounted Angulator (TMA) described by van der Heide and others (2011). This was done from an hour before to an hour after the time of low water, that is, when the areas were completely exposed and tidal movement would not affect their spatial distribution. With the TMA, using trigonometry, we were able to determine a bird’s spatial position with high accuracy (maximum prediction error of 8.7 m at 500 m; van der Heide and others 2011).

We mapped the spatial distribution of four common shorebird species: oystercatcher (Haematopus ostralegus), Eurasian curlew (Numenius arquata), bar-tailed godwit (Limosa lapponica), and black-headed gull (Chroicocephalus ridibundus). These focal species were chosen for three reasons. First, due to their body size, all four species are easy to follow and clearly visible, which prevented double counting and inaccurate positioning (van der Heide and others 2011). Second, all four species form sparse flocks, a feature that represents a degree of sensitivity to interference of conspecifics (Goss-Custard 1980; Piersma 1985). In contrast to social and interference-insensitive species, the distribution of such interference-sensitive species should mostly be determined by the distribution of food resources (Folmer and others 2010). Third, each of these species should differ in its degree of association with mussel and oyster reefs. For example, as blue mussels form a substantial part of their diet, oystercatchers tend to be highly associated with reef builders (Goss-Custard 1996). Eurasian curlew typically respond to an increased abundance of crabs and shrimps in and near reefs compared to sandy intertidal flats, but they also feed on bare mudflats (Goss-Custard and Jones 1976; Petersen and Exo 1999). The degree of association for bar-tailed godwits is probably lower, because they feed on a large variety of benthic animals often along the edge of the receding and advancing tide (Goss-Custard and others 1977; Scheiffarth 2001). Black-headed gulls feed on a large variety of prey and can be found in many different habitats (Dernedde 1994; Kubetzki and Garthe 2003).

Data Analysis

Both study sites were subdivided by Thiessen polygons (Thiessen 1911) in ArcGIS (Environmental Systems Research Institute, Redlands, California, USA). Each polygon defines a discrete area around each sampling station (both random and predetermined) in such a way that any location inside the polygon is closer to that point than to any of the neighboring points. No great differences were detected between shorebird numbers during the four tidal cycles, so data were pooled to calculated densities. Densities of each bird species (# ind. m−2) were calculated for each polygon and merged into a single master dataset that now contained data on abiotic variables (sediment organic matter, silt, and redox), biomass of all benthic species, and bird densities for each sampling station. To approach a normal distribution for analyzed variables, organic matter content was reciprocally transformed (y = 1/x), redox potential was log transformed (y = log10(x)) and all other variables were square root transformed (y = √x).

Next, we used SEM (Amos v18) to test the spatial effects of the reefs on abiotics and the possible direct and indirect effects on the distribution of macrobenthic and bird species. For each bird species, we created default models that included all potentially important causal relationships between straight-line distance to the center of the reef (calculated in ArcGIS), directional effects that may arise from strong winds or currents (calculated in ArcGIS as the deviation of each station from the north–south axis through the center of the reef), sediment conditions (organic matter, silt fraction, and redox), macrobenthos biomass, and bird density (Figure 1). These models focus on explaining shorebird distribution from information on underlying resources. Therefore, each model only included macrobenthos species that are known prey items for that particular bird species (Table 1). Apart from modeling the effect of prey density on shorebird distribution, the models also tested for possible relationships between sediment variables, distance to the reef, and bird density. Sediment conditions can, directly or indirectly, affect bird distribution (Myers and others 1980; Yates and others 1993; Johnstone and Norris 2000). Furthermore, distance to the reef might influence bird distribution because birds may be attracted to these areas in anticipation of altered sediment conditions and prey densities. In summary, all four default models include (Figure 1): (1) effect of distance and direction to the reef on sediment variables, (2) the effect of sediment variables on macrobenthos, (3) effects of macrobenthos variables on bird density, (4) direct effects of sediment variables on bird density, and (5) effect of distance to the reef on bird density.

Figure 1
figure 1

The conceptual path analysis model. Arrows depict direct effects of one variable (boxes) on another. Numbers represent specific mechanisms described in “Materials and Methods” and “Data Analysis”.

Table 1 Variables Included in the Model to Test the Default Model for Each Focal Bird Species

To test whether the identified relationships extended beyond the reefs themselves, we analyzed each model twice—once with all data points included (119 stations) and a second time with the stations inside the reefs excluded (111 stations). Models were analyzed with stepwise backward elimination of relations included in the default model (threshold significance for elimination: p < 0.05). After each elimination step, we used the χ 2 test (probability level > 0.05) to test for an adequate fit (that is, that observed data did not differ significantly from those predicted by the model), and compared the model to previous models using Akaike’s Information Criterion (AIC). Unidentified models were excluded from the results. We also excluded macrobenthic species from the model if they were not correlated with the modeled bird species, whereas sediment conditions were omitted if they were not related with either macrobenthic species or bird density. Furthermore, when abiotic or benthic variables exhibited strong significant collinearity (r > 0.4) without one explaining the other (for example, different proxy’s for sediment conditions), we only included the variable with the highest explained variation in our models. The latter was done because SEM models become notoriously unreliable when relations with very strong covariance are included (Petraitis and others 1996; Grewal and others 2004).

Results

Organic matter, silt content, and redox were all highly correlated (r-values for OM-silt, OM-redox, and silt-redox were 0.9, 0.5, and 0.5, respectively) and exhibited strong spatial gradients, with organic matter and silt increasing and redox decreasing in the direction of the reef. A map overlay of organic matter and the distribution of the four shorebird species suggest that oystercatchers, and to a lesser extent also curlews and bar-tailed godwits, tend to aggregate in these organic matter-rich areas in and around the reefs (Figure 2). In contrast, the spatial distribution by black-headed gulls appears much less affected by the presence of the reefs.

Figure 2
figure 2

Overview of the two reefs and their surrounding intertidal flats, showing Thiessen polygons (each polygon contains one sampling station), the position of the reefs (striped black areas), and the distribution of sediment organic matter content in relation to the distribution of A oystercatchers, B curlews, C bar-tailed godwits, and D black-headed gulls. Black dots represent the positions of the birds. Circles with a black dot indicate the position of the observation platforms.

Organic matter was included as a proxy for sediment conditions in the SEM models instead of silt content or redox because of its highest explained variation (R 2’s were 0.45, 0.31, and 0.43, respectively). The distributions of several macrobenthic species were strongly affected by sediment organic matter, which in turn explained a significant part of the distribution of all four shorebirds (Figure 3; Appendix 1A, B in Supplementary material). The correlations suggest that organic matter had a positive effect on the biomass of Lanice conchilega, Hediste diversicolor, Cerastoderma edule, and crustaceans (explaining 7, 12, 39, and 11% of their variance, respectively).

Figure 3
figure 3

Diagram of the SEM results with all sampling stations for oystercatchers, curlews, bar-tailed godwits, and black-headed gulls (AD) and the sampling stations inside the reefs excluded for the same bird species (EH). Arrows indicate significant direct effects. The thickness line of each arrow indicates the magnitude of the standardized path coefficient, which is presented numerically next to each path. The R 2 values adjacent to the boxes represent the total variance explained by all significant predictors (*p < 0.05; **p < 0.01; ***p < 0.001).

All default models based on the fully saturated model (Figure 1) and species-specific feeding relations (Table 1) demonstrated poor model-data fits (Table 2). After stepwise backward elimination and removal of non-significant relations, all final models demonstrated a strong fit. In contrast to the default models, final models demonstrated low Chi-square values, a probability level above 0.05 and low AIC’s (Table 2; Appendix 1 in Supplementary material). After removing the sampling stations within the reefs from the dataset, all final models still had an adequate fit and the structure of the models remained nearly identical (Table 2; Figure 3). The models including the sampling stations on the reefs yielded a slightly better fit for oystercatchers, bar-tailed godwits, and black-headed gulls, whereas the model for curlews improved after removing the reef stations.

Table 2 Model Fit Summary from SEM for the Default Model and the Final Modified Model for the Dataset with All Sampling Stations Included and for the Dataset Wherein the Sampling Stations inside the Reefs were Excluded

The final models for each bird species revealed significant correlations with macrobenthic species, but also with abiotic variables. Distance to the reef, organic matter and C. edule, were significant predictors of oystercatcher density (Figure 3A, E), with the final model explaining 62 (including local effects) to 59 (excluding local effects) % of the variance. The standardized effect of distance to the reef on oystercatcher density (−0.417 to −0.380) was stronger than the effect of organic matter (0.338 to 0.331) and biomass of C. edule (0.152 to 0.179). For curlews (51 to 44% of the variance explained), crustaceans and distance to the reef were significant predictors for both models (Figure 3B), whereas H. diversicolor was dropped in the model that excluded the reef effect (Figure 3F). The standardized effect of distance to the reef on curlew density (−0.597 to −0.617) was larger than the effect of crustacean biomass (0.195 to 0.168) and biomass of H. diversicolor (0.141, only in the model which included local effects). L. conchilega and organic matter were the two significant predictors of densities of bar-tailed godwits (Figure 3C, G). The standardized effect of organic matter on bird density (0.370 to 0.386) was larger than that of the biomass of L. conchilega (0.227 to 0.187) and the final models both explained 23% of the observed variance. Finally, for black-headed gulls, organic matter and C. edule were significant predictors of density (Figure 3D, H). The standardized effect of C. edule on black-headed gull density (0.372 to 0.405) was larger than the effect of organic matter (−0.303 to −0.289). The final models explained 9 to 8% of the variance.

Discussion

Although ecosystem engineering can determine the spatial distribution of resources (for example, Gutierrez and others 2003; van de Koppel and others 2006) and resources in turn importantly control the distribution of consumers (for example, Nachman 2006; Folmer and others 2010; Piersma 2012), the interaction between these two processes so far has rarely been examined (Olff and others 2009). Here, we demonstrate that ecosystem engineers can affect consumer-resource interactions far beyond their own physical spatial boundaries in intertidal soft-sediment systems. Reef-building bivalves like mussels and oysters cover a relatively small part of the intertidal mudflats of the Wadden Sea (±1%). Our results, however, imply that their ecological impact is much larger than their size may suggest.

We found strong spatial gradients of increasing sediment organic matter and silt fraction and decreasing redox potential in the direction of mixed mussel and oyster reefs, which in turn affected the distribution of benthic species. Moreover, distance from the reefs, sediment characteristics, and prey abundance simultaneously affected the distribution of the three studied species that have more or less specific prey requirements (oystercatchers, curlews, and bar-tailed godwits). This is most likely because the birds feed in the modified areas in anticipation of higher prey abundances. Black-headed gulls, the only species that did not cluster on and around the reefs, are versatile foragers with many diet options and this may explain why the reefs and the modified areas did not affect their spatial distribution. When the data points for the reefs themselves were excluded from the statistical analysis, the outcomes did not change, thus emphasizing the importance of the spatially extended effects of reefs. Only the ragworm H. diversicolor was excluded from the model as a predictor for the distribution of curlews. This was, however, understandable as ragworms were mostly found in muddy sediments in and around the mixed reefs.

Community structure alteration by ecosystem engineers through spatially extended habitat modification seems to occur in many different ecosystems including beaver-inhabited wetlands (Wright and others 2002) and cordgrass-inhabited cobble beaches (Bruno 2000). However, the relevance of habitat modification by ecosystem engineers on its surrounding and higher trophic levels may vary with environmental conditions. For instance, although our results show that habitat modification by reef builders can be pronounced and exceed the spatial boundaries of the reefs themselves, spatial engineering effects by the same species on rocky shores are typically more limited. In these systems, blue mussels modify environmental conditions mainly by providing structural protection for associated fauna (Thiel and Ullrich 2002; Gutierrez and others 2003). Hard substrate is already present and fine particles produced by mussels (feces and pseudofeces) are washed away by more intense hydrodynamics, resulting in more limited modifications at larger spatial scales (Thiel and Ullrich 2002). Furthermore, the effect of habitat modification by reef builders may also interact with the presence of other ecosystem engineers. For example, the tube-worm Lanince conchilega is also considered as an ecosystem engineer in soft-sediment systems, as their tubes provide substrate and facilitate the deposition of fine sediments (Friedrichs and others 2000; Zühlke 2001). Because the presence of L. conchilega is positively correlated with the abundance and richness of the benthic community (Zühlke 2001; Callaway 2006; Godet and others 2011), L. conchilega may locally enhance the engineering effect of the reefs on the benthic and shorebird community.

In our study, SEM proved to be a useful tool for disentangling the relative importance of consumer-resource interactions and spatial habitat modification by ecosystem engineers. Using stepwise backward elimination of significant relations, we obtained models with reliable fits of multiple ecologically relevant variables. The method is correlative and therefore does not provide any direct evidence. Ideally, this method should be complemented with other, more direct approaches like smaller-scale manipulative experiments. However, before the reefs established themselves 7 years ago the study area was sandy and homogeneous, and in this respect the study reported here can be regarded as experimental (but in want of detailed description of the re-establishment situation).

In conclusion, our results indicate that consumer-resource interactions can be affected by reef builders far beyond the spatial boundaries of the reefs. This implies that these reefs have a much larger ecological impact on the intertidal community than their actual size suggests, which in turn means that loss of ecosystem engineers may result in disproportionally large consequences for biodiversity values in protected intertidal areas, like the Wadden Sea. Although the Pacific oyster is an alien species that invaded the Wadden Sea in the late 1970s (Troost 2010 and references therein), recent studies showed that oyster reefs might compensate for the large loss of mussels in 1990–1991 by replacing the ecological function of blue mussel reefs (Kochmann and others 2008; Markert and others 2009; Troost 2010). Nevertheless, the effects of Pacific oysters on the intertidal community and trophic interactions should be further investigated. Overall, our study emphasizes that conservation and restoration of reef builders should be considered a crucial step in the restoration of such systems.