Benthic ecosystem functioning under climate change: modelling the bioturbation potential for benthic key species in the southern North Sea

Species data

The considered species data (presence, absence and species mean abundance per station, sampled with a Van Veen grab, 0.1 m²) were extracted from a data set collected for the EU-Project ‘Managing Fisheries to Conserve Groundfish and Benthic Invertebrate Species Diversity’ (MAFCONS) during the years 2003 and 2004. In total, 284 stations were sampled and used for the analyses. Not all countries which participated in the MAFCONS weighed the species individually. Therefore, the individual species mean wet biomass (g) were provided by Meyer et al. (2018) who compiled these data from different surveys in the study area. Our selection of the benthic key species Abra alba, Amphiura filiformis, Callianassa subterranea, Echinocardium cordatum, Goniada maculata, Nephtys hombergii, and Nucula nitidosa is based on previous studies with a focus on ‘ecosystem engineers’ and ‘bioturbation’ in the North Sea (Dauwe, Herman & Heip, 1998; Bouma et al., 2009; Braeckman et al., 2010; Birchenough et al., 2012; Braeckman et al., 2014b; Wrede et al., 2017; Meyer et al., 2018; Gutow et al., 2020). Details about the ecological and behavioural features and preferences of the species can be found in (Hartmann-Schröder, 1996; Ingle & Christiansen, 2004; Southward & Campbell, 2006; Zettler & Alf, 2021).

Species bioturbation potential

The BP_p is an index to estimate the species contribution to the biogenic mixing of sediments and was first described by Solan et al. (2004a). It is calculated as follows:

$$\mathrm{B}{\mathrm{P}}_{\mathrm{p}}=({\mathrm{B}}_{\mathrm{i}}/{\mathrm{A}}_{\mathrm{i}}{)}^{0.5}\ast {\mathrm{A}}_{\mathrm{i}}\ast {\mathrm{M}}_{\mathrm{i}}\ast {\mathrm{R}}_{\mathrm{i}}$$

B_i and A_i are mean body size (biomass in grams) and mean abundance (individuals 0.1 m²) of species (i) at a sampled station, while M_i and R_i are life traits (mobility, reworking). The latter are categorical scores according to Solan et al. (2004a) and allocated for over 1,000 macrofaunal species in Queirós et al. (2013) describing the mobility from one (for organisms that live in fixed tubes) to four (free movement via burrow systems) and reworking from one (for epifauna) to five (for regenerators). The BP_p was log₁₀ (x+1) transformed to down-scale large values. It was calculated separately for each species and sampled station. Finally, the total BP for all analysed key species (BP_t) was calculated as follows:

$\mathrm{B}{\mathrm{P}}_{\mathrm{t}}=(\mathrm{B}{\mathrm{P}}_{\mathrm{p}}1+\mathrm{B}{\mathrm{P}}_{\mathrm{p}}2+\dots +\mathrm{B}{\mathrm{P}}_{\mathrm{p}}\mathrm{n})$

The BP_t in our study is not equivalent to the BP_c in other studies, because we only consider seven selected species and not the entire community although we assume that these species represent a large proportion of the community bioturbation potential. Furthermore, the BP_t was only calculated for southern North Sea, because species with a core distribution in the northern North Sea were not included. Thus, total community BP for the entire North Sea would be misleading when only based on southern species, while modelling the BP_p needs to include the possible expansion towards northern regions.

Environmental parameters

The choice of environmental parameters is an important issue for the species distribution modelling. They can be grouped in three main categories: (1) limiting parameters that have eco-physiological influence, (2) natural and human-induced disturbances and habitat characteristics, and (3) resources, matter and energy consumed by species (Guisan & Zimmermann, 2000; Guisan & Thuiller, 2005). Six relevant environmental parameters were chosen (projected mean bottom temperature and salinity as well as depth, mud content, median grain size and peak wave stress), determining the spatial variability of the considered benthic ecosystem engineers (Callaway et al., 2002; Neumann, Ehrich & Kröncke, 2009; Reiss et al., 2010), to model and project BP_p and BP_t for the North Sea. The parameters were rasterized, with a resolution of 0.06 × 0.06 decimal degree and an extent of 60.4 north, 50.9 south, 3.1 west and 10.3 east (World Geodetic System 1984). Prior to the modelling process, raster cells covering land were excluded.

Mean bottom temperature (°C) and salinity of the simulated years 2001, 2050 and 2099 (Mathis & Pohlmann, 2014) were used to address climate change signals in the North Sea for the 21^st century. The annual means of volume-averaged temperature and salinity for 2001 are close to the climatological means for the period 1981–2000, so that the physical conditions in 2001 can also be regarded to represent mean conditions of the recent past. A global climate projection of the IPCC SRES scenario A1B was utilized to dynamically downscale potential climate change impacts for the North Sea to a meso-scale horizontal resolution of about 3 km (Mathis, Mayer & Pohlmann, 2013; Mathis & Pohlmann, 2014). Projections of mean bottom temperature for February and June were applied to consider seasonality and because the second quarter of the year was ecologically most affected by previous winter temperatures (Kröncke et al., 1998; Kröncke, Reiss & Dippner, 2013). After cold winters, a higher percentage of arctic-boreal species were found, whereas a positive North Atlantic Oscillation (NAO) index linked to warmer conditions facilitate an increase of species distributed in the southern North Sea (Kröncke, Reiss & Dippner, 2013). Simulations for the years 2001, 2050 and 2099 were chosen to project both, changes until mid-century and the end of the century to show changes over time. In general, mean bottom temperature for February and June was projected to change between −1.9 and 5.4 °C. An increase of winter temperatures (February) by 3 °C (2050) and by 3.4 °C (2099) were projected for the Danish coast, Skagerrak, central North Sea and German Bight (Supplemental Material S1). A decrease of winter temperatures by 0.5 °C (2050) was projected for the Southern Bight. During summer (June), a temperature increase by 1.8 °C (2050) and by 5.4 °C (2099) was projected for the southern North Sea (to a depth of 30–40 m) and to some degree for the Dogger Bank. A decrease of summer temperatures by 1.9 °C was projected for the Skagerrak and partly in the northern North Sea. Mean bottom salinity was projected to increase in summer by 1.0 (2050) and 1.4 (2099) for the Dutch and Belgian coast, the German Bight (only 2050) and the northern North Sea. In contrast, salinity was projected to decrease by 1.7 (2050) and 1.5 (2099) for the German Bight, the Dutch coast, the Oyster Ground and Skagerrak (Fig. S1a-i; Mathis & Pohlmann, 2014).

Depth data, with a resolution of one arc-minute (1.852 km), were derived from the General Bathymetric Charts of the Oceans (GEBCO) global bathymetry data set from the British Oceanographic Data Centre (GEBCO, 2003). Mud content and median grain size were collected during the North Sea Benthos Project (NSBP 2000) and the project Managing Fisheries to Conserve Groundfish and Benthic Invertebrate Species Diversity (MAFCONS). Data of peak wave stress, with a resolution of about 12 km, were provided by the Proudman Oceanographic Laboratory (Liverpool, UK), generated with the help of a 3-dimensional hydrodynamic model (Davies & Aldridge, 1993). Peak wave stress was calculated from a 1-year model run, covering the period September 1999 to September 2000.

More details regarding the environmental parameters can be found in (Reiss et al., 2011; Mathis & Pohlmann, 2014; Weinert et al., 2016). In contrast to bottom temperature and salinity, all other environmental parameters were kept constant over time.

The modelling approach

The spatial distribution of the BP_p for the seven species in the North Sea was modelled for the year 2001. The impact of climate change was projected for the years 2050 and 2099. For the analyses, the R environment version 3.5.1 (R Core Team, 2018) was used and applied in two steps: First, the R package ‘biomod2’ version 3.3-7 (Thuiller et al., 2009; Thuiller et al., 2015) was used to generate ensemble and consensus species distribution models, to model and project the probability of occurrence as well as presence and absence for each species. Biomod2 has the advantage, to combine several model algorithms to reduce uncertainties derived from the usage of different algorithms (Thuiller et al., 2009) and to use a consensus model based on weighted evaluation indices (Marmion et al., 2009) to separate between poor and reliable model results. Second, RF algorithm (Breiman, 2001) was chosen and the R packages ‘randomForest’ version 4.6-14 (Liaw & Wiener, 2002) and ‘caret’ version 6.0-80 (Kuhn et al., 2018) were applied to model and project the BP_p. For the second step, we followed the methodological approach of Hill et al. (2017).

Modelling and projecting the BP_p

For the first step, the R package ‘biomod2’ version 3.3-7 (Thuiller et al., 2009; Thuiller et al., 2015) was used to apply four model algorithms, which represent different approaches commonly used (i.e. regression models, classification trees and machine learning models): Generalized Linear Model (GLM) (McCullagh & Nelder, 1989), is an extension of the classical linear regression model, that allows the application of different distributed data like Gaussian, Poisson, Binomial, or Gamma. “GLMs are based on an assumed relationship (link function) between the mean of the response variable and the linear combination of the explanatory variables” (Guisan, Edwards & Hastie, 2002). Multiple Adaptive Regression Splines (MARS) (Friedman, 1991), is a flexible nonparametric regression model that combines piecewise linear functions (splines) and recursive partitioning to handle both linear and nonlinear relationship between dependent and independent variables. A forward/backward stepwise approach is used to generate the optimal model via cross-validation. Generalized Boosting Model (GBM) (Ridgeway, 1999), is an additive regression model consisting of many weak regression or classification trees, fitted in a forward, stagewise fashion (boosting) to estimate the final robust model (Elith, Leathwick & Hastie, 2008). Random Forest (RF) (Breiman, 2001; Prasad, Iverson & Liaw, 2006) is a classification and regression model, where classification trees (500–2,000) are drawn to maximum size without pruning based on resampled subsamples with replacement (bootstrapping) using a randomised subset of predictors. The average of all trees is calculated to aggregate the results. In this study, the default model options for the RF classification model were changed. The parameter ‘mtry’ (i.e., the number of variables randomly sampled at each node) was tested with values of 2, 3, 4 and 5 following Darr, Gogina & Zettler (2014) and the parameter ‘nodesize’ (i.e., the minimum size of nodes) was changed to 1, the default value for classification models (Liaw & Wiener, 2018). Thus, for each species, four model options were tested and the model with the highest evaluation index (mean AUC and TSS) was used to generate the RF species distribution model.

For model calibration and evaluation, the species occurrence data were randomly split into two subsets, 70% used for training and 30% for testing the model. Independent data were missing, thus 20 replicate runs were performed per algorithm (a sort of cross-validation) to account for algorithm variability (Thuiller et al., 2015). The area under the receiver operating characteristic (ROC) curve, often described as Area Under the Curve (AUC), and the True Skill Statistic (TSS) were used to evaluate the model results (Fielding & Bell, 1997; Allouche, Tsoar & Kadmon, 2006). Both measures are based upon the error matrix (true and false positive and negative rate). The AUC varies between 0 and 1, it is a threshold and prevalence independent index representing the relationship between the true positive rate (sensitivity) and the corresponding proportion of the false positive rate (1 – specificity). The TSS ranges between −1 and 1, it is a prevalence independent index and is defined as the sum of true positive rate (sensitivity) and true negative rate (specificity) minus 1, thus omission and commission errors (false negative rate, false positive rate) were accounted for (Allouche, Tsoar & Kadmon, 2006). Thresholds were applied for both evaluation indices (AUC ≥ 0.7 and TSS ≥ 0.4) to separate between poor and reliable replicate runs (Araújo et al., 2011). For each species only robust replicate runs were included in further analyses. The robust replicate runs were merged for each model algorithm in an ensemble model, separately for each species and year (2001, 2050 and 2099). An ensemble model is basically the mean of all projections with a single model algorithm. Thus, for each species analysed, three ensemble models (2001, 2050 and 2099) were generated with each model algorithm. Finally, robust replicate runs of all model algorithms were merged in a consensus model, again separately for each species and year. A consensus model is the weighted mean of the projections carried out with all model algorithms. The evaluation indices were weighted, so that a model result with a high evaluation index was more important in the consensus model. Finally, cut-off levels were applied to convert the consensus distribution maps of species probability of occurrence into presence/absence to adapt and apply the zero-inflated model approach from Savage, Lawrence & Squires (2015). For an explanation see the second step in the analyses. The applied cut-off levels maximised AUC and TSS, thus the most reliable distribution maps were generated (Jiménez-Valverde & Lobo, 2007; Barbet-Massin et al., 2012).

For the second step, RF algorithm was used and the R packages ‘randomForest’ and ‘caret’ were applied, to generate a separate RF regression model for each species. Thus, for each sampled station the calculated BP_p, the modelled and projected probability of species occurrence from the first step of the analyses and the values from the environmental parameters were used to model and project the distribution of the BP_p of the seven benthic key species in the North Sea under climate change. The evaluation metrics root-mean-square error (RMSE) and mean absolute error (MAE) were calculated by 10-fold cross-validation and indicate the average error of a model (Chai & Draxler, 2014). Both parameters have the same approach and measure the differences between prediction and actual observation. They indicate the average error of a model given in the same unit as the response variable, in this study the BP_p per raster cell (0.06 × 0.06 decimal degree). While the RMSE gives more weight to large residuals and is thus more sensitive to outliers in a small sample size (n < 100), the MAE gives equal weight to small and large residuals (Chai & Draxler, 2014). Both evaluation metrics can range from zero to infinite, where zero means no error of the model. The lower limit of 50 data points (sampled species presence) was defined to apply 10-fold cross-validation for the evaluation, which gives an average of five nonzero data points per fold for six out of seven species. For the burrowing shrimp C. subterranea less than 50 data points were available. Thus, 6-fold cross-validation was chosen to maintain the average of five nonzero data points per fold. A detailed description of the method is given in Hill et al. (2017).

Analyses of the bioturbation potential under climate change

The zero-inflated model approach from Savage, Lawrence & Squires (2015) was adapted and applied to generate the final BP_p distribution maps for the years 2001, 2050 and 2099. Therefore, each BP_p distribution map was masked with the correspondent binary distribution map of species presence/absence. Thus, BP_p values were assigned only for raster cells, which were modelled and projected as presence in the first step of the analysis. Vice versa, zero was assigned for raster cells, which were modelled and projected as absence. Afterwards, the BP_p values in the remaining raster cells were summed up to generate the BP_p in the North Sea for each year individually. This approach prevents over-projection and enables to generate reliable BP_p distribution maps. Finally, the BP_p distribution maps were added separately for the years 2001, 2050 and 2099 to calculate the cumulative total sum of the BP and generate the BP_t distribution maps. Bioturbation potential gain or loss was calculated by comparing the total BP_p and BP_t between the years 2001, 2050 and 2099. A general work-flow of the two-step approach can be found in Supplemental Material S2.

Model statistics

The selected species distribution models above the applied thresholds showed a good model performance with AUC values between 0.72 and 0.8 and TSS values between 0.43 and 0.56 (Table 1). They were used to generate the final consensus models of species probability of occurrence as well as the binary distribution maps of species presence and absence. The consensus models are based on all applied model algorithms (GLM, GBM, MARS, RF) for five out of seven species. For the other two, it is a consensus model based on GLM, GBM and MARS (N. hombergii), as well as an ensemble model based on GLM only (G. maculata; Table 1). For the second step of the analyses, the model performance of the BP_p distribution modelling with RF was evaluated with the parameters MAE and RMSE. Both parameters showed low values (MAE 0.06-0.15, RMSE 0.15-0.39) and the MAE was about two to three times smaller than the RMSE for the models of all species (Table 1), both indicating relatively small differences between prediction and actual observation. In addition, the MAE and RMSE were calculated in per cent for each species to put them in relation to the calculated mean BP_p per raster cell, based on the sampled data. The prediction of BP_p per raster cell was good. An error of <3% (MAE) was predicted for the models of all species and <5% and <8% (RMSE) for the models of six and one species, respectively (Table 1).

Table 1:

Overview of the observed species presence and absence, the included model algorithms, the model evaluation parameters as well as the importance of environmental parameters in the model.

Overview of the observed species presence and absence (p/a), the amount of models used to build the final consensus model with included models across all model algorithms, the included algorithms (with ¹⁾GLM, GBM MARS and RF, ²⁾GLM, ³⁾GLM, GBM and MARS) as well as the mean (m) and standard deviation (SD) of the evaluation parameters Area Under the Curve (AUC) and True Skill Statistic (TSS). The Random Forest regression model was evaluated with the indices mean absolute error (MAE) and root-mean-square error (RMSE). Both indices were also calculated in percent to present the average error of mean BP_p per raster cell. Results of the relative importance of environmental parameters (%) in the model (black = temperature June, orange = temperature February, sky blue = salinity June, yellow = depth, blue = median grain size, vermillion = mud content, reddish purple = peak wave stress).

Species	p/a	Models	AUC		TSS		MAE (%)	RMSE (%)	Variable importance (%)
			m	SD	m	SD
A. alba	53/231	551)	0.77	0.01	0.5	0.02	0.09 (2.25)	0.25 (7.38)
A. filiformis	132/152	141)	0.74	0.02	0.43	0.03	0.15 (0.28)	0.38 (0.94)
C. subterranea	35/249	681)	0.8	0.005	0.55	0.04	0.10 (0.49)	0.23 (1.30)
E. cordatum	116/168	241)	0.72	0.004	0.44	0.02	0.14 (0.08)	0.39 (0.28)
G. maculata	145/139	52)	0.73	0.01	0.43	0.03	0.06 (1.55)	0.16 (4.72)
N. hombergii	99/185	143)	0.72	0.001	0.43	0.02	0.06 (0.77)	0.15 (2.21)
N. nitidosa	53/231	771)	0.8	0.009	0.56	0.02	0.11 (1.20)	0.31 (4.07)

DOI: 10.7717/peerj.14105/table-1

BP_p and BP_t potential gain and loss

For four key species BP gain and for three key species BP loss was projected until 2099 in the entire North Sea (Fig. 1). Continuous BP_p gain until the end of the 21^st century was projected for the bivalves N. nitidosa (113%) and A. alba (198%), the bristle worm N. hombergii (176%) and the burrowing shrimp C. subterranea (246%) (Fig. 1). The predicted BP_p gain from 2050 to 2099 was even higher for A. alba and N. hombergii, as for the period from 2001 to 2050. However, together the species contributed only 12% (2001), 20% (2050) and 28% (2099) to BP_t (Fig. 2). In contrast, BP_p loss was projected for the brittle star A. filiformis (−42%), the bristle worm G. maculata (−7%) as well as the sea urchin E. cordatum (−8%) from 2001 to 2099 (Fig. 1).

Overall, BP_t for the southern North Sea only was projected to increase slightly until 2050 (+8%) and decreased slightly again until 2099 (Fig. 1). The sea urchin E. cordatum was the most important bioturbator and was projected to be responsible for over 60% of BP_t in all the analysed years (Fig. 2). The relative contribution of the sea urchin E. cordatum decreased moderately over time and accounted for 71% (2001), 66% (2050) and 62% (2099) of BP_t (Fig. 2). In contrast, the burrowing shrimp C. subterranea was predicted to account for the highest proportional gain (it contributed to 5%, 11% and 15% of BP_t modelled for 2001, 2050 and 2099, respectively) and the brittle star A. filiformis for the highest proportional loss (it contributed to 15%, 12% and 8% of BP_t) (Fig. 2). The bristle worm N. hombergii was predicted to almost double the proportional amount of BP_t from 2001 (4%) until 2099 (7%) (Fig. 2).

Spatial distribution of the bioturbation potential

Maximum values of BP_t per raster cell were projected for the southern North Sea only (Fig. 3), because species with a core distribution in the northern North Sea were not included. The highest BP_t was projected for the Oyster Ground and the southern edge of the Dogger Bank. Around this main distribution area, rather low values of BP_t were projected for the inshore regions of the southern and eastern North Sea, and the southern Bight (Fig. 3). The predictions showed, that BP_t increases moderately until 2050 (+8%) and decreased again until 2099 (−7%, Fig. 1). Nevertheless, the main distribution areas of the species BP_p were predicted to shift more and more to the north (see below), which goes along with a decrease of BP_t per raster cell in the southern North Sea (Fig. 3). This decrease was specifically pronounced in the regions with the highest BP_t, e.g., the Oyster Ground, although overall distribution of BP_t values in the southern North Sea remain relatively constant over time (Fig. 3).

The main distribution area with high values of BP_p of the different species in the southern North Sea was projected for the Oyster Ground (A. filiformis, C. subterranea, N. hombergii, and N. nitidosa), the southern Bight (A. alba, E. cordatum) and German Bight (G. maculata) (Figs. 3 and 4). Most species did not only show increase or decrease of BP values, but an expansion towards the northern and eastern North Sea (Figs. 3 and 4). Predictions of gain showed an increase of BP_p specifically in the northern North Sea for most species, but again related to a decrease of the maximum values of BP_p in the southern North Sea (Figs. 3 and 4). Predictions of BP_p loss go along with a decrease in the main distribution areas, especially for the brittle star A. filiformis (Fig. 4).

A temperature induced northward shift of southern North Sea bioturbators

Among the consequences of anthropogenic climate change, the rapid increase of sea surface temperature (SST) is likely the main factor influencing marine ecosystems to date (Hiddink, Burrows & García Molinos, 2015; Kröncke et al., 2019; Moore & Smale, 2020; Pilotto et al., 2020). Shallow coastal ecosystems such as the southern North Sea are drastically influenced by an increase in mean SST of 1.5–1.8 °C since 1950 (Beare et al., 2002). For the North Sea Mathis & Pohlmann (2014) projected an approximately linear increase in SST of about 2 °C per 100 years for the emission scenario A1B, with an increase in the near coastal zone of up to 5 °C (Schrum et al., 2016). A similar increase was found for the bottom temperatures used in this study for the southern North Sea with predominantly mixed water masses, while changes in bottom temperature in the stratified northern North Sea were somewhat less pronounced (Mathis & Pohlmann, 2014; Weinert et al., 2016). The changes in projected bottom salinity (Mathis & Pohlmann, 2014) were only moderate until 2099 (−1.7 to 1.4; Fig. S1l, S1o) and we assume that salinity changes have a minor effect on the distribution of the species and thus changes in species BP. Therefore, we focus our interpretation in the following on temperature changes as one of the most important drivers of distribution changes in the North Sea, although interactive effects with salinity are also conceivable.

In the present study, a latitudinal northward shift in BP_p was found for four out of seven species. This is in line with the results of Weinert et al. (2016), who found a latitudinal northward shift of 64% of the benthic species of the North Sea. All four species, namely C. subterranea, A. alba, N. nitidosa and N. hombergii, have in common that the BP_p in the core area of the southern North Sea decrease, while the overall spatial distribution expands and BP_p increases in the more peripheral areas of the North Sea.

In terms of direction, for abundances of the mud shrimp C. subterranea a southward shift and a habitat loss of 31% was predicted by Weinert et al. (2016). In the present study we found a decrease in BP_p in the core region (Oyster Ground, German Bight) from 2001 to 2099, but in combination with wider distribution of low BP_p intensities. This spatial expansion though leads to a remarkable overall increase in BP_p of C. subterranea of over 300% for the entire North Sea from 2001 to 2099. C. subterranea can be found with the highest biomass in areas with a mud content of >50% (Gutow et al., 2020), while higher abundances were found on a wide range of sediments (Coleman & Poore, 1980; Gutow et al., 2020). Thus, our results indicate that negative effects on the distribution and the concerning BP_p due to increasing temperatures may be compensated by the broad tolerance to environmental conditions of some species. However, this tolerance remains uncertain for C. subterranea with the strong preference of large adults for muddy sediments, which are not ideally represented in the environmental settings of our model.

The widespread predatory polychaete N. hombergii inhabits predominantly sediments containing fine sands with a mud content up to 20% in shallow waters (>30 m) and has a core distribution in the German Bight (Meißner, Darr & Rachor, 2008). Beukema (1990) stated that increasing winter temperatures may increase the abundance of N. hombergii in coastal and Wadden Sea areas, because it is negatively affected by cold winters. Indeed, we also found an increase of BP_p for the entire North Sea with warming, but a rather opposing effect on the bioturbation of N. hombergii in the southern North Sea, where BP_p decreased in the core distribution area. The temperature and salinity induced northward shift seemed to be amplified by the broad habitat requirements concerning sediment structure and adaptable feeding habits, while the decrease in the south indicates a warming beyond the species tolerance levels.

For the small bivalve species A. alba and N. nitidosa Weinert et al. (2016) predicted a habitat loss over 50% with distributional shift in north-western direction for N. nitidosa and in north-eastern direction for A. alba due to the expected temperature and salinity changes. Similarly, in the present study locally high BP_p of both species decreases, while the overall BP_p increased through a northward range expansion. They are detritus feeder and inhabit muddy and fine-grained sediments (Kröncke & Bergfeld, 2003; Van Hoey, Degraer & Vincx, 2004; Van Hoey, Vincx & Degraer, 2007; Kröncke et al., 2011). While N. nitidosa is not sensitive to cold winters, A. alba showed a decline or complete absence after cold winters such as in 1995/96. However, both species were found in higher densities after mild winters (Kröncke & Bergfeld, 2003; Reiss, Meybohm & Kröncke, 2006; Van Hoey, Vincx & Degraer, 2007), explaining the temperature-driven changes in our future projection.

An almost stable distribution with minor changes in BP_p were found for the polychaete Goniada maculata and the echinoderm Echinocardium cordatum. The latter species is actively reworking fine sediments and was found to be one of the most important bioturbators in the North Sea (Kröncke et al., 2004; Wrede et al., 2017). The relatively stable BP_p of E. cordatum with only a slight decrease is partly in agreement with results of Weinert et al. (2016), who found hardly any shift in the core distribution over time, but a loss in suitable habitat until 2099 for the species. The latter seem to be contradictory to our results, since the overall BP_p has not changed over time, but the BP_p values distributed more evenly over the shrunken distribution area of the species in 2099 (Fig. 4).

In contrast, the brittle star Amphiura filiformis revealed a loss in BP_p by 41% in particular in the south eastern North Sea. It burrows in sediments with moderate mud content and is as interface feeder in often stratified areas with low organic matter supply (Kröncke et al., 2004). Our result seems surprising since A. filiformis is an eurytherm species also distributed in southern regions, which indicates tolerance of warmer conditions. However, its core distribution is in the temperate and boreal regions (Kaschner et al., 2019) and similar to the other species with decreasing distribution in the southern North Sea, the projected warming of up to 5 °C in the southern North Sea seem to exceed the temperature tolerance of the species.

Changes in total bioturbation potential (BP_t) and ecosystem functioning

Changes in water temperature (here SST) and de-eutrophication were found to be main drivers for past long-term changes in the spatial community bioturbation potential (BP_c) variability in the southern North Sea, (Meyer & Kröncke, 2019; Meyer, Nehmer & Kröncke, 2019). While the BP_c is representing the community bioturbation potential (Solan et al., 2004a), we used and defined the total BP (BP_t) as a sum for the selected species to represent the summative response of the main bioturbators to warming. This can be justified by the finding that relatively few species account for most of the bioturbation potential for example in the German Bight (Wrede et al., 2017) and that the loss of single key bioturbators can have substantial consequences for ecosystem functioning (Lohrer, Thrush & Gibbs, 2004). Due to this focus on the selected species, diversity effects are not represented in the total BP. For example, it was shown that the importance of bioturbation and bioirrigation for biogeochemical cycling is reduced by interspecific competition for food or space (Hale et al., 2014). Thus, weak bioturbators, that were not considered in this study, could have a modifying effect on bioturbation and the resulting ecosystem functioning due to species interactions in a diverse community. However, these interactive effects wouldn’t been fully represented by the BP_c as a static trait-based indicator either.

Considering the analysed southern North Sea key bioturbators, the BP_t remained relatively stable until 2099 mainly because of opposing trends in BP_p among the species described above (Fig. 1). A high functional redundancy in the southern North Sea was already found for the BP (Meyer & Kröncke, 2019) and other biological traits (Shojaei et al., 2021). This compensational effect based on changes in the distribution of species might buffer against climate change impacts on ecosystem functions, but it remains unclear if and to what extent this is actually the case, in particular under the combined influences of other anthropogenic stressors (e.g., pollution, sand mining and dredging, bottom-contacting fishery). Some important bioturbators with a high contribution to the BP_c (based on their traits, biomass and abundance) such as A. filiformis or N. nitidosa were found to contribute comparatively little to the expressed biogeochemical cycling and sediment turnover (Wrede et al., 2017). Thus, realised and potential bioturbation might differ substantially and species might not be replaceable due to nuanced differences in activity (Hale et al., 2014; Queirós et al., 2015). Other species within the same functional guild but inhabiting southern warmer waters might actually compensate a loss in potential bioturbation when expanding their distribution range into the North Sea (Neumann et al., 2013; Zettler et al., 2018). While this aspect was not addressed in the present study, distribution modelling in general could indeed assist in investigating these large scale impacts of distributional changes and the role of non-native species (Reiss et al., 2015).

However, we found changes in the spatial distribution patterns and regional quantity of the BP_p on species level in the southern North Sea (Fig. 2), which is also reflected in the spatial changes of the BP_t. Several species showed a decrease in BP_p over time in the core distribution areas in the south (Figs. 3 and 4), while a moderate increase was observed in the peripheral and northern areas. The predicted long-term spatial changes of the present study are in line with the long-term trends in BP_c of macrofauna communities in the southern North Sea since the 1980s, which were related to SST increase mediated by the de-eutrophication processes (Meyer & Kröncke, 2019; Meyer, Nehmer & Kröncke, 2019) and possibly the decrease in primary production, e.g., (Capuzzo et al., 2018; van Beusekom et al., 2019; Desmit et al., 2020). In our study, a decline in BP_t was mainly projected for the depositional areas of the southern North Sea with predominantly fine, muddy sediments. In these less permeable sediments, bioturbation and bioirrigation by macrofauna are specifically important for the biogeochemical cycling in sediments. In contrast, in the more permeable sandy sediments, where we projected an increase of BP_t (e.g. the Dogger Bank), physical processes such as pore-water advection are much more important than biotransport (Neumann et al., 2021). Thus, BP might decline especially in areas where it is most important for biogeochemical cycling and ecosystem functioning.

Caution in the use of BP

While we found indications of bioturbation changes until the end of the 21^st century, there are some limitations in the use of theoretical BP_p as a predictor for bioturbation activities and as a predictor for expected changes in ecosystem functioning. The theoretical BP_p (Solan et al., 2004b; Queirós et al., 2013) is currently the only proxy for using existing data as a basis for a large-scale statistical modelling of bioturbation, because North Sea wide field measurements are not available. The BP_c successfully predicts the distance of particle transport (Queirós et al., 2015), while the prediction of bioturbation activity, depth, and the biodiffusion coefficient is limited. However, comparisons of the theoretical BP_c with the results from experimental approaches showed significant correlations (Gogina et al., 2017; Morys, Powilleit & Forster, 2017; Wrede et al., 2017), which supports the use of BP_c.

There are also limitations on the use of BP_c as an indicator in combination with projected changes of environmental variable and expected changes in ecosystem functioning. We did not account for changes in sediment structure or current regime, which were kept constant in the distribution modelling (Weinert et al., 2021). It is well known that some environmental parameters that can affect the spatial distribution of species in the North Sea often change simultaneously, such as temperature and food availability, which was not taken into account in our study either. Some studies though have already demonstrated temperature-related changes in the community structure and bioturbation potential of benthic organisms (Reiss, Meybohm & Kröncke, 2006; Clare et al., 2017; Meyer, Nehmer & Kröncke, 2019), which we utilize in our modelling approach. In addition to the geographical distribution of the species, temperature also regulates the metabolic rates (Brown et al., 2004). On the one hand, a combination of high temperatures and low food availability may lead to a reduced bioturbation due to the need to save energy and the concerning BP_c would be an overestimation (Maire et al., 2007, 2008). On the other hand, temporarily high temperatures can lead to higher bioturbation intensity due to increased activity, while a lower food availability may lead to increased bioturbation activity to ensure the uptake of enough organic carbon. Both processes would lead to an underestimated BP_c. Indeed bioturbation is more strongly linked to the metabolic rate than abundance due to the mechanistic relationship with respiration, feeding and moving (Cozzoli et al., 2021). Although metabolic rates were not considered in our study, the modelled BP is a product of abundance and biomass, with the latter more closely linked to the metabolism of organisms. However, the approach accounts for the direction and intensity of distributional changes of bioturbators as an important driver of changes in bioturbation (or biotransport in general) and the consecutive ecosystem functions (Cozzoli et al., 2021).