Shorebirds Affect Ecosystem Functioning on an Intertidal Mudflat

Ecosystem functioning and services have provided a rationale for conservation over the past decades. Intertidal muddy sediments, and the microphytobenthic biofilms that inhabit them, perform crucial ecosystem functions including erosion protection, nutrient cycling and carbon sequestration. It has been suggested that predation on sediment macrofauna by shorebirds may impact biofilms, and shorebirds are known to consume biofilm, potentially causing significant top-down effects on mudflat ecosystem functioning. We carried out an exclusion experiment on the Colne Estuary, Essex, to examine whether shorebird presence significantly affects sediment erodibility measured with a Cohesive Strength Meter (CSM) and microphytobenthos biomass measured using PAM fluorescence (Fo) and chlorophyll a content. We also tested for treatment effects on sediment-water nutrient fluxes [nitrate, nitrite, ammonia, phosphate and dissolved organic carbon (DOC)] during periods of both dark and light incubation. Excluding shorebirds caused statistically significant changes in regulating and provisioning ecosystem functions, including mudflat erodibility and nutrient fluxes. The presence of shorebirds lowered the sediment critical erosion threshold τcr, reduced nitrate fluxes into the sediment under illumination, lowered nitrate efflux, and reduced phosphate uptake, compared to sediments where birds were excluded. There were no significant differences in macrofauna community composition within the sediment between treatments after 45 days of bird exclusion, suggesting a direct link between shorebird presence or absence and the significant differences in biofilm-related variables. This study introduces previously unknown effects of shorebird presence on ecosystem functions within this system and highlights an area of shorebird science that could aid joint conservation and human provisioning action.


INTRODUCTION
Ecosystem functioning and ecosystem services have provided a rationale for conservation over the past decades (Cabello et al., 2012). Intertidal mudflat ecosystem functions include nutrient cycling, erosion protection and carbon sequestration, which mediate associated services (Foster et al., 2013). Intertidal flats provide natural 'soft' coastal erosion defense by reducing wave energy, lowering water velocities and thereby shear stress on the estuary bed (Spalding et al., 2014). Benthic microalgae [microphytobenthos (MPB)] form complex matrices of cells, sediments and extra polymeric substances (EPS) (Underwood and Paterson, 2003). These biofilms have a stabilizing effect on surface sediments, reducing erodibilty and aiding in the accumulation of particles and microbes (Gerbersdorf and Wieprecht, 2015). Estuarine sediments and biofilms are central components in estuarine nutrient cycles, ultimately affecting fluxes of these nutrients between land and sea (Thornton et al., 2007;Nedwell et al., 2016). Organic compounds are recycled and remineralized within sediments, particularly in coastal marine areas where nitrogen and phosphorous loads can be very high (Correll et al., 1992;Hochard et al., 2010). Nitrogen loading into marine systems can lead to eutrophication and decline in water quality, making its source and removal pathways of high interest (Burgin and Hamilton, 2007) and changes in nutrient loads can impact benthic communities (Culhane et al., 2019). MPB mediate fluxes of NO 3 − , NO 2 − , PO 4 3− and NH 4 + between the water column and sediment layers (Sundback et al., 1991;Correll et al., 1992;Feuillet-Gerard et al., 1997), contributing to this process either by direct uptake/release or by altering oxygen concentration (Sundback and Graneli, 1988). Dissolved organic carbon (DOC) may also provide an important part of both global and coastal carbon sinks (Maher and Eyre, 2010;Legge et al., 2020), making effects on DOC fluxes in this environment relevant to anthropogenic climate change effects and mitigation (McKinley et al., 2016).
Mud and sand flats are essential habitats for the survival of resident and migratory overwintering shorebirds (Burton et al., 2006), which feed primarily upon infaunal and epifaunal invertebrates (Bowgen et al., 2015). Some small sandpiper species Calidris spp. also directly consume biofilm during, or in preparation for, migration (Kuwae et al., 2008;Jardine et al., 2015). Grazing of MPB and bioturbation by macrofauna can lead to alterations in sediment erodibility and other ecosystem functions (de Deckere et al., 2001;Hale et al., 2019). This poses questions regarding the effect of biofilm removal and bioturbation by shorebirds (Mathot et al., 2018), which may have significant knock-on effects altering ecosystem functions.
Research suggests that shorebirds could have significant direct and/or indirect effects on ecosystem function, e.g., via the impacts of foraging on macrofauna and/or biofilm or disturbance and reworking of sediment (Orvain et al., 2014b;Mathot et al., 2018). In the Bay Of Fundy (BOF), semipalmated sandpipers Calidris pusilla appeared to cause an ecological cascade effect by reducing densities of their mud shrimp prey Corophium volutator, which caused biofilm proliferation, leading to an increase in sediment stability (Daborn et al., 1993). However, subsequent research in the BOF has not indicated a trophic cascade effect, possibly due to compensatory interactions by macrofauna (Hamilton et al., 2006;Cheverie et al., 2014). Trophic webs and ecosystem functioning were compared in the Marenne-Oleron Bay, France, indicating that estuarine trophic webs including shorebirds have enhanced primary productivity through increased nutrient cycling (Saint-Beat et al., 2013). Despite evidence that estuarine shorebirds may significantly alter ecosystem functioning, the majority of shorebird research has an ornithological focus and potential top down effects on ecosystem functions such as erosion defense and nutrient cycling have not yet been experimentally tested (Mathot et al., 2018). The ecology of intertidal sediments is complex, compensatory interactions can mask effects (Hamilton et al., 2006), including trophic cascades (Fahimipour et al., 2017). Manipulative experiments are a valuable tool, to be utilized alongside 'natural' or 'observational' experiments to assess possible ecological mechanisms behind processes observed at wider spatial or temporal scales (Rogers et al., 2012).
The Colne Estuary, Essex, United Kingdom is a complex of habitats featuring many sand and mudflats, protected internationally under The Conservation of Habitats and Species (Amendment) (EU Exit) Regulations 2019, for supporting over 30,000 shorebirds. Our study site within the Colne Estuary, the Fingringhoe Wick Site of Special Scientific Interest (SSSI), was a location for the six year Coastal Biodiversity and Ecosystem Service Sustainability research program (CBESS), which provides key background information on the biotic and abiotic characteristics of the site.
Changes in community composition and mudflat characteristics can be rapid, occurring over months (Sahan et al., 2007;Rosa et al., 2008;Murphy and Tolhurst, 2009) weeks (Daborn et al., 1993;Hamilton et al., 2006), days (de Deckere et al., 2001;Tolhurst et al., 2008) and even hours (Tolhurst et al., 2006a,b). We designed and carried out a two month field exclusion experiment, supplemented by laboratory measurements, to investigate shorebird effects on two ecosystem functions, namely erosion protection (using a measure of sediment erodibility as a proxy) and nutrient cycling (including nitrate, nitrite, ammonia, phosphate and DOC). We tested three hypotheses: (1) surface biofilm biomass would be significantly altered in the presence of shorebirds, (2) sediment erodibility would be significantly altered in the presence of shorebirds and (3) nutrient fluxes between the sediment and water column would be significantly different between treatments (shorebird presence and absence) with flux direction and magnitude for different nutrient species increasing with greater MPB biomass.

Description of Study Site
Fieldwork was undertaken between 20 January and 03 April 2017 on the mudflat adjacent to Geedon Saltings at Fingringhoe Wick Essex Wildlife Trust Nature Reserve, Essex, United Kingdom (grid reference TM 05065 19030). This time period covered the peak overwintering and start of the migratory periods for shorebirds in the East of England. The study location comprised an area of mudflat approximately 400 m 2 situated on the upper shore. Observations during 2016 noted flocks of dunlin Calidris alpina and knot Calidris canutus, and scattered individual redshank Tringa totanus and gray plover Pluvialis squatarola foraging at the study site on receding and incoming tides. The study location was set within a larger area of estuarine mudflat, approximately 130,000 m 2 of which could be visually surveyed for shorebird activity from a fixed point (Geedon hide).
Previous CBESS studies showed that during winter, sediment at the site is mostly silts and clays, with a very low proportion of sand (maximum 'very fine sand' content in a sample was 6.5%; coarser sand contents were all lower than this), with sediment particle size at the site predominantly <63 µm (mean 95.9 % ± 0.3 SE). Mean D 50 = 6.9 µm ± 0.2 SE; Mean D 16 = 1.9 µm ± 0.04 SE; D 95 = 66.9 µm ± 13.2 SE (D x = particle diameter representing the x% cumulative percentile) (Wood et al., 2015). Mean percentage surface sediment water content at the site is 62.3 ± 0.4 SE (Maunder and Paterson, 2015). This site lies within the polyhaline section of the estuary, with salinity ranging from 18-30, depending on freshwater flow conditions, with lower salinity during winter (Nedwell et al., 2016).
CBESS research also included sampling of fauna within the Colne estuary, demonstrating that during winter fish were absent, with only Ctenophores recorded during Fyke netting (Wood et al., 2015). Macrofauna recorded during winter CBESS research included ragworm Hediste diversicolor, mud snail Peringia ulvae, Baltic clam Macoma balthica and nematodes across a total of 22 quadrat sites, in which three samples were taken at each (Wood et al., 2015). A year-long fish sampling study carried out at two different locations along the estuary where our experiment was undertaken, found that fish were absent at all sampled sites during January, and absent from three out of five sites during February (Green et al., 2009). Where fish were present at two sites during February, total abundance (fish 100 m −2 ) was approximately 2, and less than 1 during March (Green et al., 2009).

Experiment Design
The manipulative experiment was set up on 20 January 2017 (day 0). The experimental layout was a randomized design of 20 spatial plots (Figure 1), each 1 m x 1 m, allocated to two treatment levels; control (shorebirds present in open unmanipulated plots) and exclosure (shorebirds absent), with n = 10 replicates of each treatment. Previous work in the estuary showed that spatial variability in biofilm abundance is greatest at the fine scale and small at the meter scale (Taylor et al., 2013;Nedwell et al., 2016), therefore a completely randomized design was employed to maximize statistical power of the experiment. Exclosures were bamboo frames, approximately 30 cm in height, covered on all sides (including the top) by opaque 'fruit-cage' bird exclusion netting (plastic mono-thread) with a 2 cm aperture. Exclosures prevented access to the sediment by birds, but allowed access to infauna and small fish (<2 cm width). All plots were at least three meters apart, to allow sampling from all sides and prevent plots unduly influencing each other. Exclosure and control plots were unpaired and separated by similar distances, with treatments arranged sequentially to reduce the potential for spatial bias. The exact locations of plots were selected to represent the heterogeneity within the wider mudflat. No scouring or bite marks indicating the presence of larger fish (Eggold and Motta, 1992) were found within any plots during the experiment. Plots were arranged parallel to the tide line (within a minute of immersion/emersion time of one another). Plots were situated on the upper shore, where shorebirds spend most time foraging due to the longer emersion time (Granadeiro et al., 2006). Camera footage (see below) and direct observation recorded no events of birds standing on exclosures (behavior which may otherwise have caused input of droppings into exclosure absence plots as well as control presence plots) (Schrama et al., 2013;Jauffrais et al., 2015).

Assessment of Possible Experimental Artifacts
To test the effect of the exclosures on the water flow within the study area, a 'plaster ball dissolution test' was carried out on days 17 and 18 (Cheverie et al., 2014). No significant difference was detected between plaster dissolution rates in control plots and exclosure plots (t = -1.057; df = 8; p = 0.322), demonstrating that our exclosures had no significant effects on tidal water flows in the vicinity of the mudflat surface.
Exclosure shading tests were carried out after the experiment to prevent additional mudflat disturbance, during a sunny day (cloud cover < 10%), hence resulting in an estimation of shading at the higher end of the actual range during the study period. Shading effects on Photosynthetically Active Radiation (PAR) reaching the sediment surface in exclosures were small (9.9%), and of a similar level to that in other manipulative studies in this type of environment (Cheverie et al., 2014). Further information reinforcing this conclusion is given in the discussion.
A Go-Pro HERO 4 camera fitted with a Cam-Do Blink time lapse controller mounted within a Cam-Do Solar-X enclosure (Cam-Do Solutions, 2017) was deployed to monitor bird activity within the study area for four weeks (21 February 2017 to 21 March 2017). This was mounted on a vertical pole 3.5 m above the saltmarsh at grid ref: TM 05031 19032. The camera captured a still of the plots every five minutes during daylight hours. Although species identification was not possible using captured images, numbers within the field of view were used to broadly determine whether numbers of birds using the study area were consistent with those recorded during visual surveys.
Weather data were collected during the experimental period [peak wind speed (km h −1 ), daily precipitation (hours day −1 ) and peak temperature ( • C)], and plotted against biofilm biomass (F o ) and shorebird numbers to assess potential effects of these variables on the experiment, such as extreme weather events, which can have significant effects on shorebird activity (Sutherland et al., 2012) and mudflat characteristics (Tolhurst et al., 2006b;Fagherazzi et al., 2017;Hale et al., 2019). No extreme weather events occurred during the experiment and no evidence was found of a relationship between F o and daily precipitation (hours), peak temperature ( • C) and peak wind speed (km h −1 ) during the experiment (Figures 2A,B), although the potential for delayed responses has not been assessed. However, all plots were subject to the same weather and this is not considered to be a constraint to the experiment.

Response Variables
Between day variation in mudflat characteristics have been shown to be of greater significance than within day variation (Tolhurst and Chapman, 2005), therefore repeated measures of F o were made to compensate for this effect. Table 1 shows dates and days at which sampling events took place. On 20 January 2017, immediately following plot setup, 'day 0' minimum fluorescence (F o ) measurements were taken using a pulse amplitude modulated fluorometer (PAM, Walz, Effeltrich, Germany) to determine MPB biomass (Honeywill et al., 2002). MPB are key drivers of intertidal flat properties and processes (e.g., Murphy and Tolhurst, 2009), so to determine when the full sampling event would be most likely to detect any effects we monitored F o (as a proxy for MPB biomass) on days 3, 13 and 26, as a convenient indication of treatment effects, to determine when erodibility and nutrient flux variables should be measured and to confirm that early in the experiment there were no significant differences between treatments. F o was also measured on day 45 to evaluate the effect of shorebird presence/absence on MPB biomass and associated properties, and on day 64 to determine if trends continued. A subset of 6 exclosure and 6 control plots were measured on day 3 for a total of 60 F o measurements (n = 5 in each of the 12 plots); subsequently all plots were measured, for a total of 100 F o measurements (n = 5 in each of the 20 plots) on days 13, 26, 45, and 64 to investigate how surface MPB biomass responded to shorebird presence/absence over time.
Due to the large number of measurements required in each plot during a tidal cycle and considering the impact of dewatering during the tidal cycle (Maggi et al., 2013;Orvain et al., 2014a;Fagherazzi et al., 2017) a 5 min low light partial dark adaption treatment was used prior to each PAM measurement, which is a preferred method to conventional dark adaption for the measurement of minimum fluorescence as a proxy of MPB Frontiers in Marine Science | www.frontiersin.org biomass (Jesus et al., 2006b). Sampling was carried out during periods of clear weather with little wind and no rain, at least one hour after the tide had exposed the sampling area to allow initial drying of plots. A consistent low light sampling environment was achieved using plastic 40 mm (diameter) × 60 mm (length), cylindrical opaque dark adaption chambers with a 6 mm aperture hole at the top. This also enabled in-situ sampling with the PAM fluorometer without removal of the chamber. This reduced the variation in light intensity during the measuring period. To further eliminate potential effects of varying light intensity and sediment water content during sampling events, exclosure and control plot sampling was alternated. To minimize the effect of varying light intensity and phase of vertical migration between sampling events, sampling periods were timed to cover low tides peaking as close to midday as possible.
Our experience of the site is that variability at the meter scale is low (Redzuan, 2017). Additionally, the repeated F o sampling (described above) gives further confidence that plots were not significantly different at the beginning of the experiment. All in situ mudflat variables were measured on 06 March 2017, after 45 days of shorebird exclusion, to test the effect that a period of shorebird exclusion had on selected mudflat properties. Sampling included in-situ measurements of F o (as described above), in-situ sediment critical erosion threshold (τ cr ) using a Cohesive Strength Meter (CSM) (three measurements within six plots of each treatment, total 36 measurements) (Tolhurst et al., 1999;Vardy et al., 2007) and contact coring for analysis of chlorophyll a content (three measurements within seven plots of each treatment; total 42 measurements) (Honeywill et al., 2002). Flux cores (Perspex tubes of 0.1 m diameter and approximately 0.2 m in depth) were also collected (one from each plot) for laboratory analysis of nutrients and macrofauna.
Contact cores (surface ∼2 mm) were freeze dried in the dark and chlorophyll a extracted using cold methanol over 24 h, and measured spectrophotometrically, correcting for phaeopigments (Stal et al., 1984).
Flux cores were carefully returned to the laboratory within an hour of leaving the site and immersed in seawater from the site, within oxygenated and temperature and light controlled indoor mesocosms (Thornton et al., 1999). Rubber bungs were used to ensure equal headspace volume across cores. Cores were left submerged and open to settle overnight prior to sampling on the following day. Throughout headspace water sampling, Perspex lids were tightly fitted to prevent leakage. Magnetic stirrers maintained water flow over the sediment surface. On 07 and 08 March 2017 these were sampled for sediment-water biogeochemical fluxes of nitrate, nitrite, ammonia, phosphate and dissolved organic carbon (DOC). Headspace seawater samples were taken at the beginning and end of 2 h dark and light incubation periods. Cores were left for at least one hour to adjust to light levels prior to each incubation. Sampling was completed according to general methods described by Thornton et al. (1999). Flux measurements were repeated in both light and dark conditions, using 500W halogen 'daylight' lamps to provide 'lit' conditions (500 µmol m −2 s −1 PAR) and covering mesocosms with opaque Perspex covers to provide 'dark' conditions. Water samples were analyzed for their nutrient concentrations using a Seal AA3 segmented flow Nutrient Analyzer (SEAL Analytical Inc.).
Individual cores used for nutrient flux measurements were subsequently sieved (500 µm mesh) to retain macrofauna. Macrofauna were preserved in 95% ethanol and identified to species level (where possible) using a microscope, quantified and densities (m −2 ) calculated. Through data comparison with previous work at the site (Wood et al., 2015) we were confident that sufficient sampling had been undertaken to assess potential differences in community composition between shorebird presence and absence plots.
Bird surveys began on 03 January 2017 (−17 days) and were carried out at least every two weeks (see Table 1) to monitor the level and type of use of the study area by shorebirds. Monitoring began before the experimental setup to ensure current use of the study area by shorebirds and aid in deciding the best location for the experimental plots. Surveys were carried out using the 'look-see' methodology (Bibby et al., 2000), from a fixed location (Geedon Hide; TM 05081 19170). Surveys were undertaken for at least 2 h either side of low-tide, including as much of these timeframes as possible (four hours maximum) within daylight constraints. Particular care was taken to also include visual observation of the tideline crossing the plots wherever possible. Counts of species within the surrounding visible mudflat were taken every half hour. Continual observation of the study area was made, quantifying numbers and identifying species entering presence plots throughout the surveys. Equipment included a 20−60 × 82 telescope and 10 × 42 binoculars. No birds were recorded within or on the absence plots during any of the surveys. During F o measurements, shorebird tracks were noted within all presence plots at some point during the study, indicating use of all presence plots by shorebirds. No tracks were recorded in any absence plots at any point during the experiment.

Statistical Analysis
To evaluate the effects of shorebird presence and time (days) on biofilm biomass throughout the experimental period, we used a linear mixed-effects model (plot nested in treatment) to analyze F o data with plot as a random effect and time (day) and bird presence/absence as fixed effects. This model was run using NLME package in R version 4.0.
To evaluate the effect of shorebird presence/absence on MPB biomass and sediment erodibility, F o (days 3, 13, 26, 45, and 64), chlorophyll a (from surface 2 mm) (day 45) and critical erosion threshold (day 45) data were analyzed using a mixed model, two-way nested ANOVA design with (plot nested in treatment) with plot as a random factor and shorebird presence/absence as a fixed factor, using the GMAV (1997) statistical package (University of Sydney, Australia). Although baseline data were not collected, ANOVA detects differences between treatments over and above variability among individual plots (Underwood, 1997). To counteract the issue of multiple comparisons we used Bonferroni correction testing each hypothesis at a confidence level of 0.01 (0.05/5).
To evaluate the effect of shorebird presence/absence on nutrient flux (day 45), nutrient data were analyzed using a two-way orthogonal ANOVA design with dark/light incubation and shorebird presence/absence as fixed factors, using the GMAV (1997) statistical package (University of Sydney, Australia). Where Cochran's test was significant (ammonium and phosphate), data were normalized by rank transformation and the analysis repeated. We also used reversals in flux (for example an efflux from the sediment in the absence of shorebirds becoming an influx into the sediment in the presence of shorebirds) as an indication of changes suggesting 'ecologically significant' implications for ecosystem functioning.
To assess whether shorebird presence/absence had significantly altered macroinvertebrate community structure, day 45 taxa density was analyzed using R version 3.6.1 with vegan package. Non-metric multidimensional scaling (NMDS, Bray-Curtis dissimilarity, 20 restarts) was used to visualize differences in community structure at day 45 in two dimensions (Clarke, 1993). The MDS had a stress 0.037, therefore considered an adequate representation (Clarke, 1993). Analysis of similarities (ANOSIM) was also performed to test quantitatively for differences in community structure between shorebird presence and absence.
To assess the potential for biases associated with the exclosures, plaster ball dissolution (days 17 and 18) and shading effect (post experiment) data were also analyzed using a one-way orthogonal ANOVA, using the GMAV (1997) statistical package (University of Sydney, Australia).
To evaluate shorebird pressure on the mudflat, species count data were first converted into 'bird-days, ' by calculating the sum of the number of each shorebird species present on every count, multiplied by the number of days between that and the subsequent count (Gill et al., 2001;Lewis et al., 2014). This allowed comparison of shorebird pressure on the wider mudflat. Only species considered regular foragers on mudflats and recorded foraging on the surrounding mudflat were included in this analysis; for example lapwing Vanellus vanellus and golden plover Pluvialis apricaria were removed due to their high dependence, and almost exclusive foraging, on coastal grassland and arable fields (Mason and Macdonald, 1999). Furthermore, these species were recorded roosting on the mid to low shore only during low tides, further reducing the likelihood that they contributed to any effects within the upper shore study site. To compare mudflat variables with density of species recorded in presence plots, count numbers of such species were log 10 transformed and plotted over time with mean F o in shorebird presence and absence.

Microphytobenthic Biomass
Results of the linear mixed effects model show a highly significant difference in F o (measure of MPB surface chlorophyll a) between shorebird presence and absence, with F o higher in the bird exclosure treatments. There was no significant effects of time (days) or interaction between treatment with time ( Table 3).
F o initially increased in shorebird presence and absence plots, increasing more rapidly in absence plots, peaking on day 26 before decreasing (Figure 2A). On day 3, there was no significant difference in F o between shorebird presence and absence plots, but on day 13 this difference had become significant. The largest difference was measured on day 26, when mean F o in shorebird presence and absence plots was highly significantly different ( Table 3).
The two subsequent sampling events (days 45 and 64) showed decreasing F o with progressively smaller differences between presence and absence plots. Mean F o in shorebird absence plots was still higher on day 45 but was not significantly different (Bonferroni corrected 0.01 significance level), and by day 64, F o levels were very similar between treatments (Figure 2A). There was no significant difference in chlorophyll a content (µg g −1 ) in the top ∼2 mm of sediment between presence and absence plots on day 45 ( Figure 3B).

Sediment Erodibility
To evaluate the effect of shorebirds on erosion protection, erosion threshold (τ cr ) was measured on day 45. Significantly greater erosion threshold was found in shorebird absence plots than in presence plots ( Figure 3C, Table 3).

Sediment-Water Nutrient Fluxes
There was significantly greater net nitrate influx into the sediment when shorebirds were absent compared to when they were present ( Figure 4A) and a significantly greater net nitrite efflux from the sediment into the water column when shorebirds were present (Figre 4B).
There was no significant difference in net phosphate flux between shorebird presence and absence plots. However, under lit conditions mean values changed from an influx into the sediment to a small efflux into the water column (Figure 4D), which is considered ecologically significant.
There was no significant difference in net dissolved organic carbon (DOC) flux between shorebird presence and absence plots ( Figure 4E). However, in shorebird presence during light incubation, we found a large reversal in flux direction of DOC into the sediment rather than the water column ( Figure 4E).
No significant difference in ammonium flux between the sediment and water column was found ( Figure 4C).

Macrofauna Density
To evaluate the indirect effect of shorebirds on erosion protection, nutrient cycling and carbon sequestration via changes in macrofauna density, the numbers of macrofauna were counted (from the same cores used for the nutrient measurements). Macrofauna recorded on day 45 were mud snails P. ulvae, Baltic clams Macoma balthica, midge larvae (Chironomidae), ragworms Hediste diversicolor, Arctic barrel-bubble Retusa obtusa and common cockles Cerastoderma edule. Mean densities (m −2 ) in each treatment are shown in Table 4. Raw macrofauna counts revealed presence of a single specimen of C. edule and R. obtusa in only two and three plots, respectively. H. diversicolor counts were also sparse (see Table 4). On day 26 P. ulvae was visually noted on the mudflat surface for the first time during F o sampling. Mud snails can compensate for the loss of higher predators on intertidal mudflats (Hamilton et al., 2006;Cheverie et al., 2014). This species was subsequently present within the study area during all F o sampling events, noted throughout the study site in presence and absence plots.
The non-metric Multi-Dimensional Scaling (nMDS) plot ( Figure 3D) indicated that macrofauna communities between treatments were not significantly dissimilar; a large overlap between community composition is indicated, although the spread of data points is larger in shorebird presence demonstrating larger variability in community composition. ANOSIM confirmed there was no significant difference in community composition between shorebird presence and absence plots (R = −0.038, P = 0.623).

Bird Surveys
Over the study period, 10 shorebird species were recorded using the wider mudflat, with a total of 78,811 bird days ( Table 2). Of these, three were recorded in the presence plots; C. alpina (84 bird-days), T. totanus (35 bird-days) and P. squatarola (28 birddays). Camera data indicated that numbers of shorebirds using the study area were broadly consistent with those counted during surveys. Although the image quality (due to distance from the plots) made detection of individual birds difficult, flocks were noted using the plots, often as the tideline crossed them. Flocks were noted on camera footage in and around the plots between 23 February and 5 March (day before main sampling event). The experimental plots were laid out in an area of mudflat representing approximately 0.3% of the area visually surveyed. Peak C. alpina, P. squatarola, and T. totanus numbers within experimental plots comprised approximately 0.16%, 0.35%, and 0.8% (respectively) of peak numbers within the survey area, thus within the same order of magnitude as that expected based on the areas of plots and the overall mudflat area.

DISCUSSION
Excluding shorebirds caused significant changes in regulating and provisioning ecosystem functions, including mudflat erodibility, nutrient fluxes and carbon sequestration. Effects on MPB biofilm biomass and erodibility were, however, not as predicted in our hypotheses. We suggest that these effects were driven by shorebird bioturbation of surface sediments and MPB biofilms and possible direct grazing of MPB by C. alpina.

Effects on MPB and Erodibility
Hypothesis 1 was not rejected; our linear mixed-effects model showed a highly significant difference in F o between shorebird presence and absence, with no significant interaction between other factors (see Table 3). Significantly greater MPB F o values were found in shorebird absence plots on days 13 and 26. By day 45 the difference had become less significant, to the extent of being non-significant when Bonferroni correction was applied (0.01 level). Despite this, on day 45 the difference in F o remained visually notable in the field, which is reflected in Figure 2A. These differences between treatments occurred during a period of increased shorebird activity in the study area. Despite the decline in surveyed shorebird numbers on day 28, the 83 ha −1 shorebirds present at this point was notably greater than at the beginning or end of the experiment (when numbers were 30 and 28 ha −1 , respectively) (Figure 2A). The survey visit on day 28 is also considered to be an underestimate due to the flushing of a large proportion of the foraging shorebirds on the incoming tide by a marsh harrier Circus aeruginosus. Differences in F o between shorebird presence and absence on days 3 and 64 were non-significant and occurred when shorebird numbers were smaller, suggesting that the effects found may be dependent upon shorebird density. There was no significant dissimilarity in macrofauna community structure between shorebird presence and absence plots ( Figure 3D). The present study recorded a greater diversity of species at the study site than during previous large scale work at the site (Wood et al., 2015), albeit the majority of infaunal species were present sporadically and in very low numbers (see Table 4). This validates our macrofauna sampling effort, in that we had enough replicates to detect all species known to be present, despite likely patchiness in invertebrate distributions (Van Colen, 2018). These findings differ to suggestions that a top-down ecological cascade effect driven by shorebirds can increase biofilm biomass (Daborn et al., 1993), supporting instead more recent work (Hamilton et al., 2006;Cheverie et al., 2014). Our results provide strong indication that, through bioturbation and/or grazing (and/or a yet unknown pathway), shorebirds can have a significant reductive effect on the biomass of surface MPB biofilms. Thus, shorebirds can alter key ecosystem functions such as erosion protection and nutrient cycling via direct and/or indirect effects on MPB. The increase in MPB in the absence of shorebirds concurs with results reported by Hamilton et al. (2006), where the authors acknowledge that this finding is the opposite to that expected in the event of a trophic cascade. On day 45 bulk chlorophyll a content within the surface 2 mm of sediments showed the same directional response as surface biofilm biomass and was also not significantly different. Bioturbation and grazing by macrofauna can significantly affect surface MPB biomass and resuspension (Grant and Daborn, 1994;Hagerthey et al., 2002;Harris et al., 2015); but as macrofauna were not significantly different between our shorebird presence/absence plots, and motile macrofauna could access all plots, the changes in MPB biomass are highly unlikely to have been due to macrofauna. Physical effects of birds upon primary producers is evident within many freshwater and marine environments (Cadee, 1990;Mitchell and Perrow, 1998;Nacken and Reise, 2000) and physical mixing of intertidal mud has been shown to significantly reduce chlorophyll a, F o, and colloidal carbohydrate (Tolhurst et al., 2012). It follows that physical disturbance (bioturbation) by shorebirds in our study location, through foraging (including biofilm grazing in some species) and tracking (walking), can have a significant effect upon MPB biomass and related sediment properties. These results suggest that bioturbation by shorebirds can be a more significant driver of effects on MPB than trophic cascades. Further work is required to confirm the mechanisms by which shorebirds in this part of the world reduce MPB biomass.
Hypothesis 2 was not rejected, sediment critical erosion threshold (τ cr ) was significantly smaller when shorebirds were present than when they were absent (see Figure 3C). This pattern is most likely to have been driven by both direct bioturbation during walking and feeding of shorebirds on the mudflat surface and, because MPB commonly significantly increase mudflat erosion threshold (Hale et al., 2019;Hope et al., 2020), indirectly by grazing decreasing the biomass of MPB. The exact mechanistic pathway(s) and their magnitude require further investigation. The erosion shear stresses exerted on intertidal mudflats by combined waves and tides are very variable, but commonly in the 0-1 Nm −2 range and typically below 4 Nm −2 (Christie and Dyer, 1998;Whitehouse and Mitchener, 1998). Thus, the τ cr measurements suggest that erosion would occur frequently (i.e., during most tidal cycles) in the presence of shorebirds and much less frequently in the absence of shorebirds. Given the importance of sediment erodibility for many ecosystem functions Hope et al., 2020), including nutrient fluxes and erosion protection; the effect of shorebirds on erodibility demonstrates their importance as ecosystem engineers (Passarelli et al., 2014) and their significant role in ecosystem functioning.
Although F o is widely used as a proxy for MPB biomass, it is important to acknowledge that this relationship varies depending upon the physiological state and taxonomic composition of MPB due to vertical migration of MPB (Serodio et al., 2001(Serodio et al., , 2006Serodio, 2004;Du et al., 2018). By standardizing our time of sampling within the tidal exposure period, tidal migration rhythms influencing F o were accounted for between treatments. Though changes in the relationship between F o and Chl a over time may have occurred, we found significant differences in F o between treatments at each time of sampling. Our results show the same directional response of F o and Chl a to shorebird presence, suggesting an underlying relationship in this case. Actual Chl a concentration varies vertically within the sediment depending upon factors such as MPB migration, light intensity, water content and sediment compaction (Perkins et al., 2003;Tolhurst et al., 2003;Jesus et al., 2006a;Maggi et al., 2013) and also shows temporal changes. We did not design our sampling regime to specifically focus on the F o to Chl a relationship, which requires a higher level of sampling granularity.

Effects on Nutrient Fluxes
Hypothesis 3 was not rejected; statistically significant differences in the fluxes of nitrate, nitrite and dissolved organic carbon (DOC), were found between presence and absence treatments. Orders of magnitude changes in the scale of some fluxes were observed (nitrate ∼100x, nitrite ∼10x and DOC ∼2000x).
Despite not being formally significant, the reversal of phosphate flux into/out of the sediment is considered to be ecologically significant. These results suggest that shorebirds significantly alter ecosystem functioning associated with nutrient cycling  -Beat et al., 2013;Mathot et al., 2018;Hope et al., 2020) and carbon storage (Maher and Eyre, 2010). Differences in the surface active MPB biomass (F o ) can explain the nutrient flux alterations by shorebirds. Photosynthesis and nutrient assimilation by MPB significantly affects nutrient flux rates, including nitrate (Dong et al., 2000) and phosphate (Sundback et al., 1991). Further, the EPS matrix within MPB biofilms provides additional organic matter to support heterotrophic bacteria, which reduce nitrite to nitrous oxide (Dong et al., 2002). We found evidence to suggest that the presence of shorebirds can significantly reduce nitrate uptake into intertidal sediments ( Figure 3A). The reduction of active surface MPB biofilms by shorebirds is a likely mechanism that may reduce nitrate and phosphate uptake, nitrification, coupled nitrification-denitrification, and through the reduction of extracellular organic carbon, reduce bacterial degradation rates (Thornton et al., 2007). Our findings suggest that shorebird effects on MPB can limit the drawdown of nitrate, nitrite and phosphate into sediments in an already nitrate rich estuary (Thornton et al., 2007). The observed alterations of nutrient fluxes suggest that shorebirds play a significant role in estuarine nutrient pathways, effectively controlling and engineering nutrient fluxes between the sediment and water column (Passarelli et al., 2014. Bioturbation by macrofauna is known to significantly affect nitrate and ammonia fluxes at the study site and elsewhere, through sediment reworking, ventilation and burrowing (Nizzoli et al., 2007). We suggest that bioturbation by shorebirds (Mathot et al., 2018) is likely to have contributed to the significant effects found here.
While the measured nutrients were typically characterized by a reduction in fluxes into the sediment from shorebird presence, DOC flux into sediment from shorebird presence increased significantly in lit conditions. It is possible that through the observed reduction of MPB biomass by shorebirds, competition for nutrients may have been reduced, allowing bacteria to proliferate and increase assimilation of DOC and ammonium (Amin et al., 2012). Migratory birds can also introduce bacteria to communities (Steiniger, 1969) via fecal droppings (Muller, 1965) and external tissues (Muza et al., 2000), potentially further increasing these process rates. These results indicate that changes in shorebird abundance could affect wider ecosystem functioning such as carbon sequestration and coastal biogeochemistry more broadly (Nedwell et al., 2016;Hope et al., 2020).

Secondary Effects
Use of the mid and upper shore at low tide by C. alpina, despite often being a 'tide follower' (Granadeiro et al., 2006), may have been driven by the visual cues of MPB communities on the mudflats, either as a cue for the presence of invertebrate prey or to feed upon MPB directly (Hamilton et al., 2003;Drouet et al., 2015;Jimenez et al., 2015). C. alpina is an opportunistic feeder with a broad diet (Dierschke et al., 1999) using visual and tactile foraging cues , and possibly exploited areas with high diatom biomass to maximize the breadth of feeding opportunity.
Avian guano (in particular shorebird droppings) is a potentially important source of nutrients in coastal areas (Schrama et al., 2013). It has been suggested that C. alpina droppings increase growth rate and biomass of the diatom species Entomoneis paludosa through increases in nitrogen and phosphorous input to the sediment (Jauffrais et al., 2015). However, the Colne estuary has very high nutrient loads (McMellor and Underwood, 2014;Nedwell et al., 2016) and MPB biomass was smaller, rather than larger in shorebird presence, suggesting that nutrient enrichment of biofilms by guano is not a major mechanism in this case. These findings reflect the complexity of the real-world scenario compared to laboratory studies (Jauffrais et al., 2015); in the present study shorebirds reduced MPB biomass on the upper shore. This indicates that the effects of bioturbation and/or grazing by shorebirds, which lead to alterations in ecosystem functioning, significantly outweigh the effects of nutrient input via guano in our study site.
Shorebirds significantly affect ecosystem functions (nutrient flux and erodibility), at least within the upper shore, in a temperate climate during late winter. However, these effects are likely to vary temporally and spatially (Underwood and Paterson, 1993;Gerwing et al., 2015) depending as they do upon the abundances and functioning of other organisms present (Underwood, 1994;Norazlimi and Ramli, 2014). For example, we found that shorebird effects were temporary and seasonal, restricted to an approximately one month period when shorebird density peaked at the study site (Figure 2A). This suggests that the observed phenomenon is seasonally and density dependent, reliant on sufficient density of shorebirds (which are present in larger densities during winter) to cause effects on ecosystem functioning. Similarly, compensatory grazing by the mud snail Peringia ulvae may have limited the temporal effect of shorebirds on MPB during this study, effectively resetting the state of the system as bird density declined (Hamilton et al., 2006;Cheverie et al., 2014). The collapse of the shorebird effect on F o was concomitant with the emergence of large numbers of P. ulvae. This MPB grazer was first noticed on the mudflat surface on day 26, was noted spread across the mudflat within all plots (Table 4), and can rapidly reduce the abundance and thickness of biofilms (Sahan et al., 2007). Subsequently the difference in F o between treatments steadily decreased, eventually becoming non-significant. On day 45, no significant difference between macrofaunal communities was evident. It is our interpretation that the snails had a homogenizing effect on biofilm distribution. Once the snails emerged and while birds remained, the effects of the birds became weaker. Once the birds left, continued grazing by the snails removed the residual bird effects (compensatory effect). Despite our restriction to observational evidence regarding the temporal change in numbers of P. ulvae, it is known that mudsnails can mask effects on MPB (Hamilton et al., 2006;Cheverie et al., 2014) and it is plausible that this occurred here, reducing the detectability of ecosystem function effect pathways. Here we highlight that shorebirds play a key community role in the regulation and control of ecosystem function, through inter and intraguild interactions with macrofauna and MPB with which they are intrinsically linked (Kuwae et al., 2012;Cheverie et al., 2014).
We found no evidence to suggest that macrofauna community structure differed between shorebird presence and absence, however, such effects have been detected in Canada in exclusion experiments on semipalmated sandpiper C. pusilla, where reductions in C. volutator densities were found (Hamilton et al., 2006;Cheverie et al., 2014). The differences between these studies may be due to geographic or shorebird species differences, or due to the fact that C. volutator is not present at our study site.
We also emphasize that differences in MPB surface biomass between treatments eventually became non-significant, despite shorebird exclosures remaining in-situ. We conclude therefore that shorebirds, rather than experimental artifacts, drove the measured MPB biomass changes and subsequent effects on ecosystem functions.

CONCLUSION
Here we have identified previously unknown effects of shorebirds on ecosystem functioning. Although limitations are acknowledged regarding the link between F o measurements and actual Chl a content, the end effect of shorebird presence on erodibility and nutrient fluxes was found to be significant, and a large amount of existing literature indicates that MPB are highly likely to drive this effect. The removal of shorebirds significantly increased surface biofilm F o and sediment erosion threshold. Shorebird absence was also found to affect nutrient cycling regimes and carbon sequestration on the mudflat; differences in biofilm biomass led to significant alterations in the flux of nutrients under lit conditions, including nitrate, nitrite and phosphate, all of which showed an increased flux into the sediment in the absence of shorebirds. The uptake of DOC in the light into the sediment was significantly greater in the presence of shorebirds.
The mechanism by which shorebirds reduced biofilm biomass was not experimentally tested, although the literature provides a number of possible drivers including physical disturbance (bioturbation) through tracking (walking) and foraging. Considering the presence of large numbers of C. alpina, which has been shown to consume MPB, it is plausible that direct consumption of biofilm may have contributed, but this is not confirmed. The lack of significant differences in macrofauna densities between treatments suggests that altered numbers of these invertebrates were not driving a change in bioturbation or grazing on the biofilms, and thus were not a significant driver of the measured effects.
The finite period of effects and community interactions between shorebirds, macrofauna and MPB reduce the clarity of the situation regarding consequences of declining shorebird species on coastal ecosystem functions. The work presented here indicates a potential shorebird density-dependent effect, resulting in stronger impacts on ecosystem function by birds during winter that may be 'reset' by other organisms or reduced bird densities in spring and summer. This reflects the complexity of intertidal mudflat ecosystem functions Hale et al., 2019), but is a step forward in disentangling the many factors influencing them. This research indicates that shorebirds play a significant role in the ecosystem functions provided by intertidal mudflats, including erosion protection, nutrient cycling and carbon sequestration. However, further research is required, involving longer-term, larger-scale experiments, to better understand the mechanisms behind ecosystem function regulation by shorebirds.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
JB conducted the initial experimental design and background research, performed primary data analysis and interpretation, was the primary manuscript author and led field and lab work. GU was a Ph.D. supervisor and provided quality control of the experimental design, flux core methodology, data analysis and interpretation. AP performed the laboratory data analysis of contact cores (Chl a concentrations) and some flux core samples (dissolved organic carbon) and contributed to the manuscript. RD was a Ph.D. supervisor and provided quality control of data interpretation regarding ecological community effects. TT was the Ph.D. project originator, primary supervisor and provided quality control of the experimental design, sediment properties measurements methodologies, data analysis and interpretation. GU, RD, and TT were manuscript reviewers and contributors. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the Natural Environment Research Council (grant number NE/L002582/1) and undertaken in partnership between the University of East Anglia and University of Essex within the EnvEast Doctoral Training Partnership.
Matt Cole, David Booty, Barbara Booty, Ken Hudgell, and David Smith for their assistance in the field and the lab. We thank Essex Wildlife Trust for their permission and cooperation. We thank Charlie Williams at Natural England for his cooperation and permission regarding the SSSI. We thank DH for her invaluable review and input during the publication process. Additional support was also provided by the University of East Anglia and EnvEast during fieldwork and writing.