Spatial variability in macrofaunal diet composition and grazing pressure on microphytobenthos in intertidal areas

Microphytobenthos forms an important part of the diet of macrofauna (macrozoobenthos) in many intertidal ecosystems. It is unclear, however, whether the dependence of macrofauna on microphytobenthos varies spatially within and among tidal systems. We aim (1) to assess the spatial variability in the importance of microphytobenthos in the diet of macrofauna (i.e., between and within two tidal basins and as function of elevation), (2) to quantify grazing pressure of the macrofaunal community on different potential food sources (microphytobenthos, phytoplankton and terrestrial organic material) for several sites in two tidal basins and (3) to compare microphytobenthic production and summer/autumn grazing of the total macrofaunal community and grazing pressure per feeding type, with potential microphytobenthic production estimated from rates in early spring, when grazing was low. Using a natural stable isotope approach, we identified microphytobenthos as a more important food source for macrofauna than phytoplankton and terrestrial organic material. Microphytobenthos dependency differed between tidal basins for the genera Bathyporeia (sand digger shrimp), Macoma (Baltic tellin), and Peringia (mudsnail) and for sampled individuals of all genera combined, and did not vary as function of elevation. We showed that macrofaunal grazing on microphytobenthos is quantitatively important and, in some cases, approached microphytobenthic production rates in early spring. No positive relation between microphytobenthic production in early spring and macrofaunal grazing in summer/autumn was observed. This suggests that the studied consumer‐resource interactions are coupled on a larger spatial scale (i.e., mesoscale, ≈10 to 100 km), rather than the fine (mm to m) scale.

in the diet of macrofaunal species. The relative importance of phytoplankton in the diet of the semelid bivalve Theora lubrica, which inhabits estuarine subtidal sediments, increases seaward in Gokasho Bay, Japan (Yokoyama and Ishihi 2003). Furthermore, the proportion of macroalgae in the diet of herbivores was found to increase with the availability of these macroalgae among three sub-estuaries (Olsen et al. 2011). Christianen et al. (2017) demonstrated large spatial heterogeneity on the scale of tens of kilometers in the contribution of microphytobenthos in the diet of macrofaunal species in the Wadden Sea (the Netherlands), which could, however, not be attributed to a specific factor. In estuaries, major abiotic gradients are presently associated with distance from the estuarine mouth (salinity, temperature, hydrodynamics, and sediment composition) and gradients from tidal flats to subtidal channels (elevation and hydrodynamics, and resulting depth of overlying water and sediment composition of the bed) ( Van der Wal et al. 2017). Some papers have used the spatial distribution and temporal dynamics of microphytobenthic biomass to correlatively explain patterns in macrofaunal species composition (e.g., Van Colen et al. 2008) and macrofaunal biomass (Van der Wal et al. 2008). However, it is not yet clear if and how the dependence of macrofaunal species on microphytobenthos differs between tidal basins, or varies with the aforementioned abiotic gradients.
Field experiments where macrofauna were removed showed that macrofauna can exert significant top-down control on microphytobenthos and thereby constrain microphytobenthic biomass (e.g., Weerman et al. 2011). Pratt et al. (2015) found a negative relationship between recent deposit feeding activity of the bivalve Macomona liliana and microphytobenthic biomass at a scale of 10s of meters in a subtropical Manukau Harbor, New Zealand. Grazing pressure can be quantified using 13 C-labelling . Middelburg et al. (2000) showed at a relatively sandy tidal flat in the Westerschelde estuary (southwest Netherlands) that during a period of 2.4 d half the amount of labeled microphytobenthos was removed, of which 40% was respired and 60% was resuspended. However, Middelburg et al. (2000) emphasized that the resuspended fraction may depend on weather conditions and bioturbation by macrobenthos, while the respired fraction may depend on the benthic community composition. Moerdijk-Poortvliet et al. (2018) found large (temporal) variation in the net loss of carbon fixed recently by microphytobenthos in the Oosterschelde (southwest Netherlands).
Microphytobenthic biomass shows a strong temporal variability in temperate regions, which is often seasonal (Sundbäck et al. 2000), although exceptions are found (Thornton et al. 2002).
In the Westerschelde, a microphytobenthos bloom generally occurs in March/April and disappears in May/June likely due to top-down control (Weerman et al. 2011). Recruitment of larvae (using microphytobenthos) starts approximately in May (Van Colen et al. 2008). After larval recruitment, macrofauna is relatively sessile (Herman et al. 1999). Estuary scale averages of microphytobenthic production measured in June have been linked to averages of the proportion of macrofaunal biomass that depends on microphytobenthos (macrofauna sampled in June 1 yr after microphytobenthic production was measured) . It is unclear, however, to what extent microphytobenthic production and macrofaunal grazing pressure are linked on smaller spatial scales (within tidal basins).
In this research, we aim (1) to assess the spatial variability in the importance of microphytobenthos in the diet of macrofauna using a natural stable isotope approach, (2) to quantify grazing pressure of the macrofaunal community on microphytobenthos, phytoplankton, and terrestrial organic material for several sites in two tidal basins, and (3) to quantify the (indirect) relationship between microphytobenthic production in early spring (when grazing pressure is still low) and grazing pressure of the total macrofaunal community and of separate macrofauna feeding types in summer/autumn. The study is conducted in two contrasting tidal basins in terms of species assemblages and salinity (Ysebaert et al. 2003), at contrasting elevations in the tidal zone (high and low). For the first aim, we hypothesize that the relative importance of microphytobenthos in the diet of macrofauna only varies for suspension feeders and facultative suspension/deposit feeders, depending on spatial variability in the availability of microphytobenthos and phytoplankton. Furthermore, it is hypothesized that the relative importance of microphytobenthos in the diet of macrofauna increases with elevation, as the availability of microphytobenthos biomass is typically positively linked to elevation (Van der Wal et al. 2010). For the third aim, we hypothesize that grazing pressure of the total macrofaunal community sampled in summer/autumn depends on initial microphytobenthos production in early spring. Microphytobenthos production in early spring may be a proxy for potential food availability for macrofauna in summer/autumn, as during early spring the influence of grazing by macrofauna on spatial variability in microphytobenthic production is limited (Weerman et al. 2011).

Study sites
The study was performed in the Westerschelde estuary and Oosterschelde tidal basin in the Netherlands, two mesotrophic, tide-dominated systems (Fig. 1). The Westerschelde estuary has a salinity gradient varying from mesohaline at the Dutch-Belgium border to polyhaline at the mouth. It has an intertidal surface area of approximately 70 km 2 (Van der Wal et al. 2010) and a total surface area of about 312 km 2 (22% intertidal area). In the Westerschelde, phytoplankton productivity varies from approximately 100-300 g C m −2 y −1 along the estuarine gradient (10-30‰) . Deposit feeders form an important part of the total biomass of macrofauna in the (eastern part of) the Westerschelde (Ysebaert and Herman 2002) and the biomass of suspension feeders is relatively low (41%, Herman et al. 1999).

Daggers et al.
Diet and grazing pressure of macrofauna The Oosterschelde is a polyhaline, semi-enclosed sea-arm, only closed at the mouth during storms, and with only very limited river input (Nienhuis and Smaal 1994). The Oosterschelde has an intertidal surface area of approximately 50 km 2 (Van der Wal et al. 2010) and a total surface area of around 350 km 2 (14% intertidal area). It has a phytoplankton productivity of approximately 155 g C m −2 year −1 (Smaal et al. 2013). Suspension feeders compose a relatively large part (82%) of the total macrofaunal biomass; the remaining biomass mainly consists of deposit feeders (Herman et al. 1999).

Field methods
Sampling stations were selected from a set of existing stations, from which historical macrofauna data were available. The density and biomass of the historical macrofauna data was used to calculate macrofaunal grazing pressure, along with macrofauna samples collected in this study to determine the macrofaunal diet composition using natural stable isotopes. A stratified random selection was applied within the zone − 1 to 0 m NAP ("low" elevation; where NAP is Dutch Ordnance level, approximately mean sea level) and within the zone 0 to 1 m NAP ("high" elevation). Six study sites were selected in the Westerschelde and three sites in the Oosterschelde (Fig. 1). A historical macrofauna community composition dataset was collected by the Royal Netherlands Institute for Sea Research (NIOZ) for Rijkswaterstaat yearly in summer/autumn (August-October, 2003-2012 in the BIOMON programme using a 15 cm diameter cylinder corer, up to a depth of 30 cm, and sieved on a 1 mm mesh size sieve. In this data set, in the Oosterschelde, samples were collected at the same stations each year while in the Westerschelde a fixed number of samples was collected randomly from each ecotope each year (Bouma et al. 2006). For our field campaigns (see below), 26% of our selected stations were not sampled at the exact same location as NIOZ/Rijkswaterstaat sampling stations due to limited data availability, but within a distance of 300 m. Macrofaunal species distributions are strongly related to sediment grain size distributions (Van der Wal et al. 2008) and are relatively homogeneous in fine sands at scales from 50 cm to 500 m, but more heterogeneous at distances of > 50 m in muddy sediments (Kendall and Widdicombe 1999;Ysebaert and Herman 2002 Within a plot, three sediment cores (ø 15 cm) were taken during both sampling periods up to a depth of 1 cm using a 3 cm diameter syringe for the analysis of grain size distribution and chlorophyll-a (Chl a) content, which was used to calculate microphytobenthic production. Sediment samples up to a depth of approximately 1 cm were collected with a spoon for phospholipid derived fatty acids PLFA extractions for analysis of the δ 13 C of PLFAs characteristic for diatoms. In addition, sediment samples up to a depth of approximately 1 cm were collected with a spoon for stable isotope analysis (δ 13 C and δ 15 N) of bulk sediment. Water (3-L jars) was collected at each site during low tide at the edge of the tidal flat for the analysis of the stable isotope composition of Suspended Particulate Organic Matter (SPOM) as indicator for phytoplankton. It has been demonstrated that the δ 13 C of POC in the water column (AE −26‰) deviates little from the isotopic signature of phytoplankton (− 30‰) in the middle part of the Westerschelde considered in this study (Hansweert and Zandvliet; Van den Meersche et al. 2009). Also Currin et al. (1995) report similar values of the δ 13 C of suspended particulate organic matter and phytoplankton. In the Oosterschelde, no riverine input of organic material is present and, consequently, terrestrial organic material is scarce. Therefore, the isotopic signal of SPOM was assumed to be representative for the isotopic signal of phytoplankton in the Oosterschelde. Suspended particulate organic matter (SPOM) is a frequently used indicator for phytoplankton (Galván et al. 2008). Photosynthetic activity was measured in each plot in triplicate by constructing a rapid light curve (RLC) using a Pulse Amplitude Modulation (PAM) fluorometer (Mini PAM; Heinz Walz GmbH). In late spring, species of macrofauna that were most commonly observed in the NIOZ/Rijkswaterstaat dataset were collected using three sediment cores (ø 15 cm) per plot for analysis of the diet composition using natural stable isotopes. Cores were sieved on a 1 mm mesh size sieve and stored in open jars containing water collected on site, to keep benthic species alive (to preserve the isotope signal without transporting ice into the field). The ragworm Hediste diversicolor was kept in separate jars to prevent the species from predating on other specimens. After returning to the lab, the water was removed and macrofauna samples were stored frozen (− 20 C).

Laboratory methods
Sediment samples were transported in ice and stored in − 80 C within 8 h after sampling. Sediment samples were freeze dried for minimally 48 h and the dry bulk density of sediment samples was determined gravimetrically. Of each sediment sample, 700 mg was used for pigment analysis. Ten milliliters of acetone (90%) and glass beads were added to extract photosynthetic pigments using a cell homogenizer (Braun, Type 8530220) for 20 s. Subsequently, the extract was centrifuged (2000 rpm, 3 min) at 20 C immediately after homogenization and light absorption was measured with a Specord 210 spectrophotometer. Chl a content (μg g −1 ) was calculated according to the equation described in Ritchie (2006) and also includes phaeopigments. These were converted to Chl a concentrations (mg m −2 ) using the dry bulk density. The sediment grain size distribution was determined using a Malvern (laser) particle sizer.
PLFAs were extracted from approximately 6 g of wet sediment using a modified Bligh and Dyer extraction (Boschker et al. 1999). The followed procedure is described in detail in Middelburg et al. (2000). Fatty acid methyl esters (FAME) were identified based on the comparison of the retention time of sampled FAME with internal FAME standards (12:0, 16:0 and 19:0). The isotopic composition of FAME was measured with a Varian 3400 gas chromatograph containing a Varian SPI injector, which was coupled to a Finnigan Delta S isotope ratio mass spectrometer via a type II combustion interface. Carbon isotope ratios of the PLFAs were corrected for the addition of one carbon atom in the methyl group during derivatization. The carbon isotope ratios of the 16:2ω4, 20:5ω3, and 22:6ω3 fatty acids were used as proxy for the carbon isotope composition of benthic diatoms. The fatty acids are characteristic for diatoms, dinophytes, and haptophytes (Dijkman and Kromkamp 2006) and are a suitable indicator for microphytobenthos at our study site, as microphytobenthos assemblages are generally dominated by diatoms in the Scheldt estuary (Sabbe and Vyverman 1991). Fractionation within the diatom cell was accounted for by adding 5.4 per mil to the δ 13 C of 16:2ω4, according to Schouten et al. (1998). Fractionation factors within the diatom cell of 20:5ω3 and 22:6ω3 were calculated using the ratio between δ 13 C values of 16:2ω4 and 20:5ω3 (1.07 AE 2.21), and 16:2ω4 and 22:6ω3 (2.1 AE 2.81) in the dataset used in this study. The δ 13 C of diatom cells was calculated as the weighted average of measured δ 13 C values of the three fatty acids using the relative concentrations as weights.
Water samples for determination of the isotopic composition of SPOM were filtered onto a glass fiber filter (Ø47 mm, Whatman ref no 421026). The δ 13 C and δ 15 N of bulk suspended particulate matter and bulk sediment material was measured using a Fisons elemental analyzer coupled to a Finnigan delta S isotope ratio mass spectrometer after samples were acidified using an in situ acidification method to remove inorganic carbonates (Nieuwenhuize et al. 1994).
To test whether spatial variability exists in the diet of macrofauna, 10 common species were selected as target species (see Supporting Information, Table S2) and retrieved from macrofauna samples collected at "high" and "low" sampling stations in late spring. Hereby, we ranked the species according to their occurrence in the highest number of locations, based on abundance data in the historical macrofauna dataset (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) collected by the Royal Netherlands Institute for Sea Research (NIOZ) and Rijkswaterstaat from the Westerschelde and Oosterschelde combined. It was ensured that all feeding types were represented in the species selection. Feeding types were retrieved from the historical macrofauna Daggers et al. Diet and grazing pressure of macrofauna dataset and updated using recent literature. Macoma was considered a facultative suspension/deposit feeder (Rossi et al. 2004). From each plot, three individuals of each target species present in the sample were retrieved, identified to genus level, pooled, and freeze dried for minimally 48 h. Macoma balthica was divided into three size classes: 5-10, 10-15, and > 15 mm shell length, and pooled per size class for stable isotope analysis. Macoma balthica may show a gradual shift in diet with small juveniles feeding entirely on microphytobenthos and larger individuals depending more on phytoplankton (Rossi et al. 2004). Soft parts of the bivalves were removed from the shells and the mudsnail Peringia ulvae was treated with 2 N HCL to remove inorganic carbonates. Gut contents were not removed. Herman et al. (2000) demonstrated that gut contents had a minimal influence on the bulk δ 13 C of macrofauna; even after tracer addition with high δ 13 C of the food, the influence of the gut content on the δ 13 C of the bulk organism was 1%-5% after 4 d in their study. Therefore, the influence of natural gut contents can be considered negligible (Rossi et al. 2004).

Calculation of microphytobenthos primary production
The rapid light curves measured with the PAM in early spring 2015 were fit using the model of Eilers and Peeters (1988) rewritten by Herlory et al. (2007) for all measurements above the detection limit: Ft > 200. This yielded the photosynthetic activity (α) and photosynthetic capacity (P s ). Microphytobenthic production rates were then calculated based on these two parameters, as well as incident irradiance, Chl a concentrations and mud content (% particles < 63 μm) as described in  and Daggers et al. (2019). Monthly averages of daily microphytobenthic production rates were calculated using a sediment-optical model, taking into account temporal variability in irradiance and emersion duration (as a function of tides and bathymetry); production during immersion was modeled using light attenuation values within the water column measured at the nearest monitoring stations averaged over the period 2011-2016 Daggers et al. 2019). Daily production rates were calculated over 1 month corresponding to the sampling period.

Estimation of macrofaunal diet composition
The fraction of microphytobenthos in the diet of macrofauna sampled in late spring was estimated assuming that the isotopic signature of macrofauna is a weighted average of its food sources, plus a fractionation factor. A Bayesian mixing model (MixSIAR; Stock and Semmens 2013) was used to estimate the proportional contribution of food sources to the diet of macrofauna from the δ 13 C and δ 15 N in the Oosterschelde and Westerschelde. The diet composition was estimated for each selected (most common) species per site per elevation category. δ 13 C values of microphytobenthos were determined from PLFAs extracted from sediment sampled in each plot. To obtain δ 15 N values of benthic algae, δ 15 N values of bulk sediment were used + 2‰ for the Oosterschelde and + 10‰ for the Westerschelde, using the difference between SOM (Sedimented Organic Matter) and benthic algae in the Oosterschelde and Westerschelde reported by Riera et al. (2000). The average values of the δ 15 N of the sampled bulk sediment in the Oosterschelde (6.7 AE 0.9‰, n = 18) and Westerschelde (9.1 AE 1.3‰, n = 36), respectively, closely resembled the average δ 15 N of the SOM in the Oosterschelde (7.3 AE 0.1‰, n = 3) and Westerschelde (8.6 AE 0.2‰, n = 8) reported in Riera et al. (2000). As proxy for phytoplankton, δ 13 C and δ 15 N values of bulk suspended particulate matter from the site of interest and the two nearest sites were used. δ 13 C and δ 15 N values of the macroalgal genus Ulva which dominates macroalgal spring blooms in the Oosterschelde (Rossi 2006) derived from Riera et al. (2002) were added as potential food source in the Oosterschelde. Diet compositions were calculated with and without inclusion of macroalgae in the Oosterschelde. In the Westerschelde, macroalgae predominantly occur on hard substrates and hardly contribute to the diet of macrofauna present on these structures (Riera et al. 2004). Therefore, their contribution to the diet of macrofauna on tidal flats is expected to be even lower. δ 13 C and δ 15 N values of terrestrial organic matter in the water column near Antwerp in April, May, and June were added as potential food source in the Westerschelde, as the estuary received large amounts of terrestrial organic material throughout the year (Van den Meersche et al. 2009). In the Oosterschelde, the input of terrestrial organic material is low (Nienhuis and Smaal 1994). It has been demonstrated that the contribution of saltmarsh plants to the diet of macrofauna on mudflats in front of a saltmarsh is low (Galván et al. 2008). Therefore, saltmarsh plants were not included as possible food source in the dual stable isotope analysis. We examined visually whether δ 13 C and δ 15 N values of macrofauna were in between values of δ 13 C and δ 15 N values, to check whether all possible food sources were taken into account. The Bayesian mixing model accounts for uncertainty associated with multiple resources, fractionation values and isotope signatures (Parnell et al. 2013). The analysis was performed in the R package MixSIAR (v3.1.7; Stock and Semmens 2013). In the isotope mixing model, we assumed an isotopic shift between diet and consumer of C of 0.3 AE 0.14‰ and of N of 2.2 AE 0.3‰ following McCutchan et al. (2003). Isotopic shift of N is known to vary among species, but an average isotopic shift of N of + 2.2 is supported by a number of studies measuring isotopic shift of N in marine organisms (Macko et al. 1982;Dittel et al. 1997). As we were interested in the total grazing pressure on microphytobenthos, we assumed that only one trophic level was present in the calculation of diet coefficients. The method uses Markov Chain Monte Carlo (MCMC) algorithms to calculate a posterior probability distribution of the proportion of each resource in the diet of a macrofaunal species. The default parameters for MCMC were used ("very long" option, with three parallel chains, 300,000 iterations for each chain). We used the Gelman-Rubin and Geweke Daggers et al. Diet and grazing pressure of macrofauna diagnostics to test whether the MCMC converged on the posterior distributions for all calculated diet contributions.
Estimation of macrofaunal grazing pressure Grazing pressure of the macrofaunal community (Grazing com ) on microphytobenthos, phytoplankton and terrestrial organic material was calculated for each station using a simple energy budget model for each species that was based on an allometric relation for respiration (Soetaert and Van Oevelen 2009). Mahaut et al. (1995) found the following allometric relation for biomassspecific respiration (d −1 ): 0.0174 x biomass spec -0.156 , in which biomass spec is the individual mass. The total C demand for a species is calculated by multiplying the biomass-specific respiration by the areal biomass (mg C m −2 , derived from macrofauna samples of NIOZ/Rijkswaterstaat) and subsequently dividing by the assimilation efficiency (AsE), to include the non-assimilated matter in the C demand, and by (1-NGE), in which NGE is net growth efficiency, to include biomass production in the C demand (Soetaert and Van Oevelen 2009). The microphytobenthos grazing for each species is subsequently calculated by multiplying the total C demand by the microphytobenthos diet fraction (i.e., the results from the isotope data). In summary, Assimilation efficiencies were derived from the review of AsE and NGE values for macrofauna in shallow marine ecosystems provided in Stratmann (2018), which matched for 35% on genus level, 39% on family level and for 48% on order level. For remaining species, the median value of 0.55 was used. The NGE was taken from the same review, which matched for 48% of the biomass on genus level, 48% on family level and 49% on order level. For remaining species, the median of 0.54 was used. Grazing pressure of the total macrofaunal community and of each feeding type was calculated for each station. When a species was not sampled at a particular station for stable isotope analysis, values reported by Herman et al. (2000) (24% of the biomass) or the median microphytobenthos dependency for bivalves vs. non-bivalves were used (4% of the biomass). A median value was calculated per estuary, if the food dependency differed significantly between estuaries (Supporting Information, Table S1).

Statistical data analysis
ANOVA (ANOVA and HSD posthoc Tukey test, significance level p = 0.05) was conducted on Chl a concentrations and microphytobenthos production using "estuary" (with levels Oosterschelde/Westerschelde), "elevation" (high/low), and "season" (early spring/late spring) as fixed factors and "site" as random factor. In addition, a possible interaction effect of "estuary" and "elevation" on Chl a concentrations was taken into account. We tested whether a linear correlation was present between elevation and Chl a concentrations using the Pearson product-moment correlation coefficient. ANOVA (ANOVA) was conducted on macrofaunal biomass to test for differences between estuaries, where each "site" within each estuary was considered a random factor. Species having a high contribution to the total grazing pressure were listed per site (Supporting Information, Table S3).
Differences in the proportion of microphytobenthos in the diet of individual species and the total of sampled individuals due to the categorical predictors "estuary" (with levels Oosterschelde and Westerschelde), "elevation" (with levels high and low), possible interaction effects between "estuary" and "elevation" and the random factor "site" were tested with an ANOVA followed by a HSD posthoc Tukey test (significance level p = 0.05) when significant.
Variance in grazing pressure of the macrofaunal community on food sources was tested using an ANOVA for the factors "estuary," "food source" (with levels microphytobenthos, phytoplankton, macroalgae and terrestrial organic matter) and the random factor "site," followed by a Tukey's test when significant. Linear regression tested the ability of microphytobenthic production in early spring to predict grazing pressure of (1) the total macrofaunal community on microphytobenthos and (2) grazing by macrofauna of separate feeding types on microphytobenthos.

Site characteristics Mud content and microphytobenthic biomass
The Oosterschelde contains relatively sandy tidal flats, while sediments of tidal flats in the Westerschelde have a highly varying mud content ( Supporting Information, Fig. S1). The microphytobenthic biomass, measured as Chl a, was higher in high than in low elevation plots (ANOVA, p = 0.0003, F 1,474 = 13.14, n = 486; HSD Tukey, p < 0.01; Fig. 2). Furthermore, microphytobenthic biomass differed significantly between early spring and late spring (ANOVA, F 1,474 = 6.06, p = 0.014, n = 486), although the HSD Tukey test revealed no significantly different means (HSD Tukey, p = 0.06). A significant interaction effect of elevation and estuary on Chl a concentrations was observed (F 1,474 = 5.75, p = 0.017, n = 486); Chl a concentrations were lower in the low lying intertidal areas of the Westerschelde than in the other areas (HSD Tukey, p < 0.01).
In late spring, no fluorescence signal above detection limit was measured at all sites in the Oosterschelde and at two sites in the Westerschelde (Molenplaat, MO, and Rilland, RI; Fig. 3).

Macrofaunal biomass and community composition
Macrofaunal biomass sampled by NIOZ/Rijkswaterstaat was on average 18,373 mg AFDW m −2 (σ = 33,570) and did not differ significantly between tidal basins (ANOVA, F 1,6 = 0.02, p = 0.88, n = 73; Fig. 4). Cerastoderma edule formed an important part of the total biomass comprising 8%-43% at 5 out of 9 sites (Supporting Information, Table S3). The lugworm Arenicola marina formed a significant part of the total macrofaunal biomass in the Oosterschelde and constituted 31%-38% of the total macrofaunal biomass at these sites. At 4 sites in the Westerschelde, the capitellid worm Heteromastus filiformis was among the three species with the highest total biomass relative to other species present, comprising 27%-51% of the total biomass. An overview of the actual biomass of the species is presented in Fig. 5.

Dietary proportions of macrofauna δ 13 C and δ 15 N values of food sources and macrofauna
At most sites (Rattekaai, Paulina, Middelplaat, Hellegat, Molenplaat, Waarde, and Rilland) the majority of macrofauna samples collected in late spring had an isotopic signature in between those of possible food sources (Fig. 6). At Viane, Waarde, and Rilland, the δ 15 N of some macrofaunal species was slightly heavier than available food sources. At sites in the Westerschelde, the signatures of the majority of macrofaunal species were in between those of microphytobenthos and phytoplankton and not similar to terrestrial organic matter.

Stable isotope mixing model
Including macroalgae as possible contribution to the diet of macrofauna resulted in unrealistically high proportions of macroalgae in the diet of macrofauna for the Oosterschelde (Supporting Information, Fig. S2). Observed macroalgae cover    (%) was only < 1% at Dortsman and Viane and AE 4% at Rattekaai (Supporting Information, Table S4). Therefore, macroalgae were excluded and the diet composition of macrofauna in the Oosterschelde was recalculated using microphytobenthos and phytoplankton as only possible food sources.
The Gelman-Rubin diagnostic gave acceptable (< 1.05) values for the Bayesian mixing model fit in all cases (Supporting Information, Table S5), indicating that each Markov chain converged to stable means and the same dietary proportions were obtained from each replicate Markov chain. The Geweke diagnostic exceeded acceptable values more frequently (in 56% of cases the diagnostic is equal to or > 5%), indicating that target distributions were not yet reached within the first 10% of the chain. However, plots of the running means as function of iteration showed that in model runs where the Geweke diagnostic was exceeded, stable means were always reached.

Spatial variability in microphytobenthos dependence of macrofauna
Microphytobenthos were an important food source for the majority of macrofaunal genera (fraction in the Oosterschelde: 0.70 AE 0.20; Westerschelde: 0.61 AE 0.16, Fig. 7   The contribution of microphytobenthos to the diet of macrofaunal genera did not differ between estuaries for the total of sampled individuals (ANOVA, "Estuary": F 1,5 = 5.7, p = 0.06, n = 102; Supporting Information, Table S1). |However, the contribution of microphytobenthos was higher in the Oosterschelde than in the Westerschelde for the sand digger shrimp Bathyporeia (HSD Tukey, p < 0.005), the balthic tellin Macoma (HSD Tukey, p < 0.0005), and the mudsnail Peringia (HSD Tukey, p < 0.05). No difference in the contribution of microphytobenthos to the diet of all sampled individuals or individual macrofaunal species was observed between high and low stations (Supporting Information, Table S1).

Comparison of faunal grazing on microphytobenthos with microphytobenthos production
Total grazing pressure on microphytobenthos was of the same order of magnitude as microphytobenthos production in early spring ( Fig. 9a; Oosterschelde: microphytobenthos production,μ = 235 mg C m −2 d −1 , microphytobenthos grazing,median = 312 mg C m −2 d −1 ; Westerschelde = microphytobenthos production,μ : 435 mg C m −2 d −1 , microphytobenthos grazing,median = 290 mg C m −2 d −1 ). At a number of stations, grazing pressure on microphytobenthos exceeded local microphytobenthic production in early spring (Fig. 9a,b, points above the 1:1 line). In most cases, grazing pressure was a factor 1.3 to 10 higher than local microphytobenthic production. There are, however, two stations where grazing pressure is several orders of magnitude higher than local production: Viane (Oosterschelde) and Paulinapolder (Westerschelde) (Fig. 9b). The high grazing pressure at these stations can be explained by a high biomass of suspension feeders present (Fig. 9c), that partially fed Fig. 7. Diet contribution (fraction) of microphytobenthos, phytoplankton, and terrestrial organic material to sampled macrofaunal genera in the Oosterschelde and Westerschelde. The known feeding group of each genus is indicated as SF = suspension feeders, DF = deposit feeders, SDF = surface deposit feeders, O = omnivores. Hediste and the isopod Cyathura are not displayed for the Oosterschelde and the bristleworm Tharyx and the polychaete Spionidae are not displayed, as only a small number of individuals of those genera was collected there. Boxes represent the first quantile, median, and third quantile. Whiskers extend to the largest vs. lowest value no further than 1.5 times the interquartile range. on microphytobenthos. Microphytobenthic production was exceptionally high (> 1600 mg C m −2 d −1 ) at one station, which was associated with a high average Chl a concentration (486 mg m −2 ) and high average photosynthetic capacity (82 μmol photon m −2 s −1 ). Total grazing pressure of macrofauna on microphytobenthos could not be predicted from microphytobenthos production using a linear regression in the Oosterschelde (R 2 = 0.01, p = 0.81) or Westerschelde (R 2 = 0.10, p = 0.24).

Discussion
In this study, the diet composition of macrofauna was calculated using a dual stable isotope model in contrasting tidal basins and elevations. Microphytobenthic production in early spring was linked to macrofaunal grazing pressure in summer/ autumn. Microphytobenthos was the main ultimate food source for macrofauna species in the two studied tidal systems. External material imported through rivers and streams contributed only marginally to the diet of the studied macrofaunal genera. This finding is in line with earlier stable isotope studies conducted in temperate, subtropical and tropical tidal basins (see below). The contribution of microphytobenthos in the diet of macrofauna varied between tidal systems for a number of species, but did not vary as function of elevation. Very few studies have been conducted on spatial variability in the contribution of microphytobenthos in the diet of macrofauna (Christianen et al. 2017). To our knowledge, this is the first study where the diet contribution of microphytobenthos was compared among tidal systems or linked to tidal elevation. We demonstrated that grazing pressure on microphytobenthos is of the same order of magnitude as microphytobenthic production, which likely explains the frequently observed drop in microphytobenthic biomass in early summer in temperate systems. Spatial variability in microphytobenthic production in early spring did not explain spatial variation in macrofaunal grazing pressure in summer/ autumn, which may be associated with the time difference between sampling of macrofauna and measurement of microphytobenthic production. Assuming that macrofauna communities are relatively stable over time, our results suggest that consumer-resource interactions are coupled on a larger spatial scale (i.e., mesoscale, ≈ 10 to 100 km) rather than the fine (millimeter to meter) scale.

Isotope composition of food sources and macrofauna
The relative contribution of sources to the diet of organisms can be studied using bulk isotope signatures (Yokoyama and Ishihi 2003) or compound-specific stable isotope analysis (Oxtoby et al. 2016). Bulk isotopic signatures are suitable to obtain a quantitative estimate of dietary proportions, although limitations include isotopic routing and fractionation after consumption (Federer et al. 2010). Compound-specific analyses allow the use of biomarkers with little structural modification, such as polyunsaturated fatty acids, but require detailed knowledge on tissue composition to estimate proportions of sources in consumers (Wolf et al. 2009). The δ 13 C of sampled macroalgae was in the same range as the δ 13 C of microphytobenthos and phytoplankton (Supporting Information, Tables S7 and S8) and could contribute to the diet of macrofaunal genera when macroalgae cover is high. Several studies have quantified isotopic signatures of food sources for benthos on tidal flats in the Oosterschelde and Westerschelde Riera et al. 2000;Moens et al. 2005;Van den Meersche et al. 2009) and in other intertidal areas worldwide (Kang et al. 2003). Moens et al. (2005) demonstrated that carbon isotopic signatures of the food sources SPOM and microphytobenthos are generally narrowly constrained (δ 13 C μ, microphytobenthos : −15.7 to − 14.3; δ 13 C μ, SPOM : −24 to − 20) in the Westerschelde. The range of the δ 13 C of microphytobenthos is somewhat larger in the current study. (δ 13 C μ, microphytobenthos : −18.2 to − 12.2; δ 13 C μ, SPOM : −21.7 to − 25.8). The δ 13 C of microphytobenthos in our study was negatively correlated with mud content (Pearson's r = 0.53). This could be related to a general lowering of the δ 13 C of DIC in the upstream direction in estuaries as function of salinity (Chanton and Lewis 1999) and a concomitant increase in mud content. Generally, the δ 13 C of DIC is lighter in riverine waters than in marine waters (e.g., Chanton and Lewis 1999). In our study, we also observed relatively low δ 13 C values of microphytobenthos at a relatively silty location close to the polyhaline estuary mouth (Paulinapolder: δ 13 C μ, microphytobenthos : −18.2 to − 17.7). It is likely that more estuarine phytoplankton deposits at such sheltered sites. Part of the PLFAs used in this study as microphytobenthos markers could therefore be derived from this more 13 C depleted phytoplankton (δ 13 C of 20:5ω3 characteristic for diatoms, sampled from surface water:~20-25‰; Boschker et al. 2005), explaining the lighter signatures. However, the species composition of algae in the Westerschelde intertidal flats is dominated by microphytobenthos (Sabbe and Vyverman (1991), rather than plankton species. This suggests that the relatively low microphytobenthos δ 13 C ratios at relatively muddy sites are more likely due to site dependent differences in growth conditions on isotopically depleted DIC.
The range of nitrogen isotopic signatures of microphytobenthos reported in the current study (δ 15 N μ, microphytobenthos : 7.7 to 21.2) are similar to the values reported by Moens et al. (2005) (δ 15 N μ, microphytobenthos : 7.4 to 17.7; δ 15 N μ, SPOM : 7 to 15) for the Westerschelde, whereas a smaller range of δ 15 N values of SPOM was found (δ 15 N μ, SPOM : 8.4 to 9.9). The generally high range of nitrogen isotopic signatures of microphytobenthos and SPOM found in estuaries is associated with assimilation of the relatively heavy NH 4 + by bacteria and microalgae, likely resulting from isotopic fractionation by nitrifiers and/or fractionation by phytoplankton uptake (Mariotti et al. 1984;Riera et al. 2000;Moens et al. 2005). However, the relative importance of nitrification and, therefore, availability of heavy NH 4 + for microalgae may vary among sites (Moens et al. 2005) and seasons (Rysgaard et al. 1995).
Furthermore, nitrogen in suspended matter of marine origin is known to have a higher nitrogen isotopic signature than suspended matter of terrestrial origin (δ 15 N μ, marine OM : 8 AE 1.8; δ 15 N μ, terrestrial OM : 1.5 AE 0.2; Mariotti et al. 1984). Especially in winter, when autochthonous primary production is low, terrestrial organic material may form a significant part of SPOM of estuarine systems. In summer, SPOM in the Westerschelde is dominated (> 50%) by autochtonous phytoplankton (Mariotti et al. 1984). The observed spatial variability in nitrogen isotopic signatures of SPOM in the Oosterschelde and Westerschelde emphasize the importance of sampling and determining isotopic signatures of SPOM locally.
The carbon and nitrogen isotopic signatures of macrobenthic genera used to determine the diet composition are determined by (1) the δ 13 C and δ 15 N of consumed food sources, (2) changes in the δ 13 C and δ 15 N due to differential digestion or fractionation during metabolic and assimilation processes (isotopic shift or trophic fractionation) and (3) the position of the consumer in the food chain (trophic level).
Here we have adopted an average isotopic shift of 0.3 AE 0.14‰ for carbon and of 2.2 AE 0.3‰ for nitrogen per trophic level, as proposed by McCutchan et al. (2003). However, isotopic shift varies among individual consumers (McCutchan et al. 2003) and have, to our knowledge, not been quantified for the specific intertidal macrobenthic genera considered in this study. Furthermore, we assumed that the macrobenthic genera considered in this study are separated from the considered food sources by one trophic level, which is the generally accepted position of macro-invertebrates in the food chain of intertidal ecosystems (Kuwae et al. 2012). By correcting for shifts in δ 15 N after food ingestion, the calculated diet composition values represent an "ultimate food source" diet composition. Likewise, the calculated grazing pressure represents the grazing pressure that is ultimately present on the considered food sources. It has been demonstrated that primary food sources may be transferred to macroinvertebrates via bacteria or meiofauna . As a consequence, the number of trophic levels may be underestimated and the grazing pressure on microphytobenthos and the contribution of the considered food sources to the diet of macrofauna may be slightly overestimated. Lastly, it should be noted that diet composition or even trophic level may vary over time, due to seasonal changes in food availability.

Microphytobenthos as food source for macrofauna
We showed that microphytobenthos is the dominant food source for the majority of abundant macrofaunal species in two contrasting tidal systems, including predators, suspension feeders and facultative suspension feeder/grazers. This is in accordance with earlier findings on tidal flats in the Westerschelde , in other regions with a temperate climate such as the Wadden Sea, The Netherlands (Christianen et al. 2017), the Plum Island estuary, Massachusetts (Galván et al. 2008), but also in regions with a subtropical (e.g., Kang et al. 2003) or tropical climate (Kon et al. 2015). The contribution of microphytobenthos in the diet of suspension feeders (24%-44% for Cerastoderma) was somewhat higher in this study than previously reported values (17%, Herman et al. 2000;< 5%, Christianen et al. 2017). Our study suggests that it is not correct to assume a priori that suspension feeders (almost) exclusively depend on phytoplankton. Instead, resuspended microphytobenthos may contribute significantly to the diet of suspension feeders (cf. Fig. 6).

Spatial variability in the diet of macrofauna
In our study, the contribution of microphytobenthos to the diet of macrofauna did not differ between estuaries for all sampled individuals combined. However, the contribution of microphytobenthos in the diet of the genera Bathyporeia sp., Macoma balthica and Peringia ulvae was higher in the Oosterschelde than in the Westerschelde. Our macrofaunal community dataset shows that in the Westerschelde, Bathyporeia pilosa is the most dominant species of the Bathyporeia genus, whereas in the Oosterschelde Bathyporeia pelagica and Bathyporeia elegans occur. All three species are surface deposit feeders, indicating that the composition of organic material differs between the tidal systems. This may also explain the larger contribution of microphytobenthos in the diet of Peringia ulvae (a surface grazer) in the Oosterschelde than in the Westerschelde. Macoma is well known to be a facultative surface deposit feeder and suspension feeder (Rossi et al. 2004). As phytoplankton production in spring is generally lower in the Oosterschelde than in the Westerschelde (Kromkamp and Peene 2005), this may explain the higher proportion of microphytobenthos in the diet of Macoma in the Oosterschelde. We did not observe size dependent (ontogenetic) differences in microphytobenthos dependence of Macoma as described by Rossi et al. (2004) and Herman et al. (2000). However, the possibility that a systematic error in the estuary-specific offset used to obtain δ 15 N values of microphytobenthos from the δ 15 N of the bulk sediment resulted in significant differences in the diet of macrofauna between estuaries cannot be ruled out.The contribution of microphytobenthos in the diet of macrofauna did not vary as function of elevation in the tidal zone. This finding may be counter intuitive, as microphytobenthos biomass is generally positively correlated with emersion duration (Van der Wal et al. 2010) and the studied species generally occurred at high and low elevations (Supporting Information, Table S6). However, the lack of a significant difference in the contribution of microphytobenthos in the diet of macrofauna is in line with the finding that microphytobenthos production did not differ between high and low plots, indicating that there is no difference in the amount of microphytobenthos available for macrofaunal species between high and low tidal flats. Although sample size did not allow analysis for the factor "site," some variation in the microphytobenthos diet coefficient appeared to be site-dependent (Supporting Information, Table S8). This may be associated with local differences in the availability of different food sources. For example, a sheltered tidal flat may contain a relatively large amount of deposited terrestrial organic material.

Relation between microphytobenthos production and macrofaunal grazing
We showed that grazing pressure of macrofauna on microphytobenthos in summer/autumn (i.e., calculated using the density and biomass of the macrofauna community sampled by NIOZ/Rijkswaterstaat from August-October) is in the same order of magnitude as microphytobenthic production in early spring. This suggests that microphytobenthos is strongly subject to top down control by macrofauna, as microphytobenthic production is generally highest in spring. This finding is consistent with the results of defaunation experiments in intertidal areas (Weerman et al. 2011), which showed a rapid increase in microphytobenthic biomass after removal of macrofauna. Furthermore, Savelli et al. (2018) modeled grazing of Peringia ulvae on microphytobenthos on a mudflat in Northwest France, where grazing was regulated by microphytobenthos availability and mud surface temperature, and found that grazing rates exceeded microphytobenthic production rates in spring. Labeling experiments indicated a low transfer of carbon from microphytobenthos to macrofauna Middelburg et al. 2000). However, such experiments only provide information on carbon transfer on a time scale of a few days. As carbon originating from microphytobenthos may initially be taken up by bacteria or meiofauna , carbon transfer from microphytobenthos to macrofauna may occur on longer time scales.
Microphytobenthic production in early spring did not explain local spatial variability in total macrofaunal grazing pressure in summer/autumn (August-October). This may partly be associated with the presence of large suspension feeders, which feed on non-locally produced microphytobenthos from the water column. As suspension feeders form a significant part of the total macrofaunal biomass per unit area, the link between microphytobenthos production and total macrofaunal grazing per unit area is weak. In addition, microphytobenthic production was not linked to grazing pressure by macrofauna of separate feeding types. The macrofauna community composition is highly variable among stations, leading to large variation in the total biomass of different feeding types among stations. Microphytobenthic production was not always measured at the exact same location as where samples for the determination of the macrofaunal community composition were taken (within 300 m), and there was a time difference between measurement of microphytobenthic production (2015) and sampling for the macrofauna community composition (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). This may have caused non-linearity in the link between available microphytobenthic production and the grazing pressure of a macrofaunal feeding type on microphytobenthos. Earlier studies have demonstrated that the large-scale spatial distribution of macrofaunal species is relatively stable in the Westerschelde (Ysebaert et al. 2003), albeit with variation (Ysebaert and Herman 2002). Likewise, the spatial patterns in microphytobenthos biomass are relatively consistent over time, superimposed on seasonal, and year-to-year variability ( Van der Wal et al. 2010).
It should be noted that the presented microphytobenthic primary production rates are gross primary production rates, while the actual carbon availability for the food web is represented by net primary production rates. Therefore, the presented microphytobenthic production rates are not 1:1 comparable to the presented macrofaunal grazing rates, i.e., not all of the carbon fixed by the microphytobenthos is available for consumption. In a review, Langdon (1993) showed that the biomass-specific dark respiration of diatoms is a fixed fraction of gross primary production of on average 0.06% AE 0.01%, depending on algal growth rates. Net primary production rates approach gross primary production rates under low specific growth rates (Halsey et al. 2010).
Our results confirm the key role of microphytobenthos on intertidal flats in the diet of macrofauna, which are in turn eaten by secondary consumers such as (wading) birds and fish. This highlights the pressing need to preserve these intertidal flats, as decreasing availability of microphytobenthos may have cascading effects up the estuarine food web. Furthermore, we demonstrated that the diet composition of macrofauna likely differs among estuaries, which suggests that estuary-specific food web modeling can support management of these ecosystems.