Impacts of habitat-specific benthic fishing compared to those of short-term induced variability by environmental drivers in a turbulent Baltic Sea environment

The short term impacts of fishing pressure were compared with the variability induced by environmental drivers on quantitative benthic community impact indicators. The different pressures were evaluated through comparative multifactor statistical analyses of their effects on macrofauna indicators in a Baltic Sea area with high natural disturbance. The area is exposed to a wide range of fishing intensities from long term non-fished areas to seasonal and annually frequently fished areas. Such evaluations are important for comparing the influence and short term variability of impact indicators from both demersal fisheries and the environment, including benthic community biodiversity (species richness), density (abundance in number of individuals), biomass, and average individual mean weight, with high spatio-temporal resolution. Environmental drivers include near bottom water current speed, salinity, temperature, and dissolved oxygen concentration, considering habitat specific and seasonal conditions. Demersal fishing-induced impacts were evident for all indicators. The highest fishing impacts were estimated in soft muddy and sandy habitats and in the second quarter of the year for all indicators, considering complex interactions. All environmental drivers, especially, current speed, had significant impacts on all indicators. The significant influences and short term variability caused by environmental drivers were of the same or larger order of magnitude as fishing impacts. Consequently, the short term influence of environmental drivers and seasonal differences in fishing pressure need to be considered when using quantitative benthic fishing impact indicators and identifying areas that are more or less resilient to fishing in relation to short-and long-term fisheries management plans.


Introduction
Fisheries on continental shelf areas can be considered one of the most important and widespread anthropogenic activities affecting marine ecosystems (Kaiser et al., 2002;Nielsen et al., 2018).Bottom trawling and seining are frequently used fishing methods with gears that cause physical changes in the seabed and marine benthic habitats and their community structures (Grieve et al., 2014;Eigaard et al., 2016;Depestele et al., 2016;Gerritsen et al., 2013).Several types of such effects of fishing on the benthic invertebrate habitats and communities have been documented in the literature (e.g.Jennings and Kaiser, 1998;Hall, 1999;Collie et al., 2000;Kaiser et al., 2002Kaiser et al., , 2006;;Kaiser, 2019;Hiddink et al., 2006aHiddink et al., , 2017;;Buhl-Mortensen et al., 2013;McLaverty et al., 2020;Rijnsdorp et al., 2020) as further specified below and in the Supplementary Material (SM) S.6.Studies have mostly focused on the direct and long-term effects and impacts of demersal fishing on the benthic community and habitats (e.g.Kaiser et al., 2002;Rijnsdorp et al., 2016Rijnsdorp et al., , 2020;;Clark et al., 2016;Sköld et al., 2018;Hiddink et al., 2019Hiddink et al., , 2020) ) while the impacts of fishing compared to the influence caused by the intensity and variability of environmental drivers in the short term are less well-investigated (e.g.McConnaughey and Syrjala, 2014; van Denderen et al., 2015;Szostek et al., 2016;Lambert et al., 2017).

Fishing pressure
In recent years, robust fishing pressure indicators for towed bottom gears have emerged and were initially tested for implementation in different seabed sediment types and benthic habitats using standardised methods and widely applicable approaches and perspectives (Kaiser et al., 2002;Gerritsen et al., 2013;Depestele et al., 2016;Rijnsdorp et al., 2016Rijnsdorp et al., , 2020;;Eigaard et al., 2017;Hiddink et al., 2020).The details of the indicators are described in SM S.6.It has been shown that the impact of fishing depends on the fishing pressure from the type of gear, speed, and mass and size of all the individual gear components, such as otter boards and ground ropes, as well as the fishing frequency and sensitivity of the benthos (e.g.Grieve et al., 2014;Eigaard et al., 2016;Rijnsdorp et al., 2016).

Types of impacts of fishing on benthic invertebrate communities
In the long term, fishing can induce changes in benthic invertebrate communities, which play an important role in the functioning of marine ecosystems (e.g.Collie et al., 2000;Kaiser et al., 2002;Buhl-Mortensen et al., 2013).Fisheries can affect ecosystem functioning, such as benthic-pelagic coupling, processing of organic carbon, and re-mineralisation of nutrients (e.g.Thrush et al., 2001;Olsgard et al., 2008;Rijnsdorp et al., 2016).Other studies have documented the long-term functional impacts of fisheries on benthic habitats and invertebrate communities (e.g.Bremner et al., 2003;Tillin et al., 2006;De Juan et al., 2007;Lambert et al., 2014;Sciberras et al., 2018).Bottom trawling has been shown to reduce benthic production, biomass, and species richness (Collie et al., 2000;Kaiser et al., 2006;Worm et al., 2006; van Denderen et al., 2014a;b;Sciberras et al., 2018;Hiddink et al., 2020).The response of benthic ecosystems to bottom trawling depends also on food-web interactions which have long term implications, as detailed in SM S.6 (e.g.van Denderen et al., 2013;Sköld et al., 2018).Current research and literature show that trawls may damage or kill benthic organisms, modulated by the depth of penetration and the sturdiness of the benthos.It also show that the effect on the population is modulated by the intensity of trawling and the recovery rate of the benthos, where the latter is related to the life span of the organisms (see e.g.Rijnsdorp et al., 2020).
Benthic invertebrates show differences in their vulnerability, partly because of the strong seasonality in their life cycles (e.g.Zettler et al., 2000;Gogina et al., 2014Gogina et al., , 2017) ) and the different levels of trawling throughout the year (e.g.van Denderen et al., 2014b;Bigné et al., 2018).The recovery rates of benthic invertebrates from fishing impacts have been shown to be higher in areas with high natural disturbances (Lambert et al., 2014).Accordingly, the short-term effects of fishing can vary with seasonal recruitment and settling processes for benthic invertebrates, as well as short-term seasonal and interannual variabilities herein.Szostek et al. (2016) found no significant short-term effects of scallop dredging across a gradient of fishing intensity, whereas significant short-term influences from concurrent key environmental drivers (which vary seasonally) were found on the benthic community across benthic impact indicators on different substrate types in the English Channel.Similar and same order of magnitude effects of bottom trawling and natural disturbance from tidal-bed shear stress were reported by van Denderen et al. (2015) on the composition and function of benthic communities across habitats in the North and Irish seas.Both factors caused declines in long-living, hard-bodied, and suspension-feeding organisms, and the benthic response to trawling appeared to be smaller or absent in areas exposed to high natural disturbance.Lambert et al. (2017) found that the benthic community structure, biomass, and abundance at the population level were resilient to fishing in the short term, based on experiments and quantitative analyses, and that natural temporal variation in those community metrics exceeded the effects of dredge fishing at a dynamic study site.In addition, McConnaughey and Syrjala (2014) found no significant short-term effects of bottom trawling compared to those of environmental drivers in the form of heavy storm events on soft bottom fauna density and biodiversity in previously untrawled areas in the eastern Bering Sea.As such, the short-to medium-term influence of environmental drivers and demersal fisheries on several benthic ecosystem parameters and metrics can vary considerably and be very different throughout the year.At the same time, when evaluating the impact of fishing on benthos, it is important to consider that typical fishing intensity shows clear seasonal differences (e.g.Bigné et al., 2018) which were also found in the current study.

Environmental drivers and other impacting factors
Apart from bottom trawling, marine benthic faunal communities are driven by physical and chemical environmental factors.Among the important environmental factors affecting the distribution of invertebrate communities, hydrographical conditions, such as benthic water temperature, oxygen content, depth, and current speed (e.g.Buhl-Mortensen et al., 2013;Zettler et al., 2017;van Denderen et al., 2020), suspended material (e.g.Linders et al., 2018), and chemical compounds (e.g.Bradshaw et al., 2012), highly determine the benthic community composition, structure, and functioning.In the Baltic Sea, the salinity gradient is also an important factor that determines the occurrence of benthic invertebrates (e.g.Zettler et al., 2007Zettler et al., , 2017;;Bossier et al., 2018Bossier et al., , 2020Bossier et al., , 2021)).Studies in the southern North Sea have revealed small-scale heterogeneity in bathymetry and sediment composition with alternating ridges and troughs at scales well within 1 min x 1 min grid cells (van Dijk et al., 2012;Koop et al., 2019;van der Reijden et al., 2019), and differences in marine benthic biomass and species composition have been observed between the valleys and crests of sand waves due to small-scale hydrodynamics (Ramey et al., 2009;Rijnsdorp et al., 2016).Similarly, burrows and tubes from infaunal invertebrates have been shown to affect sediment characteristics, faunal assemblage structures, and near-bed hydrodynamic conditions (Bolam and Fernandes, 2003;Braeckman et al., 2014;Dame et al., 2001;Rabaut et al., 2007).In addition, anthropogenic pressures other than fishing, such as eutrophication, sediment extraction, shipping and transport, energy platforms, pipelines, and cables, have been shown to impact the Baltic Sea benthic habitats both in the short and long term (e.g.EU-HELCOM-BalticBOOST, http://www.helcom.fi/helcom-at-work/projects/completed-projects/baltic-boost; Korpinen et al., 2018).Similar to environmental drivers, anthropogenic pressures can be subdivided into continuous and diffusive influencing pressures (e.g.eutrophication and change in climate) and more abrupt and distinct point pressures in time and space (e.g.fishery and gravel extraction).

Impact indicators for benthic effects of fisheries
The sensitivity of benthic habitats to various fishing types and pressure levels differs greatly, and impact evaluation requires standardisation of methods (e.g.Hiddink et al., 2017Hiddink et al., , 2020;;Eigaard et al., 2016Eigaard et al., , 2017)).Therefore, there has long been an interest in identifying and comparing indicators of benthic invertebrates that would show a clear response to human-induced effects, such as the Shannon-Wiener diversity index, AZTI Marine Biotic Index (AMBI), Benthic Quality Index (BQI) (e.g.Zettler et al., 2007), Pielou's evenness index, and Simpson's index (e.g.Szostek et al., 2016).However, despite several attempts (e.g.Hiddink et al., 2006bHiddink et al., , 2020;;Zettler et al., 2007;Rijnsdorp et al., 2016Rijnsdorp et al., , 2020;;Eigaard et al., 2017) current research has not yet been able to establish and agree upon a set of overall robust, generally applicable, and adequately sensitive impact indicators.Standardised indicators should ideally be applicable to the full benthic invertebrate community and different habitats across ecosystems, be exclusively quantitative, and possess thresholds according to different fishery pressure levels for all demersal fisheries across different fishery systems (e.g.Hiddink et al., 2020).It is also desirable for indicators to capture the long-and short-term impacts differently, with specific thresholds taking into consideration different levels of resilience and recovery rates, as well as the variability and dynamics caused by the influence of environmental drivers (e.g.ICES, 2019a; b).Furthermore, the efficiency of different types of impact indicators according to coverage, robustness, and sensitivity to several types of simultaneously influencing pressures and drivers has been intensively discussed (e.g.ICES, 2019a).The present study contributes to this discussion by using the findings of Hiddink et al. (2020) that benthic community abundance (density) and biomass are the strongest indicators of trawling impacts, followed by species richness, and that the effects of trawling are more pronounced on whole-community numbers and biomass than on individual taxa, as a stepping-stone.

Aim of present study
The present study investigated habitat-specific benthic impacts of fishing pressures compared to those of short-term induced variability by environmental drivers.To address this, it is essential to consider the interactions in the impacts on benthic communities between various pressures and drivers over different time scales and habitats.In this case study, we focused on the marine environments of the Western Baltic Sea with high natural disturbance levels, pronounced gradients of environmental drivers and high variability in fishing pressure.
The aim of the present study was to evaluate the ecological effects of fisheries with towed gear by considering the following impact indicators: i) biodiversity (species richness, number of species), ii) density (total number of individuals), iii) biomass, and iv) average individual mean weight in the whole benthic community.These biological indicators were tested against different levels of fishing pressure and physical environmental driver levels and conditions to compare influences and short-term variability in impact indicators from both the demersal fishery and environment.

Benthic invertebrate community data
Benthic macrofauna and other community data were sampled and mapped under the Environmental Impact Assessment (EIA) of the fixed link between Denmark and Germany in the Fehmarn Belt area of the western Baltic Sea (FEMA, 2013; see details in SM S.1).In total, 276 samples of epifauna and infauna collected with Van Veen Grabs operated from a ship were included in the present analyses (Fig. 1; Table 1; SM S.1).After collection, the samples were sieved through a 1-mm sieve.The retained material was preserved in 4% formaldehyde-seawater solution until further taxonomic identification, enumeration, and biomass estimation in the laboratory (FEMA, 2013).
The standardized repeated grab sampling during 2009 and 2010 covered different seasons (2nd and 3rd quarter), benthic habitat types, depths, and stations located throughout the study area (Fig. 1; Table 1).Repeated sampling was conducted at the same stations in the spring (second quarter) and autumn (third quarter) in both 2009 and 2010, and the samples were collected over a period of approximately one month (see further details in SM S.1), which supports the choice of a 3-month SAR prior to each sampling event as an appropriate proxy for shortterm fishing pressure (see also Section 2.2.below).The distributions of the sampling across strata are shown in Table 1 and Fig. 1 and are further detailed in SM S. 1 and FEMA (2013).The van Veen Grab sampling did not cover the coastline sampling stations, but only the deeper located stations.Here a ship could be operated for the grab sampling without the propeller disturbed the sediment to ensure undisturbed samples (Fig. 1).
The sampling area of the van Veen Grab was ~0.1 m 2 , but as the different van Veen Grabs used had slightly different observation coverage areas (ranging between 0.098 and 0.1166 m 2 ), an area conversion factor by sampling gear was calculated and used to standardise and convert to the same surface area for the density (i.e.number of individuals) and dry weight (see FEMA, 2013) per species for each sample, that is biomass and individual mean weight, but not the biodiversity parameter (see details on standardisation in SM S.2 and Table S.2).However, biodiversity is very closely and approximately standardised and comparable to the other parameters (indicators) because biodiversity per grab sampling station was used in the analyses.A full list of invertebrate species and taxonomic groups sampled by van Veen Grab was established (SM S.1, Table S.1), and for each van Veen Grab sample, the standardised benthic community parameters (indicators) estimated and analysed per sample were: i) biodiversity (BD) measured as the total number of species; ii) density (D) (i.e. total number of individuals); iii) biomass (B) as the total dry weight (see FEMA, 2013) per taxon (i.e.species or higher taxa); and iv) average individual mean weight W (=biomass/density) integrated over all species.

Fishing pressure data
The benthic sampling stations were comprised of a wide range of fishing pressures across years, seasons, and habitats, including stations with true zero fishing pressure values (Table 1).Fishing pressure (FP) data were estimated for both the benthic sampling period from 2009 to 2010, and a longer period from 2008 to 2013 (Table 1; Fig. 2) for which international VMS data were available (see below).The FP data were based on both Danish and German VMS data (vessel satellite monitoring) available for vessels with lengths measuring 15 m and above using demersal hauled gears (mainly otter board trawlers and otter board pair trawlers, but also seiners and dredgers).Fishing effort was extracted from national VMS databases and processed by the DTU Aqua (Danish fishery data) and Thünen Institute (German fishery data).The VMS raw data were coupled to fishery logbooks and sales slip data, according to unique identifiers based on position and time, by applying and combining the methodologies described in Bastardie et al. (2010), Hintzen et al. (2012), and combined in Eigaard et al., (2016Eigaard et al., ( , 2017)); Eigaard et al. (2016) included questionnaire surveys of fishermen and net makers with estimates of the dimensions of the different gears.
The relationships between gear dimensions and vessel size (e.g.trawl door spread and vessel engine power, kW) for different gear groups were used to assign quantitative information of bottom contact (e.g.width of gear) to each logbook trip using linear models to describe the relationship between gear width and vessel size akin to Eigaard et al. (2016), and the extended logbook data were combined with interpolated vessel tracks based on VMS data using the methodology described in Hintzen et al. (2012).This enabled statistical modelling of the vessel size or engine power ~ gear size relationships for different métiers (combinations of gear types and target species) to deduce the width of the sweep

Table 1
Overview of the benthic invertebrate samples and categories from Van Veen Grab sampling used in the analyses including their temporal coverage.Furthermore, Fishing Pressure (FP) expressed as categorized Swept Area Ratio categories (SAR, km swept / km2) by quarter by category is indicated, and the average, minimum and maximum FP observed at stations in given category, i.e. the FP range included in the analyses.In the lower part of the table the range (min, max, mean) of number of species (biodiversity, BD) and number of individuals (density, D) per sampling in each category.The habitat types are 1: Sublittoral mud, 2 Sublittoral sand, and 3 Sublittoral mixed sediments.of each of the VMS-interpolated fishing events taking place across the benthic sampling stations.Station-specific FP was defined as the cumulative swept area within a radius of 1000 m around each benthic sampling station divided by the total circular area (i.e.swept area ratio, SAR; sensu Eigaard et al., 2017).Fishing pressure in the previous three months before the respective benthos sample month was calculated and used as short-term FP (see also Section 2.1).Such a high spatial scale resolution is necessary to obtain accurate estimates of the fishing footprint (Amoroso et al., 2018).This spatiotemporal resolution is also appropriate in relation to the van Veen Grab sampling design with respect to the sampling frequency and spatial scale between stations.The 0-FP-areas and stations observed from 2009 to 2010 are also long-term low fishing pressure areas and stations observed for the period from 2008 to 2013 (Fig. 2).This period provides a good representation of the long-term fishing pressure, trends, and non-fished areas because haul allocation is very conservative according to fishermen's electronic clear tow information, which is also the period in which Danish, German, and Swedish VMS data are available to estimate this.An overview of the ranges (maxima and minima), means, and contrasts according to strata (year, season, habitat type) of the FP and the dependent indicators are shown in Table 1.

Physical habitat data: hydrographical data and EUNIS level 3 benthic habitat data
Environmental parameters from hydrographical data, bottom depth, and sediment physical characteristics were obtained from two databases and a physical hydrodynamic model processing data available from the Danish Meteorological Institute (DMI, www.dmi.dk).The physical data were produced by a Baltic-North Sea ocean-ice model HBM (HIROMB-BOOS Model) in the operational setup of the DMI.A biogeochemical module (ERGOM) is dynamically embedded in the HBM.The HBM is a three-dimensional, free-surface, baroclinic ocean circulation, sea ice model.The model allowed for full two-way nesting of grids with different vertical and horizontal resolutions, as well as time resolution (further details on the models and their coupling, resolution, and environmental parameters are given in SM S.3, and in Bossier et al., 2018 and references herein).The extracted hydrographical parameters estimated and analysed were near seabed (0.5 m above seabed) including temperature (t, • C), salinity (s, psu), oxygen concentration (o, mmol O 2 /m 3 where 1,0 mmol O 2 /m 3 = 0.032 mg O 2 /l), net bottom current speed (u, m/sec), and bottom depth (m).Monthly minima and maxima, as well as the daily mean values for these parameters, were extracted for all benthic sampling station positions according to the sampling time.An overview of the ranges of the physical parameters and their contrast according to strata (year, season, and habitat type) is given in Table 2, and area-specific ranges of current speeds are shown in Fig. 3.For the statistical analyses, we used the monthly minima for temperature, salinity, and oxygen concentrations, and the monthly maxima for current speed associated with the sampling month for each of the van Veen Grab samples included in the analyses.The narrow Fehmarn Belt is a highly turbulent hydraulic bottleneck for water exchange to the Baltic Sea because the volume transports between Kattegat and the Arkona Basin are usually split between the Belt Sea and the Sound with a ratio between 70% and 30% (Jakobsen and Trébuchet, 2000;She et al., 2007;Jakobsen et al., 2010) (Fig. 3).
The classification of benthic habitats for European seas is available at https://www.emodnet.eu/en/seabed-habitats.The EUNIS Level 5.3 habitat classification considers depth, sediment grain size, light, and level of disturbance by hydrodynamic forces (Davies et al., 2004).Seafloor sediment data, together with depth information, were extracted from the EUNIS level 5.3 databases, which was processed and compiled for the benthic sampling positions.Three EUNIS level 3 habitats at the location of benthic sampling were relevant (Fig. 1): 1 Sublittoral mud (A5.3), 2 Sublittoral sand (A5.2), and 3 Sublittoral mixed sediments (A5.4).Some stations did not have sediment data available and they were categorised as 4 'outside polygons' and this category only included 2 samples where FP was above 0 (not used).

Statistical analyses and model for multivariate analysis of variance with a mixed GAM model
Finally, the FP data were merged with benthic invertebrate and physical habitat data using unique identifiers of date, position, and station number to establish an integrated dataset for statistical analyses.A general additive model with mixed effects (GAMM) was developed for the statistical data analyses (Table 3; see the reasoning for model choice in SM S.4).The dependent variables on the benthic community level were the indicators BD, D, B, and W, respectively, and the explanatory variables covering anthropogenic and environmental drivers were based on each benthic sample FP, habitat type, year, season of year, and environmental hydrographical parameters, as described above.Alternative model versions were compared by indicator and tested with interaction effects (Models 1-4, Table 3) and main effects (Models 5-8, Table 3) between the explanatory variables.The interactions chosen in the models focused on identifying short-term interactions between fishing pressure, season, and habitat type as one of the main aims of this study.The models were reduced by identifying significant effects using backward elimination of insignificant model terms, and their statistical significance was estimated (Table 4; Fig. 4; see further details in SM S.4).The stepwise backward elimination of insignificant model terms was carefully checked at each step so that the remaining significant parameter estimates did not change very much in this model reduction process.This model reduction process resulted in the final reduced models representing the model estimates shown in Table 4, where only the significant variable estimates which exactly represent the reduced models are given, that is non-significant parameters have been reduced away (see also details in SM S.4).The model estimates of the dependent variable according to the impact and significance of the explanatory effects were provided (Table 4; Fig. 4), and the overall variability of the data

Table 2
Overview of ranges of the environmental parameters and their contrast according to strata.Yearly, quarterly and habitat specific minima, maxima and averages measured at all van Veen Grab benthic sampling stations.In the statistical analyses we used the monthly minima for temperature, salinity and oxygen concentrations and the monthly maxima for current speed associated to the sampling month for each of the van Veen Grab samples included in the analyses.explained by the model was estimated (Table 4).The model can handle unbalanced data, and raw input data were analysed accordingly.In addition to raw model estimates, the model estimates were further standardised (Table 4) by dividing the model estimates with the global observed ranges between observed maxima and minima in own parameter units of the quantitative explanatory variables (FP and environmental parameters).This makes the dependent variables, explanatory variables, and the model intercept and slope coefficients of fixed effects comparable in relation to the van Veen Grab benthic sampling stations in a meaningful and statistically correct way.This allows the relative differences between standardised estimates to be compared to each other for the different models for each indicator according to the different explanatory variables.We do not apply absolute estimates in any respect to e.g.calculate thresholds, explain impacts and their causality, or provide advice, but only make relative comparisons between standardised model estimates of relative directions (negative, positive) and level of magnitudes.If, for example, one wants to compare the strength of the effect of FP for sediment type 1 in quarter 1 and quarter 2 for a given model, then that can be done directly from Table 4 and Fig. 4.
Here all estimates are relative to FP in sediment type 1 quarter 1, including the significant estimate of FP for sediment type 1 in quarter 2, for e.g.model 3.If an effect is not shown then it is not significantly different from the default (Fig. 4).
For the BD and D integer (count) dependent variables, a negative binomial distribution and log link function were used with the possibility of estimating overdispersion (see also SM S.4).The negative binomial distribution model does not assume a dependency between the mean and variance.In addition, the model enables the inclusion of zero FP as real observations without any addition of constants which is very important given the actual FP ranges (see e.g.Nielsen, 2015).A positive correlation between BD and D per sample was found (see SM S.4 and Fig. S.4), thus, it was necessary to consider the total number of individuals in each sample when parameterising and designing Model 1 for BD analysis (Table 3).However, it also appeared that the higher the number of observations in the data the lower the correlation.For the continuous (non-integer) B and W dependent variables, a Tweedy distribution was used with the reasoning detailed in the SM S.4.
The mixed model design included year as a random effect and other explanatory variables as fixed effects (Table 3).A spatial component with spline (s) between the longitude and latitude was added, allowing the inclusion of spatial and station variability.To assess model performance, overdispersion was checked using residual plots, Q-Q plots were inspected for deviations from homoscedasticity and homogeneous distribution (SM S.5), and variance inflation factors were inspected to check for collinearity.If the assumption of binomial distributed variables did not hold, the overdispersion parameter (k) was estimated (McCullagh and Nelder, 1989; see also SM S.4).In the parameterisation for the R-Packages used here, the variance = m* (1-m/k), where m is the mean.Accordingly, a very large k, as estimated in some instances in the current case, gives the same variance as in the Poisson distribution and accordingly, is no problem; actually on the contrary, see e.g.Table 4 (R Core Team, 2015; Bates et al., 2015; see also SM S.4).All analyses were performed in R (R Core Team, 2015) using the mgcv package (Wood, 2011).

Results of the statistical modelling and multivariate analysis of variance with a mixed GAM model
The main results of the GAMM analyses are shown in Table 4 and Fig. 4 for the interaction (Models 1-4) and main effect models (Models 5-8).All statistical models tested are shown (Table 3) with all dependent and explanatory variables, random and fixed effects, and model input and settings.The number of observations covered, resulting deviance explained, AIC by model, and model performance statistics are shown in Table 4 and SM S.5.
For the biodiversity and density indicators, the models explained a rather high proportion of the variability in the data, i.e. 70-71% and 52-54%, respectively, while for biomass and individual mean weight, a slightly lower proportion was explained, i.e. 38-39% and 42-43%, respectively (Table 4).Both the interaction and main effects models performed well, and the residuals showed no significant trends for any indicators (SM S.5).The main effects models explain a slightly lower proportion of variability in the data and have a slightly higher AIC compared to the interaction models for all indicators, and the number of observations was the same for all models (Table 4).As there are strong significant interaction effects between FP and sediment type and season for all indicators, it is important to take these interactions into account when making a comparative evaluation of the effect of fishing pressure (Table 4; Fig. 4).

Biodiversity (species richness) indicator, BD
FP showed a significant but small negative impact on BD in the main effect model 5 (standardised estimate, std.est., − 0.008, p < 0.007; Table 4; Fig. 4), given the exposure to a relatively high FP interval with abrasions ranging from 0 to 3 events per quarter of the year (Table 1).However, in model 1, there were significant negative interaction effects between FP and habitat type and season of year, with a relatively high estimated negative impact of FP on sediment type 1 (mud; std.est.− 0.023, p < 0.003), and especially in quarter 2 for sediment type 1 (std.est.− 0.172, p < 0.006), compared to the harder bottom sediment types and quarter 3 (Table 4; Fig. 4).These negative impacts are more extensive than the observed significant positive interaction effect between FP and quarter 2 (std.est.0083, p < 0.001).Compared to environmental drivers, especially the bottom current speed (model 1 std.est.0799, p < 0.001), but also minimum salinity (model 1 std.est.0.00015, p < 0.001) had a significant positive impact on benthic invertebrate BD, that is, the higher values the higher BD (Model 1, Table 4; Fig. 4).The minimum temperature showed a small but significant negative impact on BD (model 1, std.est.− 0.00048, p < 0.001).Even though their impacts were in different directions, the environmental drivers have a similar order of magnitude effects on benthic invertebrate biodiversity to that of fishing pressure in the short term, especially, the current speed (Table 4; Fig. 4).

Density (abundance in number of individuals) indicator, D
The main effect of FP on D is significantly negative (Model 1, std.est.− 0.037, p < 0.001), especially with a strong negative impact on sediment type 2 (sand; model 1, std.est.− 0.654, p < 0.001).However, the interaction effects of FP and habitat type showed both strong negative and relatively strong positive impacts on D (Table 4; Fig. 4).There was a strong negative impact of FP on D in sediment type 2 (sand; see above) compared to harder mixed sediments and sediment type 1 with some positive impact of FP on D (mud; model 1, std.est.0.183, p < 0.001).There were no statistically significant interaction effects between FP and season on D which indicates similar impacts of FP on D in different seasons of the year.Compared to the environmental drivers, the bottom current speed did not have a significant impact on D (Fig. 4; Table 4).For the other environmental drivers, there were small, but statistically significant, positive effects of minimum bottom temperature (model 1, std.est.0.01, p < 0.001), salinity (std.est.0.003, p < 0.001), and oxygen (std.est.0.00004, p < 0.001), that is, the higher the minimum values, the higher the density.Even though there were significant impacts on D, the environmental parameters did not have higher influence on D than FP, as was found for BD (Table 4; Fig. 4).In contrast, FP had the highest impact.

Biomass indicator, B
The impact of FP on B was more complex.Similar to BD and D, FP had a significant negative impact on B as the main effect in the interaction model 3 (std.est.− 0.067, p < 0.001; Table 4; Fig. 4) but not in the main effects model 7. Similar to BD and D, there are also highly significant interaction effects for B between FP and sediment type and season of the year which have both strong negative and positive impacts on B (Table 4; Fig. 4).Furthermore, similar to BD, there was a higher estimated negative impact of FP in quarter 2 (model 3, std.est.− 0.667, p < 0.001); however, in contrast to BD, the effect of FP on B showed a significant positive interaction with sediment type 1 (mud; model 3, std.est.0.167, p < 0.001), and especially in quarter 2 (model 3; 0.963, p < 0.001), compared to the harder bottom sediment types and quarter 3 (Table 4; Fig. 4).When compared to the environmental drivers, the bottom current speed had a significantly high impact on B in consistency with that on BD, but for B, it had a negative impact (model 3; std.est.− 7.812, p < 0.001), that is, the higher the bottom current speed, the lower the biomass, which is higher than the FP impact (Fig. 4; Table 4).In accordance with D, there were small, but statistically significant, positive effects for the other environmental drivers for both the main effects (model 7) and the interaction effects models (model 3) of minimum salinity (model 3; std.est.0.004, p < 0.001) and minimum oxygen (model 3; std.est.0.00007, p < 0.001) levels, but with a slightly higher impact for the minimum bottom temperature (model 3; std.est.0.016, p < 0.001); the higher the minimum values, the higher the biomass (Table 4).

Individual Mean Weight indicator, W
FP has a statistically significant positive impact on W, as the main effect in both interaction model 4 (std.est.0.046, p < 0.002) and main effect model 8 (std.est.0.029, p < 0.017) (Table 4; Fig. 4).However, highly significant and strong negative and positive interaction effects were observed between FP and the different habitat types and seasons (Table 4; Fig. 4).There were high estimated negative impacts of FP on W in sediment type 1 (mud; model 4; std.est.− 0.178, p < 0.001), but a positive impact in sediment type 2 (sand; model 4; std.est.0.219,

Table 3
Overview of main tested statistical GAMM models sorted by indicator with BioDiversity (BD), Density (D¼N_ind), Biomass (B) and Individual Mean Weight (W) as dependent variables and fishing pressure (FP) and environmental drivers as explanatory variables, as well as with model settings included.Models 1-4 cover full interaction models with interactions between quarter, sediment type and FP, and models 5-8 cover main effects only models for the different indicators.In all models, quarter and all hydrographical factors are fixed effects (factors).The logμ i is the link function with biodiversity BD as μ i for models 1 & 5, density D (=N_ind) as μ i for model 2 & 6, Biomass B as μ i for model 3 & 7, and Individual mean weight W as μ i for model 4 & 8.The α, β, γ, δ, ϑ, ρ, τ, ω, ∇, Δ are respectively the intercept, the covariates for number of individuals (D=N_ind), minimum temperature (t_min), minimum salinity (s_min), minimum oxygen concentration (o_min), maximum current speed (u_max), fishing pressure (FP) and the factors quarter of year (quarter), habitat type (sed_type), and year.The Δ(year is independent and there is used a smoother on year s(year,bs="re") and on the spatial component s(lon,lat,k = 75) with k = 75.In the models for BD the number of individuals (D=N_ind) are included because of dependency of BD to D (see explanations in text and SM S.4).For the integer dependent variables BD and D, a negative binomial distribution was used in the models, and for the B and W continuous variables (indicators) a Tweedie distribution was used, i.e.Distribution : N i = Negbin(μ i , φ) where φ is the overdispersion parameter for models 1, 2, 5 & 6, and the Distribution : N i = Tweedie(μ i ) for models 3, 4, 7 & 8, respectively.A spatial component with spline (s) between longitude and latitude was added allowing inclusion of the spatial and station variability, i.e. the variability between observations.There has for all models been used a smoother on year s(year,bs="re") and on the spatial component s(lon,lat,k = 75) with k = 75 as we are not focusing on estimating the spatial, station specific variability and the yearly variability, but still want to allow the models to take this variability into account when estimating the short term differences between explanatory variables.The runs cover Van Veen Grab samples only, and the number of observations analysed for all models are 24565.The analyses included samples with zero fishing pressure at the benthic sampling stations and large mussels M. edulis and A. islandica.

Model Number
Model Definition (GAMM) p < 0.007) compared to course mixed sediments (Fig. 4; Table 4).No statistically significant interaction effects were found between FP and season on W, indicating similar impacts of FP in different seasons of the year on W. Similar to the D and B indicators, there were for W (=B/D) small, but statistically significant, positive impacts observed for the environmental drivers in the form of minimum bottom temperature (model 4; std.est.0.008, p < 0.001) and oxygen content (std.est.0.0003, p < 0.001), that is, the higher the minimum values, the higher the average individual weight (Table 4).However, for minimum salinity, a small, but significant negative impact on W was observed, despite the significant positive impacts of minimum salinity on both D and B. Similar to the D indicator, W was strongly negatively impacted by the maximum bottom current speed (model 4; std.est.− 1.998, p < 0.002), that is, the higher current speed, the lower the individual mean weight.Similar to BD and B, the environmental drivers in the form of maximum current speed had a higher impact on the indicators than FP did (TableIe 4; Fig. 4).
Generally, the main effect models (models 5-8) seem to be oversimplified when looking at the significant interaction effects of fishing impact with other explanatory variables on the different indicators in the interaction effects models (models 1-4; Table 4; Fig. 4).However, from the main effects models, FP still appeared to have a significant negative impact on BD (model 5; std.est.− 0.008, p < 0.007), and a relatively higher significant positive impact on W (model 8; std.est.0.029, p < 0.017), while environmental driversin particular, bottom current speedhad stronger and highly significant positive impacts on BD and negative impacts on B and W than FP (Table 4; Fig. 4).

Different types of impacts of fishery and environmental drivers and their interactions
The current study compared four main indicators of benthic invertebrates and their reactions to the short-term variability of both fishing pressure and environmental drivers.We managed to build a suite of good performance statistical models, using the extensive dataset sampled in 2009-2010.Environmental drivers and fishing pressure had a significant impact on short-term, seasonal, and habitat-specific variability.Hence, both factors and their interactions should be considered when evaluating the effects of fishing on the benthic community in the long term (see also conclusion).This is especially relevant in turbulent marine environments with high natural disturbance and gradients in fishing pressure, such as in the Fehmarn Belt area of the western Baltic Sea.The region is an important fishing area with area -fishing pressures ranging from none to three fishing events per quarter of the year, depending on season and habitat type, and also includes areas with no fishing activity from both short-and long-term perspectives (ICES, 2019a; b; 2020; this study).In addition, the Fehmarn Belt is the main inflow and outflow corridor between the brackish Baltic Sea and the northeast Atlantic and is exposed to relatively high levels and ranges of physical drivers (FEMA, 2013;Snoeijs-Leijonmalm et al., 2017;current study).
On the one hand, fishing had a negative short-term impact on benthic invertebrate biodiversity in the area, especially in muddy sediments and in the second quarter of the year, compared to the coarser sediment habitats.In addition, fishery negatively impacted density, with the strongest effect on soft sandy habitats.In contrast, fishing affected invertebrate density positively, but less pronounced, in soft muddy sediments.Similarly, fisheries negatively affected biomass in the short term, especially in the second quarter of the year, but also showed pronounced positive impacts in muddy sediments and in the second quarter of the year.The individual mean weight was affected in both muddy and sandy habitats, but had negative impacts on mud and positive impacts on sand compared to hard sediment.As for density, there were equal impacts of fishing over the year for the individual mean  weight, that is, there were no seasonal differences across habitats.
On the other hand, we found that environmental drivers had highly significant short-term impacts on all indicators.Even though in the same order of magnitude, the environmental drivers impacted biodiversity, biomass, and individual mean weight more than fishing pressure, especially current speed, while they had less pronounced effects on the density, although significant.Minimum salinity, temperature, and oxygen concentrations showed positive effects on the density, biomass, and biodiversity indicators, (the higher the minima, the higher the indicator's value), except for minimum temperature on biodiversity and minimum salinity on individual mean weight with significant negative impacts.The bottom current speed showed highly significant short-term negative impacts on the density and biomass of the benthic community, and had substantial positive and highly significant effects on its biodiversity.Consequently, the higher the current speed.the higher the biodiversity, and the lower the biomass and individual mean weight.

Environmental drivers of short-term habitat-specific differences in the benthic community
These results are in accordance with the significant short-term seasonal and habitat-specific differences found in the area for the same indicators and environmental drivers (e.g.FEMA, 2013;Zettler et al., 2000Zettler et al., , 2017;;Gogina et al., 2014Gogina et al., , 2017)).It has been found that seasonal dynamics are more pronounced in the shallow waters of the Fehmarn Belt region, for example at 6 m depth west off Fehmarn.Here, the soft bottom is inhabited by an infauna community characterised by a common cockle Cerastoderma, however, in spring, characteristic algae (mostly Halosiphon tomentosa, which disappears in summer and autumn) attract a typical epifauna Gammarus community.Another phenomenon is that winter storms severely disturb the sediment in shallow areas, causing the abundance of bivalves such as Mya arenaria, Cerastoderma edule, and Limecola balthica to decrease in the following spring (FEMA, 2013).Accordingly, biodiversity is expected to increase on soft bottoms in spring which makes the negative fishing impacts found in the current study on biodiversity even more pronounced.
Although in the deeper waters of the Fehmarn Belt, the benthic invertebrate community was relatively constant over time, even in deep waters a certain degree of temporal variability was observed.In the central parts inhabited by the so-called "Arctica community" on muddy sediments, a slightly higher number of species, abundance, and biomass of most species are expected in autumn than in spring (FEMA, 2013).Species dominating the biomass mostly reproduce in summer and autumn and exhibit planktonic larval life stages afterwards.Hence, the lowest recruitments were observed during the winter months (Gogina et al., 2014).We found significantly higher biodiversity and density and lower biomass in spring than in autumn in general.Additionally, the environmental drivers in the form of current speed and fishing pressure had strong opposite effects on biodiversity and biomass, complicating these interactions.
Furthermore, contributing to the seasonality and complexity, seasonal hypoxia events that occur predominantly in the deeper areas of the study region negatively affect the diversity and density of soft-bottom fauna (Zettler et al., 2000).Such events are usually recorded in late summer and autumn, resulting in the retreat of species and low-diversity communities observed during this season.Recolonisation occurs within a relatively short time, leading to the general outcome that long-living species are substituted by short-life span Polychaeta species (Zettler  et al., 2000).
However, for many habitats in the southwestern Baltic Sea, where bivalves dominate the biomass, interannual fluctuations in biomass are not as high as fluctuations in abundance and species richness.For community bioturbation potential (BPc), an indicator of benthic faunal function based on bioturbation that combines information on both abundance and biomass, differences between stations in the deeper habitats of the southwestern Baltic Sea were found to be more pronounced than those between seasons within one station (Gogina et al., 2017).
The seasonal and habitat-specific environmentally driven dynamics documented above in the benthic invertebrate community of the current study area may cause the sensitivity to FP and environmental drivers to be very different and change during different seasons in different habitats for the impact indicators studied.That is, the directly observed main effect value levels (maxima, minima and means) of fishing pressure and environmental drivers may not directly explain the levels and directions of seasonal differences because they depends on strong interactions with the habitat type and other factors as well, and on potential different seasonal habitat-specific sensitivities and threshold levels for impacts according to the indicators investigated.

Habitat-specific short-term effects of fishing pressure
Direct, short-term impacts from fishery on the benthic community have been documented to depend on fishing gear, sediment type, and vulnerability of the organisms (e.g.Jennings and Kaiser, 1998;Hall, 1999;Kaiser et al., 2002;Sciberras et al., 2018).More broadly, long-term impacts and indirect effects on the broader biological community structure and biodiversity have been found, beyond the habitats exposed to bottom gears (e.g.Kaiser et al., 2002Kaiser et al., , 2006;;Bremner et al., 2003;Collie et al., 2000;Worm et al., 2006;Hiddink et al., 2016aHiddink et al., , 2017Hiddink et al., , 2019;;Kaiser, 2019;Rijnsdorp et al., 2020;McConnaughey et al., 2020).However, these studies did not consider the interaction effects or directly compare the co-occurring pressures and impacts of environmental drivers in the benthic environment.Some of these studies also lack robust statistical analyses of quantitative based indicators of benthic fishing pressure impacts.
Our analyses showed highly significant influences and short-term variability in the impacts caused by environmental drivers.Here, the bottom current speed impacts were of the same order of magnitude or larger as the fishing impacts.This is supported by Szostek et al. (2016), who found no significant short-term effects of scallop dredging across a wide gradient of fishing intensities using multivariate covariance and principal component analyses.In contrast, they found that seabed temperature and tidal bed shear stress could have induced substantial short-term effects on most benthic impact indicators.The bed shear stress was defined as a function of the maximum predicted tidal current and the bed friction coefficient, and reflected physical disturbance of the seabed.Generally, it is a strong predictor of sediment type, benthic biomass, and production (Szostek et al., 2016;Wildish and Peer, 1983;van Rijn, 1993;van Dijk et al., 2012;van Denderen et al., 2015).Szostek et al. (2016) concluded that it was not possible to demonstrate the effect of scallop fishing within the spatial limits of the king scallop dredge fishery in the English Channel.However, historical dredge fishing could have already altered the benthic communities.Only resilient species have remained, or fishing disturbance has had no impact over and above natural physical disturbance within the fishery in the natural turbulent channel environment (Szostek et al., 2016).Similarly, short-term demersal trawling effects on benthic biodiversity and biomass were absent in a sandy and previously untrawled area in the eastern Bering Sea compared to impacts from an episodic storm event, as demonstrated by a robust Before-After Control-Impact (BACI) statistical design and multivariate analysis of variance (McConnaughey and Syrjala, 2014).It has also been shown that natural temporal variation in community metrics can exceed the effects of dredge fishing at a dynamic study site in the Irish Sea (Lambert et al., 2017).Here it was found that recovery times following scallop dredging on sand and gravel were longer than for beam trawling or otter trawling, with scallop dredging having the greatest negative impact across all habitat types.Typically, fishing is likely to have a prolonged negative effect on communities and habitats in areas that experience low natural disturbance but reduced effects in more dynamic habitats.Lambert et al. (2017) concluded that natural temporal variation in community metrics overall exceeded the effects of fishing in a highly dynamic study site, suggesting that an acute level of disturbance (fished over six times) would match the level of natural variation.These studies suggest that water turbulence could have an overall greater effect on macrofauna than the bottom trawls in turbulent sea areas, such as in arctic shallow-water sandy habitats.In contrast to short term perturbation, previous studies in the arctic Bering Sea area have shown that long-term chronic effects could be significant in densely fished areas on benthic community biomass, diversity, niche breadth, and mean body size (McConnaughey et al., 2000(McConnaughey et al., , 2005)).Wildish and Peer (1983) demonstrated the impacts of tidal current speed on the production of benthic macrofauna in the lower Bay of Fundy.Additionally, van Denderen et al. (2015) reported similar effects of bottom trawling and natural disturbance on the composition and function of benthic communities across habitats in the North and Irish seas.They examined the effects of trawl and natural (tidal-bed shear stress) disturbances on benthic communities over gradients of commercial bottom trawling efforts using a trait-based approach according to life-history strategies or community function characteristics.They found that trawl and natural disturbances affect benthic communities in similar ways, causing declines in long-living, hard-bodied, and suspension-feeding organisms.Furthermore, they found that the benthic response to trawling appeared to be smaller or absent in areas exposed to high natural disturbances.In contrast, in areas with low bed shear stress, responses to trawling were detected resulting in community compositions comparable to those in areas subject to high natural disturbance, with communities being composed of either small-sized deposit-feeding animals or mobile scavengers and predators (van Denderen et al., 2015).
Hence, the current study and literature show that in the short-term, responses to a gradient in the intensity of fishing disturbance can be influenced by the extent to which benthic communities are preconditioned to disturbance by environmental drivers.In addition, a study by van Denderen et al. (2020) indicated strong interactions between fishing impacts and hydrographical conditions when analysing the effects of bottom trawling and hypoxia on benthic communities.The greatest imbalance between environmental and fishing disturbances is likely in muddy substrates (Collie et al., 2000;Diesing et al., 2013;Rijnsdorp et al., 2020).In the present study, we investigated a wide range of demersal-hauled gear fishing pressures, including areas where no fishery has been conducted on a small scale.The fishing pressure in the Fehmarn Belt around the benthic sampling stations from 2009 to 2010 represents the typical fishing pressure ranges and patterns occurring in recent years in the area according to trawl and seine effort allocation in the covered types of benthic habitats and depth strata, ranging from none to three fishing events per quarter of the year.We conclude from this observation that similar long-term effects, as described, for example, in the English Channel and Irish Sea fishery are possible here.However, when including non-fished areas in our analysis, we found that the benthic communities were impacted in the same order of magnitude or more by environmental drivers compared to fishing pressure, considering yearly and seasonal variability and interactions in influencing factors.Therefore, the benthic community in this turbulent sea area could also be influenced in the long-term by continuous interactions with natural and anthropic fluctuations.For a more detailed discussion of the impacts of fishing compared to environmental disturbance, we refer to SM S.6.

Fishing impacts differ between benthic habitat types and life history traits
In their reviews and meta-analyses of the effects of towed bottom gears on marine benthic communities, Collie et al. (2000) and Kaiser et al. (2006) showed that mortality, overall abundance, biodiversity, and species composition vary with habitat type and benthic species group as a result of fishing impact by exposure to and passage of a trawl, and that mortality also varies according to the type of fishing gear (trawl).The highest impacts were observed on sessile epifauna invertebrate species due to scallop dredging, that is hard bottom habitats, followed by severe impacts from beam trawls in softer sandy habitats, and finally otter trawls in soft muddy benthic habitats (see more details in SM S.6).Similar results were obtained in later studies by Sciberras et al. (2018), Kaiser (2019), Rijnsdorp et al. (2020), andMcConnaughey et al. (2020).Overall, it has been concluded that areas which are fished more than three times a year appear to remain in a permanently different state (Collie et al., 2000).The analyses by Collie et al. (2000) indicated that there is a difference in benthic vulnerability between habitats, with benthic invertebrates located in coarse, poorly consolidated sediments being less adversely affected by trawling than benthic invertebrates living in muddy habitats or in and on gravel.This is supported in later studies by e.g.Hiddink et al. (2017), McConnaughey et al. (2020) and Rijnsdorp et al. (2020), and also in current study, which found generally highest fishing impacts of otter trawls and Danish seiners in soft-bottom habitats.
The impact of the gear on the benthic community has recently been quantified based on biological traits such as longevity in relation to vulnerability of the benthic community to trawling, recovery rates, and the ecological function (Rijnsdorp et al., 2016(Rijnsdorp et al., , 2020;;Hiddink et al., 2017;2019).Several studies have suggested that large, long-lived organisms, such as corals, sponges, and especially erect organisms taller than twenty centimeter, will disappear from areas that are regularly trawled (see SM S.6), and that the benthic community in those areas will become dominated by fast-growing, short-lived species (e.g.Bremner et al., 2003;Tillin et al., 2006;De Juan et al., 2007;Buhl-Mortensen et al., 2013, 2016;Jørgensen et al., 2016).Other studies have shown that in areas with low trawling intensity, the abundance of filter-feeding organisms, attached epifauna, and larger animals was higher, whereas areas with high trawling intensity showed a higher abundance of mobile species, infaunal invertebrates, deposit feeders, and burrowing scavengers (Tillin et al., 2006;De Juan et al., 2007).It is noteworthy from our study that some large long-lived epifaunal organisms in the benthic community survive different fishing pressures in the long term, such as the large mussels M. edulis and A. islandica, where the larger specimens of the latter reach an age of 100 years or more and were also abundant in intensively fished areas in our study.Accordingly, general conclusions on the impacts of fishing in highly natural turbulent environments need to be made with caution.

Conclusions
Our results indicate that the short-term influences of environmental drivers and fishing pressure must be considered simultaneously and their seasonal and habitat-specific differences should be acknowledged when evaluating and using quantitative benthic fishing impact indicators.Using averages across seasons, e.g.yearly averages, in longterm impact studies may not be very informative.These results are based on a large number of observations and repeated sampling of benthic habitat data, including high contrasts in fishing pressure and environmental driver ranges.The evaluations are based on robust comparative integrated multifactor statistical analyses of relevant quantitative indicators (e.g.Hiddink et al., 2020).Significant short-term fishing impacts were found with complex interactions across seasons and habitats, while natural disturbance and variability induced by environmental conditions in this turbulent area were at least as impacting, and had significant effects on the quantitative biodiversity, density, biomass, and individual mean weight indicators in the short term.Such interactions need to be considered when evaluating and identifying areas that are more or less resilient to bottom trawling with regard to fisheries management plans in both the short and long term.

Fig. 2 .
Fig. 2. Fishing pressure by Danish and German vessels (>= 15 m length) fishing with towed gears (trawls, seines, dredges) in the Fehmarn Belt area, expressed as categorized annual Swept Area Ratio categories (SAR, km swept / km 2 ) for the year (a) 2009 and (b) 2010, (c) Quarter 2 covering both 2009-2010, (d) Quarter 3 covering both 2009-2010, and as a monthly average over the long term period 2008-2013 for (e) and (f) from 0 to 1.5 (km swept / km 2 ) or more, where (f) gives the benthic sampling station specific subset of the grid cells with long term (2008-2013) fishing pressure.The Fehmarn Belt benthic sampling stations are included in the map as well ('+' symbols).

Fig. 3 .
Fig. 3. Monthly averaged current speed at 17 m, representative for the deeper currents in the Western Baltic region, where the mixed layer depth is typically 10-15 m off coast, for (a) February 2009 (b) December 2009 (c) February 2010 and (d) October 2010.( Source: Data set BALTICSEA_REANALYSIS_PHY_003_011 from Copernicus Marine Service, 2021, marine.copernicus.eu).

Fig. 4 .
Fig. 4. Summary plots of statistical model estimates with std.errors where only the significant results and variables of significant interactions are included.The upper panel (models 1-4) show the interaction effects models, the center panel the main effects models (models 5-8), and the lower panel show results from the interaction models 1-4 with standardised estimates (Std.Est.).The darker colored groups (blue) include the FP variable, while the lighter (yellow) indicate the other variables.

Table 4
Summary of results from Models 1-8 with statistical analyses of the biodiversity (BD, Model 1 & 5), density (D, Model 2 & 6), biomass (B, Model 3 & 7) and the individual mean weight (W, Model 5 & 8) indicators as dependent variables according to explanatory variables including fishing pressure (FP), environmental drivers, habitat type, year and quarter.The table gives family and link function of model, parametric coefficients, relative model parameter estimates and standardized parameter estimates for the FP and environmental drivers, certainty of estimates and significance level.Only the significant results and variables of significant interactions are included, i.e.only estimates and parameter values are given for the significant parameters and effects in the final reduced models (see also SM S.4).Here the final reduced models are given and nonsignificant parameters have been reduced away (empty table cells).The model performance with respect to deviance (the proportion of the variability) in the data explained by the model and the model AIC are given.The runs cover Van Veen Grab samples only, and the number of observations analysed for all models are 24565.The analyses included samples with zero fishing pressure at the benthic sampling stations and large mussels M. edulis and A. islandica.