Using laboratory experiments to characterize ecosystem-functioning under various natural biomasses, we previously demonstrated that stimulatory ragworm bioturbation effect on sediment oxygen uptake varies proportionally with population biomass in sandy and muddy sediments of tidal flat in the lower region of the estuary (Fang et al., 2021). Because macrobenthos activity can vary with salinity (Verdelhos et al., 2015; but see Murray et al., 2017), a third experiment (presented here) with ragworms collected in the upper region of the Western Scheldt (Supplementary Material 1) was undertaken following the same methods reported in Fang et al. (2021). Briefly, total sediment oxygen uptake was measured in 11 plexiglass microcosms (height 12.2 cm, inner diameter 3.6 cm) filled with 1 mm sieved sediment to a depth of 6 cm, topped with seawater (salinity = 19.4) and seeded with variable densities of differently sized individuals (Supplementary Material 2). Sediment oxygen uptake was calculated from the decline in dissolved oxygen concentrations over time; these were measured continuously in the airtight closed microcosms using calibrated oxygen optodes (Pyroscience OXROB10). The ragworm contribution to sediment ecosystem functioning was calculated for each microcosm by subtracting the oxygen uptake in one microcosm without ragworms. Analysis of covariance (ANCOVA) was used to compare the ragworm effect on ecosystem functioning between sediment types while controlling for biomass. Ragworm population biomass was normalized using log transformation in all experiments to meet assumptions of the homogeneity of regression slopes and the homogeneity of variances.

Four hydrodynamic variables, namely: (1) yearly averages of maximal tidal current velocity (m⋅s^–1); (2) salinity; (3) diel variation in salinity; and (4) inundation time (% of submersion during a full tidal cycle) – all known determinants of ragworm distribution (Ysebaert et al., 2003) were predicted across a 20 m triangular spatial grid overlaying detailed bathymetric maps of the Western Scheldt for 1955 and 2010 using the 2Dh TELEMAC hydrodynamic model (Smolders et al., 2013). Subsequently these hydrodynamic variables were extracted from the 2010 hydrodynamic scenario for 3051 locations sampled between 2006 and 2011 in the Western Scheldt (BIS dataset, NIOZ-Yerseke Monitor Taskforce, Figure 1C). Ragworm biomass distribution in the estuary was then predicted based on the observed biomass in the sampled locations and using the four physical variables extracted by the 2Dh TELEMAC hydrodynamic model as explanatory variables. To this end we used quantile regression, a preferred method for species distribution modeling when only a subset of possibly limiting governing variables are available (Cade and Noon, 2003). We predicted ragworm biomass using the upper (0.95) quantile to describe change in niche potential in response to environmental drivers, after Vaz et al. (2008). The distribution model was fitted up to second degree interaction terms between explanatory variables and subsequently simplified using a stepwise bi-directional elimination procedure (Venables and Ripley, 2002; Supplementary Material 3). The explanatory variable submersion time, expressed as a fraction between 0 and 1, was linearized via arcsin transformation. Similarly to Cozzoli et al. (2014), to validate our forecast for the 95th quantile of the H. diversicolor biomass distribution, the whole dataset was iteratively sampled with replacement (9999 iterations). Due to sampling with replacement, some observations are repeated and others remain unpicked. The model was fitted on the sampled observations (training dataset) and used to predict the unpicked ones (validation dataset). The distribution of the coefficients estimated through the iterations of the pseudoreplication procedure was used to estimate the significance and confidence intervals of the estimates of the distribution model. The predicted values were discretized in 10 logarithmic classes, for which the corresponding 95th sample quantile of the validation data was calculated. To finally asses the validity of the model, observed and predicted quantiles were plotted against each other and checked for linear correlation (Supplementary Material 4).

The relationship between ecosystem functioning and ragworm biomass (Stimulated oxygen uptake = aBiomass^b, where a is the quantified coefficient and b the allometric scaling exponent; see section “Results”) was used to convert the ragworm population effect on sediment oxygen uptake and constrained to the inundation time (i.e., the period of the tidal phase when bioturbation can stimulate sediment oxygen uptake from the water column). The predicted biomass distribution and ragworm-stimulated oxygen uptake maps were generated for 1955 and 2010 to evaluate the implications of deepening, dredging and disposal practices that increasingly occurred since 1955 (De Vriend et al., 2011). The pseudoreplication test previously described and used to validate the biomass distribution model requires a large training and validation dataset. Consequently, in this study we were able to validate the biomass distribution model only. A much greater number of field observations should be collected to properly validate the predicted ecosystem functioning scenario.

The biomass-dependent stimulation of sediment oxygen uptake by ragworms did not vary between sediment types (ANCOVA: biomass effect F₁,₂₄, p < 0.001; sediment type effect F₂,₂₄, p = 0.24; biomass × sediment type interaction F₂,₂₄, p = 0.40), despite differences in granulometry, organic matter content or salinity (Supplementary Material 1) and a combined stimulatory effect per gram biomass increased sediment metabolism to the power 0.71 (Figure 2) across all environmental contexts. This uniformity likely results from the fact that burrow ventilation by ragworms is restricted to the lumen and a few mm in the burrow wall (Nielsen et al., 2004) and suggests that the ragworm bioturbation effect on sediment ecosystem functioning is proportional to its energetic requirements corroborating the allometric scaling exponent of ∼0.75 for the relation between body size and metabolism that is evidenced across taxonomic groups (Brown et al., 2004). Cozzoli et al. (2018) further support this theory by demonstrating that sediment resuspension can be ascribed to the overall metabolic rate of bioturbating populations.

Maximum ragworm biomass occurred in the most elevated intertidal areas (Figures 3A,B) corroborating this species’ preference for low-energy environments, i.e., low current velocities and low inundation time (Ysebaert et al., 2002). These conditions are usually associated to fine sediments, i.e., mud or silt, and high organic matter sediment load (Cozzoli et al., 2013). In contrast, ragworm-mediated sediment oxygen uptake peaked at intermediate surface elevations (Figures 3C,D) due to the opposing patterns in inundation time and ragworm biomass along the intertidal elevation gradient, that both positively relate to stimulated ecosystem functioning.

Anthropogenically induced changes in fluvial geomorphology between 1955 and 2010 (Figure 3E) increased seawater penetration into the estuary, modified local daily salinity fluctuations, and altered the inundation time and current velocities in the intertidal region (Figures 4A–D). Overall, channel deepening and erosion of lower intertidal areas led to the permanent inundation of 34 km², accounting for the direct loss of 216 tons of potential intertidal ragworm biomass, or 22% of the potential biomass achievable in 1955 (Figure 3F). This was, in part, offset by creation of 20 km² of new intertidal habitat (i.e., previously fully submersed habitat) in close proximity to the shoals (Figure 3E) and potentially inhabited by 85 tons, or 13% of ragworm biomass present in 2010. The overall decrease in intertidal habitat surface was stronger in the marine and intermediate sectors of the estuary (respectively, shrinking of 19 and 17% of the intertidal surface present in 1955) than in the riverine one (shrinking of 9%) (Figure 3E).

The 73 km² of prevailing intertidal habitat underwent a distinct increase in average maximum current velocity, especially in the most elevated intertidal habitats (1955; 0.16 ± 0.10 SD m⋅s^–1 vs. 2010; 0.31 ± 0.11 SD m⋅s^–1; Figure 4 and Supplementary Material 5), reducing their suitability for ragworms. Similarly, newly developed intertidal habitats in 2010 had a higher maximal current velocity as compared to the lost habitat (gained habitat; 0.43 ± 0.12 SD m⋅s^–1 vs. lost habitat; 0.11 ± 0.12 SD m⋅s^–1), and experience a longer submersion time (+11%) in comparison to the lost habitats. The change in hydrodynamics was more drastic for the marine and intermediate portion of the estuary (respectively, average increases of 74 and 84% in maximum current velocity between 1955 and 2010) than for the riverine sector (increase of 25% only) (Figure 4B and Supplementary Material 5). Differently, the riverine sector of the estuary experienced the strongest increase in daily salinity variation, with an increase of 73% between 1955 and 2010 in the intertidal area (Figure 4D and Supplementary Material 5).

Collectively, direct habitat loss and change in habitat suitability over time reduced ragworm population biomass by 32% between 1955 and 2010, causing a predicted reduction in ecosystem functioning by 28% across the estuary (1955; 204 Kmol⋅day^–1 vs. 2010; 147 Kmol⋅day^–1)(Figure 3F). The joint effect of larger intertidal habitat loss and stronger increase in tidal current velocity made so that the strongest reduction in ragworm population biomass was predicted in the marine and central sectors of the estuary (respectively decrease in potential ragworm biomass of 43 and 38% from 1955 to 2010) than in the riverine one (decrease in potential ragworm biomass of 20%) (Figure 4E and Supplementary Material 5). This implies a predicted reduction in ecosystem functioning of 36 and 34% in the marine and central sectors of the estuary, respectively. The decrease in ecosystem functioning in the riverine sector was of 16% only (Figure 4F and Supplementary Material 5).

In addition to basin-wide estimates, our approach hindcasted changes in biomass at a smaller scale, with the implications for ecosystem functioning primarily depending on the local interplay of biomass and inundation time. For example, geomorphological changes generally reduced habitat suitability for ragworms on the river banks whereas implications for biomass and ecosystem functioning were more variable on the shoals with, e.g., predicted gains in habitat suitability toward the brackish part of the estuary where reductions in inundation time co-occurred with lowered current velocities (Figures 4E,F).