Linking marine benthic biodiversity and ecosystem functions related to carbon cycling in a continental mud depocenter

The importance of carbon storage in continental seafloor sediments is increasingly recognized, yet the role of benthic macrofaunal biodiversity in the regulation of these processes remains poorly understood. Benthic macrofauna contributes to organic carbon cycling through respiration and secondary production, while sediment reworking (bioturbation) and ventilation (bioirrigation) by infauna promote the redistribution and remineralization of organic matter in sediments. Here, we investigated how benthic community structure, functional traits, and the relationship between biodiversity and calculated potential ecosystem functioning related to carbon cycling vary along environmental gradients in muddy sediments of the southeastern North Sea. Based on 171 macrofaunal taxa collected from 50 stations, hierarchical clustering based on abundance data categorized the macrofaunal communities into four distinct structural groups. Distance-based redundancy analysis (dbRDA) revealed that community functional composition was primarily structured by bottom hydrodynamics (salinity and bottom shear stress) and sediment characteristics (sorting, skewness, kurtosis, and D50). Biomass-weighted functional trait vectors, mapped onto a non-metric multidimensional scaling (nMDS) space, showed a clear structural shift with communities in the deep clusters being dominated by mobile biodiffusors and subsurface filter feeders, whereas shallow clusters were characterized by less mobile, surface-modifying bivalves and polychaetes. These contrasting patterns led to pronounced differences in calculated potential ecosystem functioning: bioturbation and bioirrigation potentials were significantly higher in deeper communities, whereas the calculated community secondary production and respiration were higher in shallow communities.

Across all stations, community secondary production and respiration increased with taxonomic and functional richness, while bioturbation and bioirrigation potentials were negatively related to biodiversity and community evenness. Interaction regressions and simple slope analyses revealed that environmental gradients significantly modulated the biodiversity–ecosystem functioning (BEF) relationships. Our findings demonstrate that environmental gradients fundamentally shape both benthic community structure and the link between biodiversity and essential ecosystem functions related to carbon dynamics in fine-grained sediment and organic matter in shelf sea systems.

1 Introduction

The current global decline in biodiversity and species abundances (Dirzo et al., 2014) is expected to cause a loss of essential ecosystem processes and functions (Oliver et al., 2015). Therefore, understanding how ecosystem functioning is linked to biodiversity is crucial for explaining ongoing ecological changes and for the development of sound strategies for management and adaptation. The connection between biodiversity and ecosystem functioning has gained significant attention in ecology and environmental sciences (Loreau et al., 2001; Tilman et al., 2014). Over the last 2 decades, hypotheses regarding biodiversity–ecosystem function (BEF) relationships have been increasingly addressed in marine ecosystems. Fundamental ecosystem functions often examined in this context are community productivity and respiration, typically measured in terms of biomass and oxygen consumption, respectively.

The link between productivity, biomass, and biodiversity is often debated. In coral reef systems, both fish community biomass and production increase with the introduction of new species (Mora et al., 2011; Duffy et al., 2016, 2017). Similarly, benthic meiofaunal biodiversity is positively related to biomass and to secondary production (Danovaro et al., 2008). However, weak or even negative relationships between richness and biomass have also been found. For example, abundance, biomass, and productivity were independent of variations in species richness in demersal fish (Greenstreet et al., 2012) and benthic invertebrate (Heip et al., 1992) communities of the North Sea, indicating that species richness and productivity are not necessarily linked to each other. Maureaud et al. (2019) suggested that community evenness rather than species richness controls ecosystem functions based on the observation that fish biomass in European seas tended to be higher in assemblages dominated by a few generalist species capable of exploiting both benthic and pelagic resources. Similarly, the relationship of community respiration and species richness varies across systems. For instance, Bolam et al. (2002) found that macrofaunal species richness and biomass were significantly correlated with oxygen consumption, yet the strength of this link appeared to depend rather on single species, such as predatory species with high energetic demands (Arndt and Schiedek, 1997), than on species richness. Taken together, these studies demonstrate that the relationship between biodiversity, productivity, and respiration is not universal, but rather dependent on the composition and structure of communities and on environmental conditions.

The relationship between biodiversity and ecosystem functions is not static but may vary in strength and direction along environmental gradients. In marine benthic systems, gradients in depth, sediment properties, and hydrodynamic forcing can alter community composition and functional trait distributions (van der Wal et al., 2017; van Son et al., 2013), thus regulating how biodiversity translates into ecosystem functioning. For example, the relationship between benthic macrofauna diversity and functions related to carbon fixation and mineralization depended on the levels of nutrients and turbidity in coastal ecosystems of New Zealand (Gammal et al., 2023). Such dependency can also result from environmental gradients shaping single-species trait expression and macrofaunal behavior. For instance, burrowing species behave differently in response to variations in sediment characteristics (Needham et al., 2011), with benthic macrofauna exhibiting faster burrowing and more intense bioturbation in soft, low-bulk-density sediments (Wiesebron et al., 2021). Although numerous studies have examined patterns of biodiversity along environmental gradients, the extent to which environmental gradients shape BEF relationships remains unclear. Addressing this gap, we examine how BEF relationships in marine benthic macrofauna vary across sedimentary and hydrodynamic gradients.

In marine benthic habitats, BEF relationships are largely mediated by macrofaunal communities inhabiting the seafloor, which play a pivotal role in ecosystem functioning through their influence on sedimentary processes (Lohrer et al., 2010; Braeckman et al., 2014). Benthic invertebrates strongly affect sedimentary biogeochemical processes through sediment reworking (bioturbation) and ventilation (bioirrigation), thus affecting oxygen consumption, carbon remineralization, element cycling, and bacterial activity in sediments (Aller, 1994; Wenzhöfer and Glud, 2004; Glud, 2008; Kristensen et al., 2012; Gilbertson et al., 2012; Laverock et al., 2014). However, the proportional contribution of the macrofauna to these processes depends decisively also on local sediment characteristics and hydrodynamic forcing (Bernard et al., 2019). Given the ecological importance of benthic macrofauna in driving sedimentary biogeochemical processes, it is essential to understand what determines the variability of these ecosystem functions. Biodiversity, and particularly the functional composition of benthic communities, is thought to play a key role in mediating these processes. Theoretical studies suggest that, in addition to taxonomic diversity, functional diversity (i.e., the combination of functional groups within a community) controls relationships between biodiversity and ecosystem functioning (Tilman et al., 1997; Hooper and Vitousek, 1997; Loreau, 1998). Cadotte et al. (2011) suggested that functional diversity explains ecosystem functioning even better than species richness. As the functional traits of organisms determine their resource utilization, ecosystems with higher functional diversity are expected to feature enhanced ecosystem functioning (Hooper et al., 2005). Moreover, biological trait structure reflects local-scale environmental conditions. This approach helps differentiate communities and thus provides better ecological insight than methods considering the taxonomic composition only (Bremner et al., 2003). For example, trait analysis has been applied in biomonitoring in macroinvertebrate communities in aquatic ecosystems for decades (e.g. Charvet et al., 1998), highlighting the added value of functional approaches when linking biodiversity to ecosystem functioning.

The benthic fauna of the southeastern North Sea has been studied extensively over the past decades with regard to community structure, species' abundances, biomass, and diversity (Künitzer et al., 1992; Kröncke et al., 2004; Reiss et al., 2010), and anthropogenic and climatic influences have been identified as drivers of long-term changes of the benthic system (Reid and Edwards, 2001; Kröncke et al., 2011; Shojaei et al., 2016; Meyer et al., 2018). However, community structure alone is often not effective in determining ecosystem changes (Loreau, 2010). Accordingly, functional ecology approaches have been applied to describe functional patterns and biogeochemical interactions. Solan et al. (2004) established a theoretical measure of bioturbation, while Queirós et al. (2013) developed a classification for macrofauna species depending on feeding modes and behavior to estimate the community level of bioturbation potential (BP_C). Building on this method, Wrede et al. (2017) identified key bioturbators for the German Bight and their impacts on nutrient fluxes. Wrede et al. (2018) proposed a similar method to estimate the bioirrigation potential (IP_C) of benthic communities based on burrow type, feeding type, and the depth of water pumping during bioirrigation (i.e. injection pocket depth), and applied this approach to model the oxygen and nutrient flux dynamics across the sediment–water interface (Wrede et al., 2019; Neumann et al., 2021). When applied to extensive datasets on benthic community composition, these measures can allow for exploring the relationship between structural diversity and important ecosystem functions.

We explored the relation between biodiversity and ecosystem functioning for the benthic macrofauna in an extensive area of muddy sediments in the southeastern North Sea. This area has been identified as a major depocenter for fine-grained material with considerable potential for benthic carbon sequestration (Müller et al., 2025; Wei et al., 2025). Recent geochemical investigations have revealed the distinctive role of mud depocenters on sandy shelves as potential hotspots for carbon burial and their potential relevance in the context of global climate change (Wei et al., 2025). Fine-grained sediments accumulate primarily off the mouths of major rivers (Müller et al., 2025), driven by specific hydrodynamic conditions that promote sediment trapping and retention (Chen et al., 2025). As a result, these areas exhibit rapid organic matter deposition and exceptionally high carbon burial efficiency, highlighting the biogeochemical significance of local depocenters (Hebbeln et al., 2003; Wei et al., 2025; Müller et al., 2025). While the geochemical importance of the muddy sediments as a carbon sink is well established, a substantial portion of research historically focused on physical and sedimentological processes. Growing attention has been directed toward understanding how benthic macrofauna community composition varies spatially within depocenters (e.g., Kauppi et al., 2017; Bonifacio et al., 2018; Lamarque et al., 2022, 2025). Despite these valuable insights from various marine systems, characterizing the specific spatial variations and the relative contributions of benthic macrofauna in regulating carbon cycling within localized shelf-sea depocenters requires further resolution. Therefore, our study on the relationship between biodiversity and ecosystem functioning specifically focusses on functions of the benthic macrofauna related to carbon cycling.

Benthic macrofauna contributes to carbon cycling in multiple, potentially antagonistic ways by incorporating organic carbon into secondary production, respiratory release of carbon as CO₂, and modifying sediment structure through bioturbation and bioirrigation. Through these activities, macrofauna redistribute organic matter in sediments and either promote or impede its remineralization, and thus affect the balance between carbon burial and release. Specifically, our objectives are: (1) to characterize the taxonomic and functional trait composition of the benthic macrofauna community across the muddy study region, (2) to identify environmental factors shaping the benthic macrofaunal community composition, and (3) to explore the relationships between biodiversity and ecosystem functions related to biogeochemical processes and their variations along environmental gradients.

2 Materials and methods

2.1 Study area and sampling

The study area was located in the muddy regions of the southeastern North Sea, a shallow shelf sea with water depths ranging from 20 to 60 m off the intertidal Wadden Sea (Bockelmann et al., 2018). Sea bottom temperatures typically range from ∼3 °C in winter to ∼17 °C in summer (Elliott and Clarke, 1991). The region is strongly influenced by continental freshwater input, resulting in coastal salinities below 30 and increasing offshore to 31–33 near the island of Helgoland (∼50 km offshore). Strong tidal and wind forcing permanently mixes waters in the areas shallower than 20 m (Otto et al., 1990; Huthnance, 1991), while mediating periodic and non-periodic stratification events in deeper zones (Burchard and Hetland, 2010; Chegini et al., 2020), and generally drives an anti-clockwise residual circulation along the coast (Kopte et al., 2022). A persistent frontal system along the 30 m depth contour separates coastal from offshore waters of the German Bight (Krause et al., 1986; Dippner, 1993). Sediments range from extensive muddy and sandy deposits to localized areas of coarse sand and glacial debris, particularly around Helgoland and off the island of Sylt (Diesing et al., 2006; Michaelis et al., 2019). A key geomorphological feature is the paleo-valley of the River Elbe, extending northwestward from the present-day river mouth and characterized by elevated mud content and enhanced organic matter input from river discharge (Figge, 1981).

Benthic macrofauna was sampled during expedition HE625 of the R/V Heincke in July 2023 at 50 stations spread evenly over the muddy sediment region in southeast of Helgoland at water depths of 7 to 51 m (Fig. 1).

Figure 1 Sampling stations clustered by the composition of the benthic community as quantified in terms of ash-free dry mass.

At each station, two replicate samples were taken with a van Veen grab with a sampling area of 0.1 m² and a minimum penetration depth of 10 cm. The first sample was sieved through a 1000 µm mesh to retain the macrofauna. The sample was fixed in a 4 % formalin–seawater solution buffered with sodium tetraborate and stored for at least 4 weeks before analysis. All macrofaunal organisms were counted and identified to the species level if possible. The taxonomy was matched against the World Register of Marine Species (WoRMS Editorial Board, 2021). For each taxon, the total wet weight (including shells) was measured with a precision of 0.001 g. From the second grab, a subsample of the upper 6 cm of sediment was taken using a coring tube with a diameter of 4.5 cm, and frozen at −20 °C for later analysis. After thawing, the grain-size distribution was determined for a sub-sample based on laser-diffraction granulometry using a CILAS 1180L particle size analyser. Another sub-sample of 40 g of sediment were dried, weighed, and incinerated at 500 °C for 5 h to estimate the organic content (%) as weight loss on combustion. The bottom-water temperature was taken as the sediment temperature measured in each sample. We acknowledge that utilizing only a single sample per station may introduce uncertainty in our measures of biodiversity and increase the probability of missing rare or patchily distributed taxa. Instead of increasing the level of replication at the single stations, we decided for a broad coverage of the entire study area with a single sample at each station in order to capture the more pronounced variations in the infauna communities along environmental gradients. We expected the gradients in the infauna communities across the sampling area to allow for a more robust analysis of the variations in the relationship between biodiversity and ecosystem functioning than local variations among replicate samples at single stations.

Spatial fields of maximum tidal bottom shear stress were obtained from the barotropic FESOM-C setup for the North Sea, with a focus on the German Bight. Model resolution ranges from 50 to 100 m in the study area. Tidal forcing at the open boundaries was prescribed using the TPXO9 atlas (Egbert and Erofeeva, 2002). FESOM-C is the coastal configuration of the global Finite-volumE Sea ice Ocean Model (Androsov et al., 2019) and has been extensively tested and validated for the North Sea in both idealized and realistic simulations (Fofonova et al., 2019, 2021; Kuznetsov et al., 2020; Sprong et al., 2020; Sidorenko et al., 2025). In addition, spatial data on bottom salinity were derived from the NEMO model for the Northwest European Shelf, which incorporates data assimilation and has a horizontal resolution of 1.5 km (Tonani et al., 2019). Seasonal bottom salinity was averaged over the multi-year period 2019–2022, and summer values were selected for subsequent analyses to match the seasonal salinity conditions with the timing of our macrofauna sampling, which also took place in summer.

2.2 Estimation of ecosystem functions

2.2.1 Calculation of benthic community secondary production and respiration

In a single case, a large individual of the ocean quahog Arctica islandica (195 g wet weight) was excluded from the dataset as an exceptional outlier, because the biomass-focused analyses would otherwise have been disproportionately dominated by this observation. Macrofauna biomass (in g wet weight) was first converted into energy units (Joules, J) using taxon-specific conversion factors provided by Brey et al. (2010). Annual secondary production (P) of the macrofaunal community was estimated based on three continuous variables: body mass, water depth, bottom-water temperature. The bottom-water temperature was taken as the sediment temperature measured in each sample (ranging from 15 to 18 °C; with calculations for each species applied to the corresponding temperature of its sampling station), and 17 categorical parameters, including five taxonomic groups, seven lifestyles, four environmental types, and the state of exploitation, using the Artificial Neural Network (ANN) model developed by Brey (2012), carried out via the BenthicPro R package (Andresen and Brey, 2018). In parallel, community respiration (R) was estimated using the respiration ANN model (Brey, 2010), based on an extensive database on respiration rates of aquatic invertebrates (approximately 23 000 measurements across more than 900 species from over 450 references). This model incorporates three continuous variables (body mass, temperature, and water depth), and 25 categorical parameters, including 19 taxonomic groups, three mobility modes, carnivory, visibility, and starvation state. Both community secondary production and respiration are converted into the unit of $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

2.2.2 Calculation of community bioturbation and bioirrigation potentials

The community level bioturbation potential (BP_C) was calculated following Solan et al. (2004) and Queirós et al. (2013):

$$\begin{array}{}\text{(1)}& {\text{BP}}_{\mathrm{C}}={\sum }_{i=\mathrm{1}}^n(B_i/A_i{)}^{\mathrm{0.5}}\times A_i\times M_i\times R_i,\end{array}$$

where B_i and A_i represent the biomass and abundance of species i within an assemblage, while M_i and R_i are categorical scores reflecting the mobility and reworking mode of species i.

The species and community irrigation potential IP_C was calculated following Wrede et al. (2018):

$$\begin{array}{}\text{(2)}& {\text{IP}}_{\mathrm{C}}={\sum }_{i=\mathrm{1}}^n(B_i/A_i{)}^{\mathrm{0.75}}\times A_i\times {\text{Bt}}_i\times {\text{FT}}_i\times {\text{ID}}_i,\end{array}$$

where B_i and A_i are again the biomass and abundance of species i, and Bt_i, FT_i, and ID_i are categorical scores describing burrow type, feeding type, and injection pocket depth, respectively.

Table S1 in the Supplement shows the categorical scores used for the calculation of the bioturbation and bioirrigation potential derived from Queirós et al. (2013) and Wrede et al. (2018). Species trait information was compiled from the CRITTERBASE data information system of the Alfred Wegener Institute (https://critterbase.awi.de, last access: 1 February 2025; Teschke et al., 2022). When trait data were unavailable, we followed the procedures proposed by Queirós et al. (2013), which suggest that genetically and morphologically similar taxa likely exhibit similar functional traits. Thus, functional scores were assigned based on the closest related taxon.

Both bioturbation and bioirrigation potentials represent relative (dimensionless), trait-based indices that integrate species abundance, biomass, and functional characteristics, rather than direct measurements of sediment reworking or chemical exchange measured from experiments. Although these indices inevitably overlook potential behavioral changes of the macrofauna across different sediments, macro-scale studies indicate that spatial fluctuations in abundance and biomass exert a greater influence on the final scores than local intra-specific variations (Queirós et al., 2013). Therefore, using standardized trait scores ensures comparable estimations of community functional potential and captures the major ecological signals across the study area.

2.3 Estimation of diversity indices

2.3.1 Taxonomic diversity

Taxonomic diversity of an assemblage was determined under the framework of Hill numbers (Hill, 1973), or the effective number of equally abundant species in a hypothetical community, which guarantees that the replication principle of ecology must be obeyed (Chao et al., 2010). The general equation of the Hill number, ^qD, is expressed as a function of the qth power sum of the relative species abundance,

$$\begin{array}{}\text{(3)}& {}^{q}D={\left({\sum }_{i=\mathrm{1}}^S{p}_i^q\right)}^{\mathrm{1}/(\mathrm{1}-q)}.\end{array}$$

The parameter q determines the sensitivity of the measure to the relative abundance of species. When q=0, ⁰D simply describes species richness, which counts all species equally. ¹D is mathematically equivalent to the exponential of the Shannon entropy, referred to as Shannon diversity. The measure for q=2 disproportionately favours dominant species that have higher relative species abundance, referred to as the Simpson index (Chao et al., 2014a). Here, we employed these three most commonly used diversity indices into the following analysis. Community evenness was calculated as Hill's diversity number evenness (E_1,0), which is represented as the ratio ${}^{\mathrm{1}}D{/}^{\mathrm{0}}D$ (Ricotta and Avena, 2003). The calculations were performed with iNEXT R package (Hsieh et al., 2025; Chao et al., 2014b).

2.3.2 Functional diversity

The functional diversity of macrofauna communities was calculated based on the species functional traits included in the calculation of BP_C and IP_C. Two functional diversity indices were selected to assess the relationship between functional diversity and ecosystem functions. Functional richness (FR_ic) was determined as the volume of the multidimensional convex hull defined by the traits of all species in an assemblage (Villéger et al., 2008). When only nominal and ordinal traits were available, FR_ic was calculated as the number of distinct trait value combinations in a community. The index was computed using the dbFD function in the FD R package (Laliberté et al., 2014). Functional redundancy (FR_ed) was calculated to demonstrate the degree to which multiple species within a community perform similar ecological functions. We followed the method described by Ricotta et al. (2016). First, community-level functional uniqueness (U) was quantified as:

$$\begin{array}{}\text{(4)}& U=Q/D.\end{array}$$

Here, Q expresses Rao's quadratic entropy (Rao, 1982), which denotes the average functional dissimilarity among individuals in an assemblage, while D is Simpson's diversity index, reflecting the maximum probability that two randomly selected individuals are dissimilar in their functional roles. Thus, U expresses the ratio of interspecific functional dissimilarity within a community. FR_ed was then calculated as the complement of uniqueness:

$$\begin{array}{}\text{(5)}& {\text{FR}}_{\text{ed}}=\mathrm{1}-Q/D=\mathrm{1}-U.\end{array}$$

2.4 Statistical analysis

Macrofauna biomass (wet weight) was converted to ash-free dry mass (AFDM) to accurately depict metabolically active tissue, using taxon-specific conversion factors provided by Brey et al. (2010). The resulting data were fourth-root transformed to reduce the dominance of high-biomass taxa and emphasize the contribution of low-biomass taxa. A Bray-Curtis dissimilarity matrix was calculated based on the transformed data and analyzed using hierarchical clustering with Ward's minimum variance method. To evaluate whether the structural differences captured by our biomass-based clustering led to variations in total community biomass, we conducted a one-way Analysis of Variance (ANOVA) with “Cluster group” as a four-level categorical factor. The ANOVA tested for significant differences in total community biomass (mean AFDM) across the four identified infauna clusters. Prior to modelling, normality and homoscedasticity of the biomass data were verified using the Shapiro–Wilk test and Levene's test, respectively. A distance-based redundancy analysis (dbRDA) was performed to identify the environmental variables that best explained biomass-based community dissimilarity. The initially selected environmental variables comprised water depth, modelled bottom salinity, bottom shear stress, and sediment-derived properties, including organic carbon content (OC), compositional percentages (sand, silt, and clay), grain size (D10, D50, D90), and statistical distribution measures (sorting, skewness, and kurtosis). To avoid severe multicollinearity within the dbRDA, a variable reduction protocol was implemented based on correlation analysis and Variance Inflation Factors (VIF). First, compositional grain-size variables (e.g., proportions of sand, silt, and clay that sum to 100 %) were excluded to avoid inherent compositional dependency and collinearity. Second, sediment compositional percentages (sand, silt, and clay) were omitted because they were highly collinear with median grain size (D50). Furthermore, depth was excluded because it was strongly correlated with both bottom salinity and bottom shear stress. Finally, VIF values were calculated for all remaining predictors to ensure that all selected variables dropped well below the conservative threshold of 5 (Zuur et al., 2010). Consequently, the optimized dbRDA model retained seven independent environmental drivers: sorting, skewness, kurtosis, D50, OC, salinity, and bottom shear stress.

The macrofaunal communities from each sampling station were mapped onto a single, unified non-metric multidimensional scaling (nMDS) ordination space based on their species composition. To evaluate functional trait distributions, we first calculated the Community Weighted Mean (CWM) for each functional trait modality across stations, using biomass as the weighting variable. These biomass-weighted functional traits were then fit onto the community-based nMDS space as passive vectors. The direction and length of each vector indicate its correlation strength with the ordination axes.

Calculated ecosystem functions, including bioturbation and bioirrigation potentials, community secondary production, and respiration, were compared among clusters using Kruskal–Wallis tests, followed by post-hoc pairwise comparisons to identify specific differences. Finally, Pearson's correlation was used to assess the relationships between taxonomic and functional diversity indices and ecosystem functions. Significant BEF correlations were further examined with interaction regression analyses to test whether these relationships varied along environmental gradients. When significant interactions were found, a simple slopes analysis was applied to evaluate changes in the relationships along the gradients. This approach allowed us to assess how environmental variables modulate BEF relationships. All statistical analyses were conducted in R version 4.4.3 (R Core Team, 2025) with package vegan (Oksanen et al., 2025) and interactions (Long, 2024).

3 Results

3.1 Characterization of benthic macrofauna communities

A total of 171 taxa, belonging to 9 phyla, 14 classes, 31 orders, and 89 families was identified from 31 086 individuals. Ash-free dry biomass across all stations ranged from 0.388 to 10.307 g per 0.1 m² (Fig. S1 in the Supplement). The Ward's cluster analysis based on ash-free dry mass matrices grouped the sampled communities into four distinct clusters. An ANOVA showed no statistically significant difference in AFDM between clusters (p=0.09, Table S2). A clear spatial pattern emerged when the clusters were projected into the study area (Fig. 1). Clusters 1 and 2 (represented by red and green symbols, respectively) dominated the deeper western region of the study area, whereas Clusters 3 and 4 (yellow and purple symbols, respectively) were mainly located in the shallower eastern region. All four clusters spatially overlapped in a transitional zone at a water depth of 20–30 m.

Communities from deeper areas were dominated in terms of AFDM by the ophiuroid Amphiura filiformis (Fig. S1), which contributed on average (± SD) 58.8 % ± 23.2 % to the total macrofauna biomass. On 5 out of 20 stations of Cluster 1, the decapod Upogebia deltaura contributed more than a quarter of the total biomass. In Cluster 2, the echinoid Echinocardium cordatum contributed over 20 % of the biomass at 6 out of 9 stations. Communities in the shallow areas of the muddy sediment region (Clusters 3 and 4) were mostly dominated by bivalves (50.0 % ± 19.9 % of the biomass), followed by polychaetes (26.2 % ± 14.8 % of the biomass). Additionally, Actiniaria (9.7 % ± 5.4 % of the biomass) were present at all stations within Cluster 3.

3.2 Environmental drivers of community clustering

Environmental variables exhibited spatial gradients with distinct patterns across the study area (Fig. S2). Higher values of D50 and D90 were observed in the deeper regions of the study area, consistent with higher sand content and lower proportions of silt and clay. Salinity exhibited an east–west gradient, with average values ranging from 25.7 to 32.3 and higher salinities occurring at deeper sites. In contrast, elevated bottom shear stress was stronger at shallower stations. Among the statistical measures of sediment grain-size distributions, kurtosis displayed a clear spatial pattern, with higher values associated with deeper sites.

In the dbRDA model, the selected environmental factors explained 44.1 % of the overall variation in community structure, with the first three dbRDA axes accounting for 64.5 %, 17.0 %, and 6.4 % of the explained variation, respectively (Fig. 2 and Table S3). A permutation test (n=9999) indicated that the overall model was statistically significant (pseudo-F=4.73, p<0.01). Projection of the four clusters onto the first and third dbRDA axes differentiated Cluster 1 from Cluster 3 and 4. It reflected a gradient in the community composition primarily along the East–West axis. Sorting, skewness, kurtosis, D50, salinity and bottom shear stress were identified as significant predictors of community composition (p<0.01, Table S3). The first dbRDA axis indicated that communities in shallow regions were strongly associated with higher bottom shear stress, whereas communities in Clusters 1 and 2, located in deeper areas, were more strongly associated with higher salinity, suggesting that these two environmental gradients shaped the benthic community. Greater kurtosis was particularly associated with Cluster 2, while sediments in Cluster 4 were distinguished by greater sorting, indicating that sediment characteristics significantly promoted the distinction between these two clusters.

Figure 2 Distance-based redundancy analysis (dbRDA) of community ash-free dry mass. Note that axes 1 and 3 are shown with axis 1 reversed to reflect the spatial orientation of the sampling area.

3.3 Functional trait summary

Functional traits involved in the calculation of bioturbation and bioirrigation potentials were projected onto the nMDS ordination based on their correlation with the ordination axes (Fig. 3), revealing a consistent spatial separation of macrofaunal communities across the study area (stress=0.14) despite some overlap between the clusters. Deeper stations tended to cluster on the left side of the plot, while those in the shallower region were clustered on the right. Within both depth zones, the communities were divided into two subgroups, with Clusters 1 and 2 situated in the upper-left and lower-left quadrants, respectively, while Clusters 3 and 4 fell in the upper-right and lower-right quadrants.

Figure 3 Non-metric multidimensional scaling (nMDS) ordination based on square-root transformed species ash-free dry mass. Functional traits (burrow type, sediment particle reworking mode, mobility, feeding type, and injection pocket depth) are projected as vectors according to their correlations with the ordination axes (stress=0.14). Red arrows indicate traits with significant correlations.

Trait modalities positioned at the extreme ends of the ordination axes with significant correlation are summarized in Table 1, providing insight into how functional differences resembled spatial community patterns. For the bioturbation trait describing sediment particle reworking, biodiffusors were positively associated with deeper stations, while surficial modifiers aligned more closely with the shallow communities. In terms of mobility, the four clusters demonstrated distinct trait modalities. Communities in the deeper area were distinguished by organisms with higher mobility. Cluster 1 was dominated by species that move freely within burrow systems, whereas Cluster 2 comprised slow but freely moving taxa. Taxa with lower mobility were associated with the shallower clusters. Cluster 3 was linked to taxa with limited movement, and Cluster 4 was further connected to tube-dwelling organisms that remain fixed within their tubes.

Table 1 Summary of trait modalities positioned at the extreme ends of the ordination axes with significant correlation.

For the bioirrigation trait that describes burrow types, organisms creating open irrigation systems were positively correlated with deeper stations, while epifaunal species or species exhibiting internal body irrigation (i.e. siphons) aligned more closely with shallow communities. Regarding feeding types, deeper communities were primarily composed of subsurface filter feeders, such as the dominant ophiuroids in those clusters. Cluster 4 tended toward deposit feeders, whereas Cluster 3 in the shallow region was represented by surface filter feeders, such as the dominant bivalves in these assemblages. Deeper stations were identified by species that produce injection pockets of 5–10 cm depth. The two shallower clusters showed contrasting patterns. Cluster 4 contained species that form shallow pockets of 0–5 cm depth whereas Cluster 3 was dominated by organisms that form injection pockets deeper than 10 cm.

3.4 Indices of ecosystem functions

Species richness across all communities ranged from 15–57 (Fig. 4a). The highest species richness was observed at a station of Cluster 3, while the species richness was lowest at a station of Cluster 4. Stations located in the transitional zone between 20 and 30 m water depth generally showed a higher species richness and evenness (Fig. 4d), indicating a more even macrofaunal composition in these communities not excessively dominated by single species. Communities at stations in deeper sections showed relatively lower values in both Shannon (Fig. 4b) and Simpson diversity (Fig. 4c).

Figure 4 Calculated taxonomic and functional biodiversity indices for each station. The darker color and larger size of the symbols represent higher estimated values.

The functional diversity indices showed similar spatial patterns as the structural indices. Functional richness was highest around 20–30 m water depth (Fig. 4e), likely reflecting the high species richness in this zone. The functional redundancy was highest at stations located in deeper regions (Fig. 4f), indicating that functional dissimilarity among individuals of that assemblage was low.

Community bioturbation potential averaged 284.0 ± 120.9 in the deep clusters and 151.2 ± 88.4 in the shallow clusters (Fig. 5a). Similarly, community bioirrigation potential averaged 223.2 ± 99.1 in the deep clusters and 86.6 ± 55.5 in the shallow clusters (Fig. 5b), indicating that both ecosystem functions occurred at rates roughly 2 to 2.5 times higher in the deep clusters. Correspondingly, the two deep clusters showed significantly higher bioturbation (H=18.3, p<0.001) and bioirrigation (H=25.3, p<0.001) potentials compared to the shallow clusters.

Figure 5 Calculated ecosystem functioning indices for each station. The darker color and larger size of the symbols represent higher estimated values. The unit of production and respiration is $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$.

Mean community secondary production was 107.5±69.0 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the deep clusters and 135.9±75.0 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the shallow clusters (Fig. 5c). Mean community respiration was 347.6±194.2 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the deep clusters and 505.9 ± 256.2 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the shallow clusters (Fig. 5d). Accordingly, respiration was significantly higher in the shallow clusters than in the deep clusters (H=4.98, p=0.025), while secondary production showed the same trend but was not statistically significant (H=1.92, p=0.166).

The bioturbation potential was significantly negatively correlated with Shannon diversity and Simpson diversity (both $r=-\mathrm{0.29}$) and with community evenness ($r=-\mathrm{0.55}$) (Table 2). Similarly, a negative correlation was found between the bioirrigation potential and the Shannon diversity ($r=-\mathrm{0.33}$), Simpson diversity ($r=-\mathrm{0.31}$) and community evenness ($r=-\mathrm{0.51}$) (Table 2). Community secondary production and respiration were significantly positively correlated with species richness (r=0.54 and 0.55), Shannon diversity (r=0.36 and 0.39), and Simpson diversity (r=0.32 and 0.34). Both functions had significant positive relationships with functional richness (r=0.54 and 0.56) but were negatively connected with functional redundancy ($r=-\mathrm{0.30}$ and −0.28).

Table 2 Summary of the correlation coefficients between diversity indices and derived ecosystem functions. P-values are indicated in parentheses. Bold style indicates statistically significant results (p<0.05).

3.5 BEF relationship modulated by environmental variables

The interaction analyses revealed that bottom shear stress and salinity strongly modulated the relationships between diversity indices and ecosystem functions (Table 3). Under conditions of lower bottom shear stress and higher salinity (i.e., in habitats of the deep clusters), both BP_C and IP_C were negatively correlated with Shannon and Simpson diversity. Under conditions of higher bottom shear stress and lower salinity (i.e., in habitats of the shallow clusters), these relationships became positive, albeit statistically not significant. The slope of the negative relationships between BP_C (or IP_C) and diversity indices further varied along the gradients of several environmental variables, including bottom shear stress, salinity, grain size (D50 and D90), and sediment clay content. Furthermore, the relationships between both secondary production and respiration with functional redundancy were also modulated by bottom shear stress and salinity. They were negative in habitats with higher salinity and low bottom shear stress but became positive under low salinity and high bottom shear stress. Additionally, this relationship was affected by sediment kurtosis. In habitats with high kurtosis, secondary production and respiration decreased significantly with increasing functional redundancy. Other significantly correlated BEF relationships were not modulated by environmental variables.

Table 3 Summary of the effects of environmental variable interactions on the relationships between diversity indices and ecosystem functions. The moderator values correspond to the selected points along the environmental gradient: −1 SD (1 standard deviation below the mean), Mean, and +1 SD (1 standard deviation above the mean). “+”, “−”, and “0” represent positive, negative, and no correlations, respectively. Significant effects (p < 0.05) are highlighted in bold. Detailed statistics are shown in Table S4, and the corresponding interaction plots are provided in Fig. S2.

4 Discussion

The benthic macroinvertebrate communities in an important depocenter for fine grained and carbon-rich sediment in the southeastern North Sea were investigated with high taxonomic and functional resolution and extensive spatial coverage. Biomass-based clustering revealed four distinct community clusters structured primarily along an east–west gradient across the area. Mobility and feeding strategies of benthic invertebrates showed consistent relationships with sediment properties, linking deeper clusters to subsurface-feeding and mobile taxa, while communities in shallower study area were associated with deposit and filter feeders and tube dwellers. The community potential for bioturbation and bioirrigation was negatively linked to benthic taxonomic biodiversity in the study region, whereas secondary production and respiration increased with taxonomic and functional richness, indicating that structural and functional diversity drive energy-related processes. Bottom shear stress and salinity likely modulated the observed biodiversity–ecosystem function relationships, disentangling the complex interplay among diversity, ecosystem functioning, and environmental variables, and demonstrating how benthic macrofaunal diversity regulates ecosystem functions linked to the cycling of organic matter and, thus, carbon in marine sediments.

4.1 Community structure and functional trait composition

The cluster analysis revealed contrasting trait compositions between communities in deeper and shallower parts of the study region. In the deeper parts, the clusters were dominated by biodiffusors with high mobility, such as the brittle star A. filiformis and the echinoid E. cordatum, which actively mix surface sediments over short distances (Kristensen et al., 2012). In general, biodiffusors can create extensive and homogeneous oxygenated zones in adjacent sediments, stimulating microbial aerobic respiration and meiofaunal growth (Michaud et al., 2005). For example, in laboratory experiments, A. filiformis at an average density of about 550 $\mathrm{ind}.{\mathrm{m}}^{-\mathrm{2}}$ mixed sediments with an intensity of approximately 35 cm² yr⁻¹ (Gilbert et al., 2003). Amphiura filiformis actively transports sediment along its arms both downward to the mouth opening (3–4 cm depth) and upward to the sediment surface (Persson and Rosenberg, 2003). Because of this efficient sediment reworking activity and due to its high dominance in terms of abundance (max. 245 ind. 0.1 m⁻²) and biomass, A. filiformis essentially contributes to the overall community bioturbation potential in deeper parts of the study area. In contrast, the communities in the shallower parts of the region were dominated by bivalves (e.g., Abra spp.) and polychaetes (e.g., Nereididae). Experimentally derived mixing rates under controlled incubation at 10 °C in surface sediments from the Gullmarsfjord (northeastern North Sea) indicate that Abra nitida exhibits a much lower baseline sediment-mixing intensity of approximately 0.5 cm² yr⁻¹ (Lindqvist et al., 2016).

Nereidid polychaetes dominating these shallow communities, such as Eunereis longissima and Nereis elitoralis, have been characterized as surface deposit feeders that feed preferentially from the openings of their built tubes, resulting in discretely motile behavior (Fauchald and Jumars, 1979). Similarly, observations of Alitta virens in large glass aquaria kept at 12 °C in both coarse sand and organic-rich mud from Conception Bay South, Newfoundland, Canada, demonstrated that these nereidids consistently perform shallow burrowing and near-surface deposit feeding across the entire 90 d observation period (Herringshaw et al., 2010).

In terms of traits related to bioirrigation, the dominant species in the deeper communities were primarily sub-surface filter feeders that create open burrow systems and relatively deep injection depths, potentially enabling intense bioirrigation and sediment ventilation associated with enhanced benthic biogeochemical activity. For example, A. filiformis contributes substantially to the total oxygen flux into soft sediments not only through its own respiration but also via diffusive processes across the additional sediment surface of its burrow walls. On slopes (65–90 m water depth) of the Skagerrak (western Sweden), A. filiformis accounted for at least 80 % of the total oxygen flux (Vopel et al., 2003).

The shallow communities were dominated by disturbance-tolerant surface modifiers with limited mobility (e.g., Abra alba, Ensis leei, Notomastus latericeus) and tube dwellers (e.g., Owenia fusiformis, Lanice conchilega, Lagis koreni), featuring as surface filter feeders. This trait is consistent with the influence of bottom shear stress, as suspension feeders tend to be associated with relatively higher flow velocities and enhanced particle delivery (Riisgård and Larsen, 1995), whereas subsurface suspension feeders primarily rely on their own pumping and irrigation behavior and may be less dependent on currents. These organisms primarily redistribute particles at the sediment surface (Bernard et al., 2019) and create surficial burrows, leading to lower bioturbation and bioirrigation (Christensen et al., 2000). For instance, deposit feeders ventilate primarily for the sole purpose of respiration but not for food acquisition (Renz et al., 2018). Therefore, continuous ventilation is not evident (Christensen et al., 2000) and less frequent compared to suspension-feeding organisms (Shull, 2009). The sand mason worm, L. conchilega, is known for its strong bioirrigation activity (Alves et al., 2017) and its pronounced influence on benthic respiration, nutrient release, and denitrification (Braeckman et al., 2010). However, the species was present only at a few shallow stations of the study area where it contributed at most 12 % to the benthic biomass. Traits associated with tube building can alter near-bed flow (e.g., Luckenbach, 1986) and support sediment stabilization (Van Hoey et al., 2008) and organic matter retention (Bailey-Brock, 1984), thereby promoting burial rather than remineralization of organic matter. Altogether, this clear contrast in traits between deeper and shallower communities provides insight into how functional traits and the derived indicators (i.e. bioturbation and bioirrigation) can be used to reflect the potential of benthic communities with regard to essential ecosystem functions related to carbon sequestration and sediment biogeochemical processes. In particular, macrofauna-driven particle reworking and porewater exchange can enhance the transport of oxygen and electron acceptors such as nitrate into deeper sediment layers (Aller, 1982), which affect microbial aerobic respiration (Glud, 2008) and nitrate reduction (Laverock et al., 2011). These processes regulate the degradation and remineralization of organic matter and, eventually, influence the balance between mineralization and long-term sequestration of carbon in marine sediments.

4.2 Environmental drivers of community variations

Sediment sorting, skewness, kurtosis, D50, salinity, and bottom shear stress clearly differentiated the four clusters in the dbRDA. Although water depth was excluded prior to modelling due to severe multi-collinearity with salinity and bottom shear stress, these remaining hydrodynamic and sediment characteristics captured the overall depth-related environmental gradient that shapes the benthic communities. The clusters in shallower parts of the muddy area were associated with high bottom shear stress. Bottom shear stress represents the frictional force exerted by currents on the seabed, which can lead to erosion, transport and redeposition of sediments. These processes induce major disturbance in soft-sediment habitats, altering sediment structure and causing mortality and displacement of resident fauna (Thistle, 1981; Sousa, 1984; Dernie et al., 2003). For example, Rowden et al. (1998) reported that the abundance of A. filiformis was negatively related to sediment surface current velocity, likely because of intense sediment erosion (Jago et al., 1993). Furthermore, elevated flow velocities and suspended particle concentrations can induce switches among feeding modes. For instance, A. filiformis can shift from deposit feeding to passive suspension feeding under active hydrodynamic conditions (Loo et al., 1996). This behavioral shift reduces their direct mechanical interaction with the seabed, which might further suppress community-level bioturbation intensity. The combined analysis of functional traits and environmental factors indicated that bottom shear stress, along with its profound influence on shaping sediment structure (e.g., grain size and kurtosis), suppressed the abundance of strong bioturbators in the shallow cluster, which ultimately led to depressed bioturbation and bioirrigation potentials.

Salinity is an important environmental factor structuring the macrofauna communities (Kröncke et al., 2011) as species differ markedly in their physiological tolerance to salinity. The muddy sediment region in the southeastern North Sea is located in a transition zone with a pronounced salinity front between saline North Sea water and relatively fresh coastal waters resulting from river inputs (e.g., Schrum, 1997; Chegini et al., 2020). The two deeper clusters, which had greater biomass but not higher species richness, were associated with higher salinity conditions. This pattern is consistent with findings from other coastal and estuarine systems. For example, Hampel et al. (2009) reported that macrofauna density, biomass, and diversity clearly differed between polyhaline (∼25) and euhaline (∼30.8) zones of the Schelde estuary (The Netherlands and Belgium), with higher salinity habitats supporting greater biomass but not higher species richness. Similarly, in Paranaguá Bay, Brazil, benthic invertebrate composition and abundance were strongly shaped by poly- and euhaline salinity gradients (Netto and Lana, 1996), although overall richness did not differ markedly within this salinity range. Together, these findings indicate that salinity acts as an environmental filter that constrains species occurrence according to physiological tolerance limits, resulting in shaping community composition and biomass rather than species richness. In our study, higher-salinity parts of the muddy sediment supported greater macrofaunal biomass without a corresponding increase in species richness.

The association between salinity and ecosystem functioning observed here does not imply a direct physiological control of salinity on bioturbation or bioirrigation. Instead, salinity indirectly influences ecosystem functioning through its role in structuring community composition and functional trait distributions. In the Baltic Sea, salinity was shown to be the most important environmental predictor for the presence of some groups of bioturbating Bivalvia and Polychaeta, especially A. alba, L. koreni and Nephtys caeca (Gogina et al., 2017). These taxa were also dominant in our dataset and critically contributed to bioturbation and bioirrigation. This indicates that higher ecosystem functioning at higher salinities emerges from the presence and dominance of key bioturbating species rather than from direct salinity effects on organismal behavior.

4.3 Biodiversity and ecosystem functioning

In general, no baseline values exist for bioturbation (BP_C) and bioirrigation (IP_C) potentials since these indices are dependent on community structure and therefore vary among communities. Our estimated values are of the same order of magnitude as those reported by Wrede et al. (2019), who also sampled muddy habitats in the German Bight. Although their study applied experimental treatments using artificial monocultures of organisms with different bioturbation behaviors, including A. filiformis and Nucula nitidosa, the resulting BP_C and IP_C values are comparable. The bivalve N. nitidosa contributes substantially to the total macrofauna abundance in the shallow clusters. The community potentials for both BP_C and IP_C were negatively correlated with taxonomic diversity but were not correlated with functional diversity. Additionally, community evenness was also negatively correlated with functional diversity. Together, these patterns indicate that BP_C and IP_C were disproportionately shaped by a few dominant functional groups or single dominant species rather than by overall trait richness within the community. The importance of keystone bioturbators, such as A. filiformis, in controlling bioturbation processes has previously been established at a larger spatial scale for the entire German Bight (Wrede et al., 2017; Weinert et al., 2022). Such dominance reflects a selection effect, in which ecosystem functioning is facilitated by the presence of particular highly efficient species rather than by diverse communities as a whole (Loreau, 2000).

The strength of the negative relationship between BP_C and IP_C and diversity varied across environmental gradients. In deeper parts of the muddy sediment region, where fewer powerful bioturbating species were present, the negative relationship was particularly pronounced. In contrast, this negative relationship became weaker or even slightly positive, though not statistically significant, in shallower areas. These spatial differences suggest that the presence and influence of key species depend strongly on local environmental conditions. Several environmental variables, including sediment grain size and clay content, appeared to modulate the relationship between BP_C and IP_C and community evenness, indicating that community composition and ecosystem functioning are essentially shaped by strong environmental filtering through sediment composition (van Dalfsen et al., 2000; Gray and Elliott, 2009). In particular, the filtering effect of the environment appears to constrain which dominant bioturbators are present and how effectively they perform, ultimately determining the strength of the BEF relationship in the muddy sediment.

Our estimations of macrofaunal secondary production are consistent with values reported for the North Sea along the UK continental shelf (Bolam and Eggleton, 2014), but are slightly lower, most likely due to differences in model implementation. Previous studies have shown that Brey's empirical model can underestimate secondary production, particularly for large-sized molluscs (Petracco et al., 2012; Beukema and Dekker, 2022). However, empirical approaches that incorporate population and environmental parameters have been shown to perform satisfactorily when applied at the community level (Dolbeth et al., 2005). Community respiration was about three times higher than secondary production in terms of amount carbon processed, revealing that a substantial fraction of assimilated carbon is allocated to maintenance metabolism rather than biomass production, which is a typical pattern of communities dominated by long-lived, taxa with low population turnover (e.g. A. filiformis; Muus, 1981). Model-derived macrofaunal respiration was well in the range but slightly above sediment oxygen consumption measured empirically in adjacent areas of the southeastern North Sea (Provoost et al., 2013; Oehler et al., 2015). This slight overestimation likely occurred because the ANN models used summer bottom water temperatures instead of annual means. These higher temperature inputs could accelerate the calculated metabolic rates, artificially inflating the estimated annual production and respiration (Brey, 2010). A mismatch between biomass-modeled respiration and measured total oxygen consumption in the sediment has also been reported in for a sandy coast (Charbonnier et al., 2016), estuaries (Herman et al., 1999), and a continental shelf (Wijsman et al., 1999). Biomass-based models estimate the potential community respiration derived from average metabolic rates, whereas in situ measurements may reflect potentially constrained conditions, under which macrofaunal activity and oxygen consumption may be suppressed (Horn et al., 2021). Consequently, direct agreement between modelled respiration and measured oxygen fluxes is not expected, although respiration remains a fundamental indicator of benthic metabolic demand and carbon cycling (Glud, 2008).

Experimental and modelling studies have predicted that secondary production and respiration are functions of biodiversity (McGrady-Steed et al., 1997; Duffy et al., 2017). In agreement with these predictions, secondary production and respiration of the benthic macrofauna community of the study region were both positively related to species richness. This is consistent with findings from other coastal ecosystems such as sandy shorelines, seagrass meadows, soft-sediment, and hard-bottom habitats (Rodil et al., 2020, 2022) as well as for the entire North Sea including a wide range of benthic habitats (Clare et al., 2022). Although species richness captures the number of taxa present, it does not directly display the functional roles species perform. However, increases in taxonomic richness often expand the range of trait combinations within a community (Petchey and Gaston, 2002; Hooper et al., 2005). Therefore, species richness and functional richness are typically strongly linked with each other. Single-replicate sampling might theoretically lead to a slight underestimation of functional richness by potentially missing rare species. However, given the environmental homogeneity of the study area, most rare taxa are likely to support common and redundant functions rather than unique ones (Mouillot et al., 2013). This close link implies that richness influences ecosystem functioning mainly by enhancing functional trait diversity, rather than through species number alone. Mechanistically, higher functional richness promotes greater niche complementarity (Loreau, 2000), as species with different trait combinations exploit distinct resources and increasing overall resource use efficiency (Wilson, 1999). Functional richness has been empirically demonstrated to enhance ecosystem productivity in several coastal marine systems (Dolbeth et al., 2015; Wang et al., 2024). Consistent with this mechanism, both secondary production and respiration in the muddy sediment area increased with functional richness, supporting the idea that trait diversity enables species to contribute additively to ecosystem functioning (Petchey, 2003). Thus, our results suggest that greater functional trait diversity promotes complementarity among species, enabling the community to sustain ecosystem functioning.

Functional redundancy (FR_ed) was negatively correlated with both community respiration and secondary production in the study area, although these relationships were relatively weak. The results from experimental studies suggest a positive effect of functional redundancy on ecosystem functions. For example, a mesocosm experiment revealed that secondary production of a seagrass bed in Chesapeake Bay at the US Atlantic coast was strongly depressed by the removal of a non-redundant amphipod species (Duffy et al., 2001). However, simplified ecological experiments involving a limited number of species may easily overestimate the effects of species removal and conclude that the loss of single, functionally non-redundant species leads to decline in ecosystem functions (Rosenfeld, 2002). A negative correlation between FR_ed and an ecosystem function may suggest that the overlap in functional niches and the associated lack of complementarity may intensify interspecific competition (Tilman et al., 2014), limit resource partitioning, and reduce overall energetic efficiency (Hooper et al., 2005). Community secondary production and respiration are not only measures of overall energy balance, but they are also useful in understanding various aspects of benthic dynamics and the impacts of environmental change (Dolbeth et al., 2012). Accordingly, understanding the relationship between functional redundancy and these ecosystem functions can be essential for evaluating the impacts of environmental variability.

Similar to functional richness, the relationship between FR_ed and secondary production and respiration was modulated by environmental variables, including bottom shear stress and salinity. Additionally, sediment kurtosis modulated the relationship with the slope shifting from negative to positive from the deep to the shallow stations of the study area. Kurtosis describes the “peakedness” of the sediment grain size distribution. Higher kurtosis typically characterizes low-energy depositional settings where sedimentation is stable and grain size distributions are narrowly peaked (Sahu, 1964). In contrast, strong hydrodynamic forcing frequently resuspends and mixes particles of different sizes, leading to poorer sorting and a more flattened grain size distribution (Wright and Nittrouer, 1995). The higher kurtosis at the deeper stations affected the negative relationship between FR_ed and secondary production (or respiration), whereas the effect became insignificant in areas characterized by low kurtosis. In previous studies, macrofauna dynamics could be explained by variations in kurtosis (Schrijvers et al., 1995; Lercari and Defeo, 2003; Bolam et al., 2008). Kurtosis could perform as a proxy for habitat stability, which, in turn, influences the strength of functional traits to translate into ecosystem processes. In deeper, high-kurtosis sediments where conditions are more stable, environmental filtering favours species with similar trait combinations, leading to functionally redundant communities, which may prioritize persistence (Biggs et al., 2020) rather than higher performance in ecosystem functioning. Therefore, the observed patterns indicate a potential trade-off between functional stability and ecosystem functioning efficiency, where highly redundant communities are more resilient but less conducive to ecosystem functioning under stable conditions.

5 Conclusion

Our study demonstrates that benthic macroinvertebrate communities in a depositional carbon sink in the southeastern North Sea are structured taxonomically and functionally by environmental gradients, even within an apparently homogeneous muddy sediment region. Energy-related ecosystem functions, such as secondary production and respiration, were positively associated with taxonomic and functional richness, while bioturbation and bioirrigation were dominated by a few key species, reflecting that different ecosystem functions are influenced by biodiversity through distinct pathways. The strength of these relationships varied along environmental gradients, including bottom shear stress and salinity, highlighting the role of local environmental conditions in shaping community structure and function.

Differences in species traits between deep and shallow clusters provide further insight into how BEF relationship are related to environmental gradients. Taxa found at deeper stations were characterized by traits associated with high bioturbation and bioirrigation potential. In addition, the co-occurrence of high biomass and production indicates that taxa in deeper clusters may effectively convert organic matter into benthic biomass, which might imply that carbon and nutrient cycling in these habitats are strongly biologically driven (Birchenough et al., 2012). In contrast, dominant taxa in shallower clusters exhibited traits associated with low bioturbation and bioirrigation potential, implying that carbon dynamics in these areas are largely governed by physical processes. In addition, we observed positive relationships between both taxonomic and functional diversity and community secondary production and respiration. Together, these patterns suggest that changes in biodiversity may have important consequences for benthic ecosystem functioning and, ultimately, for carbon storage capacity (Ramalho et al., 2020). Specifically, our calculations show that carbon processing in our study region is split between distinct ecosystem functions: shallower communities display high secondary production and respiration, whereas deeper communities are characterized rather by elevated bioturbation and bioirrigation potentials.

Functional diversity and trait-based indices offer mechanistic understanding of how community composition affects ecosystem functioning, making them useful tools for monitoring and managing continental margin systems and depositional hotspots like the muddy sediment region of the southeastern North Sea. However, a direct measurement that would allow for linking macrofauna community to carbon sequestration is still missing. Therefore, addressing this gap will be crucial for understanding how benthic communities regulate carbon cycling. The preservation of functional diversity is thus essential for both sustaining benthic ecosystem functioning and enabling these coastal systems to act as carbon reservoirs under continuous environmental change.