Hindcasting Ecosystem Functioning Change in an Anthropogenized Estuary: Implications for an Era of Global Change

Understanding how altered hydrodynamics related to climate change and anthropogenic modifications affect ecosystem integrity of shallow coastal soft-sediment environments requires a sound integration of how species populations influence ecosystem functioning across heterogeneous spatial scales. Here, we hindcasted how intertidal habitat loss and altered hydrodynamic regimes between 1955 and 2010 associated with geomorphological change to accommodate expansion in anthropogenic activities in the Western Scheldt altered spatial patterns and basin-wide estimates of ecosystem functioning. To this end we combined an empirically derived metabolic model for the effect of the common ragworm Hediste diversicolor on sediment biogeochemistry (measured as sediment oxygen uptake) with a hydrodynamic and population biomass distribution model. Our integrative modeling approach predicted an overall decrease by 304 tons in ragworm biomass between 1955 and 2010, accounting for a reduction by 28% in stimulated sediment oxygen uptake at the landscape scale. Local gains or losses in habitat suitability and ecosystem functioning were primarily driven by changes in maximal current velocities and inundation regimes resulting from deepening, dredging and disposal practices. By looking into the past, we have demonstrated how hydro- and morphodynamic changes affect soft-sediment ecology and highlight the applicability of the integrative framework to upscale anticipated population effects on ecosystem functioning.

Understanding how altered hydrodynamics related to climate change and anthropogenic modifications affect ecosystem integrity of shallow coastal softsediment environments requires a sound integration of how species populations influence ecosystem functioning across heterogeneous spatial scales. Here, we hindcasted how intertidal habitat loss and altered hydrodynamic regimes between 1955 and 2010 associated with geomorphological change to accommodate expansion in anthropogenic activities in the Western Scheldt altered spatial patterns and basin-wide estimates of ecosystem functioning. To this end we combined an empirically derived metabolic model for the effect of the common ragworm Hediste diversicolor on sediment biogeochemistry (measured as sediment oxygen uptake) with a hydrodynamic and population biomass distribution model. Our integrative modeling approach predicted an overall decrease by 304 tons in ragworm biomass between 1955 and 2010, accounting for a reduction by 28% in stimulated sediment oxygen uptake at the landscape scale. Local gains or losses in habitat suitability and ecosystem functioning were primarily driven by changes in maximal current velocities and inundation regimes resulting from deepening, dredging and disposal practices. By looking into the past, we have demonstrated how hydro-and morphodynamic changes affect soft-sediment ecology and highlight the applicability of the integrative framework to upscale anticipated population effects on ecosystem functioning.

INTRODUCTION
Coastal ecosystems are at the frontline of environmental change. Expanding human activities, accelerating sea level rise, and changing wind climates modify hydrodynamic patterns that define the delicate balance in physical exchange processes that underpin the ecology of these shallow ecosystems (Birchenough et al., 2015;Khojasteh et al., 2021). Consequently, worldwide loss in FIGURE 1 | Study area: the Western Scheldt estuary (Netherlands) (A); model organism: Hediste diversicolor, common ragworm (B); sampling points: the colored dots show the sampling points used for oxygen uptake measurements (C). The smaller black dots mark the location of the 3051 benthic samples collected between 2006 and 2010 and used to fit the H. diversicolor potential biomass distribution model. Coordinates are in Amersfoort -EPSG:4289. ecosystem integrity during the Anthropocene has reduced some of the critical services provided by the biodiversity of these highly productive ecosystems to society (Barbier et al., 2011 and references therein). According to the definition of Reiss et al. (2009), the concept of "ecosystem functioning" describes the combined effects of the natural processes that sustain an ecosystem and it is determined by the interplay of abiotic (physical and chemical) and/or biotic factors. A key ecosystem function of coastal ecosystems is the efficient processing of nutrient inputs from terrestrial sources, which could otherwise cause eutrophication and disrupt marine food webs. Soft-sediment dwelling animals (further referred to as macrobenthos, typically including crustaceans, bivalves, and bristle worms) play pivotal roles in the biogeochemical cycling of estuaries and coastal ecosystems by channeling nutrients and carbon assimilated into their biomass to higher trophic levels such as wading birds and fish. Additionally, these burrowing animals alter biogeochemistry at the sediment-water interface via respiration, excretion, and indirectly through their sediment reworking and ventilation activities (i.e., bioturbation; Kristensen et al., 2012) that importantly stimulate the microbial mineralization of sedimentary organic matter (Snelgrove et al., 2018).
Multiple anthropogenic pressures increasingly affect the water column and sediment processes involved in macrobenthosmediated biogeochemical cycling (Knights et al., 2015). Key examples are dredging and disposal practices that induce bathymetric, hydrodynamic, and sedimentological changes of the seabed that affect the distribution and bioturbation effects of macrobenthos (e.g., Mermillod-Blondin and Rosenberg, 2006;Cozzoli et al., 2017). Accordingly, quantifying the implications of anthropogenic change to shallow coastal and estuarine sediments for the delivery of ecosystem services, requires a sound integration of the contributions of macrobenthos to biogeochemical cycling, and sediment ecosystem functioning in general. Although the functional role of macrobenthos can exhibit a strong context-dependency across environmental gradients (e.g., Cassidy et al., 2020), logistical constraints typically limit replication of empirical studies across the landscape that can impede spatial upscaling of functioning at the ecosystem level and beyond. However, when the functional role of populations across spatial gradients is well understood, a promising way to estimate larger-scale functional implications of changing environments is the integration of species ecosystem functioning models with predictive population distribution models across environmental gradients, even when the underlying mechanisms are complex and difficult to unravel (e.g., Ellis et al., 2006). Such integrated modeling framework can extrapolate changes on population level to ecosystem functioning across heterogeneous spatial scales, or can make predictions for different time periods. Hydrodynamic variables are amongst the most important determinants for macrobenthos distributions (Ysebaert et al., 2002), but they are rarely measured with full spatial coverage therefore limiting landscape applications, let alone that the temporal changes that occurred during the Anthropocene has been recorded. Coupled hydro-and morphodynamic models fill this gap as they describe water motion, sediment transport and bed-level changes by numerically solving a coupled set of mathematical equations (Smolders et al., 2013).
To demonstrate this approach we integrated a metabolically supported species-ecosystem functioning model with species distribution predictions to assess the spatial and temporal variation in ecosystem functioning in the Western Scheldt (Netherlands) as an exemplar study ( Figure 1A). The Western Scheldt is a coastal plain estuary of high ecological and economic importance protected under the EU Water Framework and the Floods Directive. The estuary has underwent severe geomorphological changes since the 1950s (De Vriend et al., 2011) to accommodate safe navigation of larger ships to the port of Antwerp, the second-largest marine harbor in Europe. Ecosystem functioning was measured using sediment oxygen uptake, which is a good proxy for total sediment metabolism (Canfield et al., 1993) and includes the direct and indirect contributions of macrobenthos to biogeochemical cycling (Glud, 2008). We specifically considered the contribution to sediment metabolism of the benthic ragworm Hediste diversicolor ( Figure 1B); a species widely distributed across the coasts of the Americas and Europe including the intertidal river banks and shoals of the Western Scheldt (Ysebaert et al., 2003) and the bioturbation activities of which importantly stimulate sediment biogeochemical cycling (e.g., Kristensen, 2001). We first established an ecosystem functioning model by empirically quantifying the contribution of ragworm population biomass to stimulated sediment oxygen uptake across a diversity of sediment types. Then we extrapolated stimulated ecosystem functioning across the spatially heterogeneous landscape by integrating the ragworm ecosystem functioning model with a ragworm biomass distribution model based on hydrodynamic predictor variables. Lastly, we estimated spatio-temporal variability in ecosystem functioning between 1955 and 2010 as a result of anthropogenically driven loss and/or change in habitat suitability related to hydrodynamical modifications.

Ragworm -Ecosystem Functioning Model
Using laboratory experiments to characterize ecosystemfunctioning 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.

Upscaling Ragworm -Ecosystem Functioning Effects
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.

Ragworm-Ecosystem Functioning Model
The biomass-dependent stimulation of sediment oxygen uptake by ragworms did not vary between sediment types (ANCOVA: biomass effect F 1 , 24 , p < 0.001; sediment type effect F 2 , 24 , p = 0.24; biomass × sediment type interaction F 2 , 24 , 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.

Habitat Change and Ecosystem Functioning
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 . 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 2 , 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 2 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 2 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).

CONCLUSION
By combining quantitative estimates of functioning under varying environmental contexts, spatio-temporal data describing hydrodynamic conditions, and modeling tools, we have shown how small-scale studies can be upscaled from a population perspective to the larger spatial (landscape) scales necessary to underpin effective ecosystem-based management in the current era of global change. The integrative modeling framework can be applied to predict functional implications of hydrodynamic change resulting from, e.g., alternative dredge disposal strategies, sea level rise, or river flooding regimes. By looking in the past in the specific case study presented here, we demonstrated how particularly the ultimate balance between maximal current velocities and inundation time in the intertidal habitat will define faunal-mediated functioning in future estuaries and soft-sediment shorelines. Finally, the allometric scaling of ecosystem functioning with population biomass corroborates the metabolic theory of ecology, supporting the possibly wider application of the integrated framework to assess ecosystem functioning change.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
CV and XF conceived and designed the study, analyzed the experimental data, and wrote the first draft of the manuscript. XF ran the experiment. SS developed the hydrodynamic model. XF and FC modeled the species and ecosystem functioning distributions. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was funded by a MARES Ph.D. grant to XF (2012-1720/001-001-EMJD). The research leading to results presented in this publication was carried out with infrastructure funded by EMBRC Belgium -FWO project GOH3817N. Extra funding for this project was obtained from the Special Research Fund (BOF) from Ghent University through GOA-project 01G02617.