Trawling effects on biogeochemical processes are mediated by fauna in high-energy biogenic-reef-inhabited coastal sediments

Dynamic, sandy environments are generally less vulnerable to mechanical stress compared to silty, low-energy habitats. Biogenic reef communities, however, may provide an exception to this. This study explores the physical, biological, and biogeochemical effects of bottom trawl fishing on a coastal ecosystem dominated by the tube-building polychaete Lanice conchilega. Two specific gear types, both used to exploit North Sea sole (Solea solea), were compared: electric pulse trawls and tickler chain-rigged beam trawls. We detected a ∼1 cm bathymetric deepening after trawling associated with significant losses in benthic chlorophyll a caused by both fishing gears. Tickler chain trawls significantly reduced sediment oxygen consumption (57 %), total organic matter mineralization (56 %), denitrification (61 %), nitrification (60 %), and total benthos densities (52 %), while pulse trawls had no statistically significant impact on these parameters. Before trawling, significant relationships could be found between L. conchilega and very fine sand fractions, oxygen and nitrate fluxes, macrobenthos densities, and species richness; however, the trawl disturbances from both gears disrupted these connections. Our results suggest a reduced average effect for pulse compared to tickler chain beam trawls for several ecological and biogeochemical characteristics, though their impact was still significant for L. conchilega and associated species. This study also suggests that faunal-mediated ecosystem functions in habitats dominated by L. conchilega may be sensitive to relatively shallow sediment penetration from trawl gears and should be considered when assessing habitat vulnerability.

1 Introduction

Marine sediment ecosystems are important drivers for benthic pelagic coupling (Mermillod-Blondin and Rosenberg, 2006) and carbon mineralization (Seiter et al., 2005). Biogenic reef habitats in particular, such as those created by the tube-building polychaete Lanice conchilega are able to capture substantial amounts of organic carbon compared to adjacent sediments (De Smet et al., 2016a), resulting in dense and biologically diverse benthic communities (Van Hoey et al., 2008). This species may be found in various North Sea habitats but generally prefers shallow muddy and fine sands where it may exhibit densities of up to 1000 individuals per square meter (Van Hoey et al., 2008). Bottom trawling activities are a common occurrence in North Sea habitats dominated by L. conchilega (Rabaut, 2009) and can present a considerable threat to most biogenic reef ecosystems (Kaiser et al., 2002; Fariñas-Franco et al., 2018). Biological resilience to bottom trawl impacts have, however, been linked to both L. conchilega reefs (Rabaut et al., 2008) and the dynamic sandy habitats they are often found in (Kaiser et al., 2006; Sciberras et al., 2016; Rijnsdorp et al., 2020a). Consequences of trawling on more holistic ecosystem processes such as benthic pelagic coupling, however, have not been thoroughly examined in these systems.

Bottom fishing affects benthic carbon cycling by displacing bottom-dwelling organisms (Hiddink et al., 2017; Sciberras et al., 2018), altering sediment structures leading to the impairment of biogeochemical processes (Depestele et al., 2018; Ferguson et al., 2020), and resuspending fresh organic material from the sediment surface, which can result in reduced mineralization of organic carbon in seabed sediments (Tiano et al., 2019; Paradis et al., 2020; Bradshaw et al., 2021; Morys et al., 2021). Conversely, trawling may also increase benthic carbon mineralization through depositional effects or if refractory organic fractions are mixed with and reactivated by fresh organic matter (van de Velde et al., 2018). Trawl-induced mortality to benthic organisms can simplify benthic communities and can lead to a loss in functional diversity (Tillin et al., 2006; Tiano et al., 2020). The removal (significant or partial) of an abundant ecosystem engineer such as L. conchilega and associated fauna would likely have major impacts on biogeochemical and community dynamics (Rabaut et al., 2008; Van Hoey et al., 2008; De Smet et al., 2016a; Foshtomi et al., 2018).

The severity of bottom trawl impact depends on the size and type of fishing gear used and speed at which it drags over the seabed (Eigaard et al., 2016; Depestele et al., 2016, 2018; Rijnsdorp et al., 2020a, b). In the southern North Sea, beam trawls rigged with tickler chains mechanically stimulate flatfish up from the seafloor, facilitating their capture. An alternative fishing gear for flatfish called “pulse trawls” has been used to expose the seabed to electrical fields, causing the immobilization of target fish for oncoming bottom nets (Soetaert et al., 2015b). Previous comparisons conducted in the Frisian Front, a relatively homogenous habitat known for high abundances and diversity of benthic organisms (Dauwe et al., 1998; Dewick et al., 2002), have found a reduced impact of pulse compared to tickler-rigged trawls on physical (Depestele et al., 2018) and biogeochemical parameters (Tiano et al., 2019), but impacts on benthic communities were similar for both gears (Tiano et al., 2020). Dynamic nearshore habitats, in contrast, may provide a more challenging environment for detecting significant anthropogenic alterations.

Regarding sensitivity to direct impacts and recovery potential, coarse-sediment ecosystems characterized by high levels of natural disturbance typically display more resistance and resilience to bottom fishing compared to muddy and more hydrodynamically static ecosystems (Kaiser et al., 2002; Bolam et al., 2014; Sciberras et al., 2016; Rijnsdorp et al., 2020a). Macrobenthos found in sandy, well-sorted sediments tend to be less sensitive to trawl disturbance than those in muddy, poorly sorted habitats (Bolam et al., 2014; van Denderen et al., 2015). Communities dominated by three-dimensional biological features, such biogenic reefs, seem to be an exception to this and can be vulnerable to trawl disturbance regardless of hydrodynamic energy level or substrate type (Kaiser et al., 2002; Grabowski et al., 2014).

The current study uses a combination of water column, sonar, sediment profile imagery, box core sampling, sediment incubation, and porewater nutrient extraction techniques to provide a comprehensive assessment of bottom trawling on ecosystem functioning. This integrated data collection allows us to explore the links between pure mechanical trawling effects, disturbance-induced changes to faunal-mediated biogeochemical functions, and their collective impact on benthic pelagic coupling. Our main objective was to measure the acute ecosystem effects of electric pulse and tickler chain beam trawls in a high-energy coastal system. Furthermore, the high abundances of Lanice conchilega in the region allowed us to quantify the effect of this species on ecological and biogeochemical parameters and how this may change due to fishing disturbance.

Figure 1 Geographic location of the study site (a). Visual descriptions for experimental plots in Pulse-ON (A–C), Pulse-OFF (D–F), and Tickler (G–I) trawl treatments, untrawled adjacent stations (AD1, AD2, AD3), sediment profile imagery (SPI) locations, and box core sampling locations (b). Experimental plots were sectioned into northwest, central, and southeast stations.

2 Methods

2.1 Study area

This study took place in shallow subtidal waters (8–10 m depth), 10 km away from the Netherlands' coast close to the Dutch and Belgian EEZ (exclusive economic zone) border (Fig. 1). This region, known as the “Vlakte van de Raan”, is an area known for its subtidal sandbanks which create a mosaic of habitat areas of varying ecological value (Degraer and Hostens, 2016; Pecceu et al., 2016). The experimental areas were characterized by “fine sand” with an average median grain size of 200.5 µm and an average chlorophyll a content (0–10 cm depth) of 13.1 µg g⁻¹ (Fig. S1 in the Supplement). The exact location for this study within the Vlakte van de Raan (Fig. 1) was chosen due to its relatively low levels of bottom trawl disturbance (no flatfish or shrimp trawl disturbance in 1 year as verified with vessel monitoring data) and the high macrobenthic densities and biomass (Craeymeersch et al., 2006).

2.2 Experimental trawling design

Three experimental treatments were conducted: (1) tickler chain beam trawling (hereafter referred to as “Tickler” treatment or tickler trawls) and (2) pulse trawling with and (3) without electricity turned on (hereafter referred to as “Pulse-ON” and “Pulse-OFF” treatment). Each treatment was repeated three times in spatially separated areas (Pulse-ON = A, B, C; Pulse-OFF = D, E, F, Tickler = G, H, I), resulting in a total of nine experimental plots of 50 by 300 m (15 000 m² per plot; Fig. 1). The locations and orientation of the treatment plots were chosen based on their position within a bathymetric depression. Plots and untrawled adjacent sites were designated to represent northwest, central, and southeast areas within this depression and were spatially separated by at least 100 m to minimize effects of trawl-induced sediment resuspension and deposition. Individual plots were also subdivided by three different sampling stations representing the northwest, central, and southeast locations within the plot. Three additional untrawled stations located adjacent to the treatment plots (AD1, AD2, AD3) were chosen within the experimental area (Fig. 1). Samples taken from AD stations deviated within 80 m of a single sampling point.

Pulse trawls (Delmaco; Soetaert et al., 2015b) were deployed from commercial fishing vessel (FV) TH6 Johanna Cornelia (24 m) and tickler trawls from FV YE238 Geertruida (22 m). The vessels were relatively small trawlers (<221 kW engine) and allowed to fish within 12 nmi of the coast. Both pulse trawl (∼1.1 tonne in air) and tickler trawl (∼2.4 tonne in air) gears measured 4 m in width and had one trawl shoe at either side of its beam. The fishing vessels made a total of six trawl passes over each experimental plot at a speed of ∼5 knots (pulse) and ∼6 knots (tickler). Experimental fishing with the pulse trawl was performed on 12 June 2018, while the tickler trawl was used on 15 of June 2018 (Table S1 in the Supplement). All experimental equipment was deployed from the research vessel (RV) Simon Stevin.

2.3 Water column turbidity

To measure turbidity during trawling events, the RV deployed an SBE 25plus Sealogger (CTD) with two optical backscatter point sensors (OBSs) while being suspended continuously at 5 m depth during experimental trawling activity. For each treatment plot measured (A, B, C, D, E, G, and H), the RV strategically positioned itself against the current and the ongoing trawling ∼80 m adjacent to the corresponding plots (due to logistic constraints, plots F and I were not measured).

2.4 Bathymetry and backscatter

Bathymetry and backscatter data were collected using a Kongsberg EM2040 single head multi-beam echo sounder (MBES) mounted on the RV (Depestele et al., 2018). All data were collected with an MBES frequency of 300 kHz, a speed of 8 knots and at an orientation parallel to the longest side of the experimental plots to avoid backscatter bias from the survey line directions (Lurton et al., 2018). Bathymetrical and backscatter data were collected before (T0) and after (T1) trawling (Table S1). The RV followed the FV at an approximate distance of 300 m during one single haul in each of the experimental plots to register the bathymetrical data. After six passes in each experimental plot, the RV collected MBES backscatter data within 12 h.

A high-resolution (0.5×0.5 m) digital elevation model of the seabed was created from the MBES data using Qimera (QPS). Trawl tracks were located from the GPS-derived position of the vessel path and visual detection using backscatter and bathymetrical information. Transects of 12 m width were positioned every 20 m along the trawl track in GIS software (N=74 for Pulse-ON, N=100 for Pulse-OFF, and N=94 for Tickler; Fig. S2). For the transects with visible trawl marks, the mean water depth was calculated for five locations inside and outside the trawl tracks. The water depth was estimated at these same locations at T0. The difference between mean water depths inside and outside the track was used to assess if the trawl caused bathymetrical changes in the track.

A 1×1 m backscatter mosaic was created in FMGeocoder Toolbox (QPS FMGT). Backscatter data were sampled by using transects for comparison of T0 and T1 and using all backscatter data in each experimental plot. Values outside 1.5 times the interquartile range were removed as outliers (1.5 % of all values).

2.5 Sediment profile imagery

Sediment profile imagery (SPI) was used to characterize surface and subsurface sediment features using a Nikon D7100 digital camera (Ocean Imaging Systems, North Falmouth, MA, USA). The SPI is a device (lowered onto the sediment from the RV) that captures images of the sediment–water interface in a vertical cross-section and can provide information on benthic fauna, physical disturbance, and chemical features (e.g., redox depth) within the upper sediment (Rosenburg et al., 2003). Three replicate images were taken at each sampling station (three stations per treatment plot; Fig. 1). A total of 180 SPI photos were captured before (T0) and after (T1) trawling in the treatment plots (A–I) and also at the beginning (T0) and end (T1) of the experimental period for untrawled adjacent stations (AD1–3; Table S1). Information for prism penetration depth (depth at which the camera penetrates the sediment) and visual observations of fauna, sediment characteristics (appearance of fine or coarse sediment, sand ripples, etc.), and oxidation, expressed as percent anoxic sediment (assessed by sediment color (black = anoxic)), were obtained using SpiArcBase software (Romero-Ramires et al., 2013). Individual SPI images can be found in Sect. S4.

2.6 Benthic samples

Samples for sediment, macrofauna, porewater nutrients, and biogeochemical incubations were collected using cylindrical NIOZ (Royal Netherlands Institute for Sea Research) box corers (50 cm height, 32 cm diameter; Netherlands Institute for Sea Research, Texel). For time steps T0 and T1 for the fishing treatments, three replicate box cores were taken inside individual plots (northwest, central, and southeast locations). Samples were collected at untrawled adjacent stations AD1 and AD2 at the beginning (T0) and end (T1) of the experimental period; however, AD3 samples were only collected at the T1 time step. This amounted to 69 total box cores collected and processed for this study. During experimental fishing activity, box cores in trawled sediments were collected within 0.5 h after six consecutive trawling passes were made in each treatment plot. Detailed sampling information can be found in Table S1.

2.6.1 Sediment parameters

Two out of the three box cores were subsampled for sediment parameters from each trawl treatment plot and untrawled adjacent station. Sediment samples were subdivided at 0–0.5, 0.5–1, 1–1.5, 1.5–2, 2–2.5, 2.5–3, 3–4, 4–5, 5–6, 6–8, and 8–10 cm slices taken from 3.5 cm inner-diameter sub-cores. All samples were stored in a −20 ^∘C freezer before freeze-drying and sieving (1 mm) for particle analysis. Sediment grain sizes were obtained via laser diffraction using a Malvern Mastersizer 2000 (McCave et al., 1986). Chlorophyll a (chl a) pigments were extracted using acetone and analyzed with UV spectrophotometry (Ritchie, 2006).

2.6.2 Porewater nutrients

Porewater was collected in the same box cores used for sediment parameters using 10 cm diameter sub-cores with vertical sampling ports. Porewater was extracted at 0, 1, 2, 3, 5, 7, and 10 cm depth using rhizon samplers (Rhizosphere Research Products; Seeberg-Elverfeldt et al., 2005; Dickens et al., 2007; Shobolt, 2010). Porewater samples were added to 10 mL polystyrene vials and stored in −20 ^∘C prior to analysis. To acquire nutrients, samples were thawed and analyzed using a SEAL QuAAtro segmented flow analyzer (Jodo et al., 1992) to determine concentrations of ammonium (NH${}_{\mathrm{4}}^{+})$, nitrite (NO${}_{\mathrm{2}}^{-})$, nitrate (NO${}_{\mathrm{3}}^{-})$, phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-})$, and silicate (Si(OH)₄).

2.6.3 Benthic biogeochemical fluxes

Incubation cores (15 cm inner diameter by 30 cm height) were collected from all box core samples (three times per treatment plot and reference station) with overlying water. These were subsequently placed in a water bath filled with surface seawater on board the RV and aerated ∼6 h to acclimatize. Magnetic stirrers, fixed to the lid of incubation cores, ensured homogenized water column conditions. After the acclimatization period, cores were sealed from air contact for 4 h to measure rates of sediment community oxygen consumption (SCOC) in the dark. During this period, oxygen concentrations were monitored every 30 s using optode sensors placed in the overlying water of the incubation cores (FireSting-O₂, PyroScience). At the end of the O₂ incubations, cores were re-opened and aerated to remain at above 90 % O₂ saturation for the ongoing nutrient incubations. Water samples of 10 mL were collected at the overlying water of incubation cores at 0, 4, 8, and 16 h time intervals to obtain flux information on NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{2}}^{-}$, NO${}_{\mathrm{3}}^{-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and Si(OH)₄. Water samples for nutrient fluxes were stored and analyzed using the same methods described for porewater nutrients. Flux estimates were obtained by fitting linear regressions on nutrient concentration changes over time and multiplying the regression coefficient with the height of the overlying water (Tiano et al., 2019).

2.6.4 Mass budget modeling

Rates for total mineralization of organic matter, denitrification, nitrification, and the percentages of nitrogen and phosphorous removal were estimated using the measured sediment–water exchange fluxes of O₂, NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{2}}^{-}$, NO${}_{\mathrm{3}}^{-}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ calculated from the core incubations. These flux values were used to create an integrated mass budget of the solutes within the sediment column (Soetaert et al., 2001; Braeckman et al., 2010). A linear inverse model is implemented that estimates the rates of nitrification, denitrification, total mineralization, and nutrient removal using the R package “limsolve” (Soetaert et al., 2009). More detailed methodology can be found in Sect. S5.

2.6.5 Macrofauna collection

After the flux measurements, incubation cores were rinsed through a 1 mm sieve to collect the benthic macrofauna inside them. Faunal samples were preserved in 8 % formalin seawater on RV Sorting, and the identification of preserved macrofauna samples took place at the Flanders Institute for Agricultural and Fisheries Research (ILVO) to obtain measurements for taxon densities and biomass (blotted wet weights) following the accredited procedure BELAC ISO 17025 norm (ILVO-DIER-ANIMALAB; certificate no: BELAC T-315).

2.7 Statistical analysis

All figures and statistical analyses were carried out using R (R Core Team, 2021).

2.7.1 Water column

Differences between water column turbidity during Pulse-ON, Pulse-OFF, and Tickler treatments were assessed with one way ANOVA and Tukey HSD pairwise tests. Assumptions of normality and homogeneity of variances were verified using Shapiro–Wilk and Levene's tests. If parametric assumptions were violated, data underwent a natural log-transformation prior to analysis.

2.7.2 MBES information

Bathymetrical differences between water depths inside and outside the track locations ${\mathrm{\Delta }}_{i,T,G}$ were calculated for each transect i at each time interval T (T0 = before and T1 = after trawling) for each trawl gear type G following Eq. (1). A paired Wilcoxon signed-rank test was used to test if delta values were significantly different between T1 and T0 for each gear type. Significant positive delta values indicated that the water depth inside the trawl track was higher than outside the trawl track.

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }}_{i,T,G}={\left({\text{water depth}}_{\text{inside}}-{\text{water depth}}_{\text{outside track location}}\right)}_{i,T,G}\end{array}$$

The indention I of each transect i was calculated by subtracting the delta value at T0 from the delta value at T1 following Eq. (2).

$$\begin{array}{}\text{(2)}& I_{i,G}={\left({\mathrm{\Delta }}_{T\mathrm{1}}-{\mathrm{\Delta }}_{T\mathrm{0}}\right)}_{i,G}\end{array}$$

A resulting positive value indicated a deepening of the trawl track. The indention of the trawl was compared between trawl types using a Kruskal–Wallis rank sum test.

Statistical differences in backscatter values were compared for the combination of time and trawl types using experimental plots as replicates in a Kruskal–Wallis rank sum test. The experimental plots were also compared using the Kruskal–Wallis rank sum test. A Dunn test with Benjamini–Hochberg method was applied for pairwise comparisons between sites.

2.7.3 SPI and benthic sample analysis

Linear mixed effects models (LMMs) were used to investigate significant differences before and after fishing (T0–T1) for SPI measurements, sediment parameters, porewater nutrient concentrations, biogeochemical fluxes, mass budget model results, and ecological characteristics (macrobenthos densities, biomass, species richness) using the lmer function in the R package “lme4” (Bates et al., 2015). For a given variable and treatment (example: oxygen flux and tickler treatment), a full model was specified using the “time step” (T0–T1) as a fixed effect variable, “temperature” as a co-variate, and “station” as a random effect variable. A reduced model was created to consider only the random effect variable (station). The full and reduced models were tested against each other using a partial F test. This approach was taken to assess the effect of fishing disturbance (time step = before and after statistical comparison) with respect to temperature while minimizing spatial autocorrelation between stations.

2.7.4 Impact of Lanice conchilega

Robust regression analysis, using the R package “robustbase” (Maechler et al., 2020), was conducted to investigate the effect of L. conchilega densities and biomass on physical, biogeochemical, and faunal parameters in trawled and untrawled samples. This alternative to least squares techniques is an iterative analysis which down-weights the influence of outliers on regression coefficients and is less affected by violations of linear regression assumptions (Koller and Stahel, 2011). To reduce possible bias from the spatial differences in L. conchilega densities/biomass between treatment plots, all pre-trawled and post-trawled data were pooled together for this analysis.

Figure 2 Water column turbidity during the time of trawling for Pulse-ON plots (a–c), Pulse-OFF plot (d), and Tickler plots (e–f).

3 Results

3.1 Water column conditions during trawling

During the experimental period, the water column exhibited mixed thermal conditions (no thermocline) with temperatures ranging between 15.8 and 17.0 ^∘C. Turbidity measured by the OBS sensors detected significantly higher amounts of water column SPM (suspended particulate matter) concentrations (mg L⁻¹) during tickler trawling (mean ± standard deviation = 62.5 ± 17.1, p<0.001) compared to Pulse-ON (19.6 ± 5.3) and Pulse-OFF treatments (13.3 ± 5.8; Fig. 2). Data from a “Meetnet Vlaamse Banken” buoy (located ∼7 km southwest of the experimental location) showed current speeds ranging from 0.08 to 1.21 m s⁻¹ during pulse trawling and 0.23 to 0.85 m s⁻¹ during tickler trawling. On the day before both pulse trawl treatments, water currents reached a maximum speed of 1.23 m s⁻¹ compared to 1.5 m s⁻¹ the day before the tickler treatment (Fig. S3).

3.2 Trawl effects on bathymetry

Delta values were significantly different before and after trawling (T0–T1) for each trawl type (Pulse-ON: p<0.05; Pulse-OFF: p<0.001; Tickler: p<0.001; Fig. S4). The track locations were deepened, on average, by 1 cm after trawling (Pulse-ON: 0.99 ± 3.92; Pulse-OFF: 1.17 ± 1.84; Tickler: 1.11 ± 2.74). The indention of the trawl track was not significantly different between trawl treatments.

3.3 Trawl effects on backscatter

The mean backscatter values of the track locations in all treatments exhibited reduced values after one trawl pass (ΔT0–T1: Pulse-ON = −0.20; Pulse-OFF = −0.23; Tickler = −0.30), though there was no statistical difference between the three trawl types (Table S2). Comparison of mean backscatter values over entire plots (after six trawl passes) showed decibel values ranging between −26.64 and −25.40 (Table S2). The mean backscatter strength did not differ by more than 2 dB between the treatments at T0 and T1 and was not significantly different between any of the time–treatment interactions (Fig. S5).

Figure 3 Examples of sediment profile imagery from the study showing a dense Lanice conchilega tubeworm reef (a), bioirrigation-induced sediment oxidation (b), and evidence of trawl disturbance: flattened tubeworms and sediment (c).

Table 1 Results for sediment and biogeochemical parameters between T0 and T1 for Tickler, Pulse-ON, and Pulse-OFF treatments (mean ± SD). Nitrogen (N) and phosphorus (P) removal are presented as the percentage of inorganic N or P removed from the sediments after being formed from organic matter mineralization.

Bold signifies that T1 is significantly different compared to T0. ${}^{*}p<\mathrm{0.05}$; ${}^{**}p<\mathrm{0.01}$; ${}^{***}p<\mathrm{0.001}$.

3.4 Sediment profile image analysis

L. conchilega tubeworms were identified in 72 % of images taken and were often associated with heterogeneous patterns of oxidized sediment (Fig. 3b). Mechanical disturbance was evident in certain T1 images showing “flattened” tubes and/or sediment (Fig. 3c). The percentage of anoxic sediment was significantly lower in T1 images from the tickler treatment compared to T0 images (p<0.01) but did not differ significantly from T0 to T1 in other trawl treatments (Table 1). Anoxic (dark) sediment was significantly higher in images from the untrawled AD1 taken at the end of the experimental period (T1) compared to the beginning (T0; p<0.01). The sediment penetration depth of the SPI prism did not differ significantly between T0 and T1 for any treatment (Table 1).

Table 2 Robust regression coefficients (slope, p value, multiple R squared) showing the relationship of Lanice conchilega densities and biomass versus sedimentary faunal and biogeochemical flux parameters in trawled and untrawled sediments.

Bolded values show significant relationships.

3.5 Trawl effects on sediment parameters

Very fine sand fractions (62.5–125 µm) and chl a content in the sediment column (0–10 cm) were significantly lower in Tickler T1 compared to T0 samples (p<0.01; p<0.001; Table 2). For the Pulse-ON treatment, T1 silt fractions and chl a content were significantly lower than before trawling (p<0.01; Table 2). No significant differences in sediment parameters between T0 and T1 were found in the Pulse-OFF treatment. Particle size distributions of surface sediments (0–1 cm depth) between T0 and T1 for the Tickler and Pulse-ON treatments were characterized by a shift towards larger grain size fractions after trawling (Fig. 4).

Figure 4 Particle size distributions for the surface sediments (top 1 cm) in experimentally trawled samples. Solid lines represent the distribution of sediment size fractions in Pulse-ON, Pulse-OFF, and Tickler plots before trawling (T0), and dashed lines represent distributions after trawl disturbances (T1).

3.6 Trawl effects on biogeochemical fluxes and nutrient profiles

Average SCOC rates (mmol O₂ m⁻² d⁻¹) were significantly reduced between T0 (172±109) and T1 (74±40) for the Tickler treatment (p<0.001; Fig. 6a; Table 1). No other consistent or significant patterns were found for nutrient fluxes between T0 and T1 fishing treatments. Average and absolute values for nutrient fluxes can be found in Tables 1 and S5.

Average porewater nitrate concentrations within the sediment column (0–10 cm) were significantly higher after trawling in Tickler and Pulse-OFF treatments compared to pre-trawled sediments (Tickler, p<0.01; Pulse-OFF, p<0.01; Table 1). Ammonium concentrations in the sediment significantly increased from T0 to T1 (p<0.01) after tickler trawling (Table 1). Mean phosphate concentrations were elevated after trawling for all experimental treatments, but this was significant only in Pulse-ON (p<0.01) and Pulse-OFF (p<0.001; Table 1).

3.7 Mass budget calculations

The Tickler treatment exhibited significantly lower values for total mineralization (p<0.001), denitrification (p<0.05), and nitrification (p<0.05) in T1 compared to T0 samples (Table 1). For all fishing treatments, the percentage of nitrogen removed after mineralization ranged between 80 % and 95 %, while sedimentary phosphorus removal varied between 91 % and 99 %; however, no changes between T0 and T1 were statistically significant for individual fishing treatments.

3.8 Trawl effects on macrofauna density and macrobenthos composition

A total of 69 taxonomical groups were identified in this study, with L. conchilega dominating benthic macrobenthos densities (54 %). The bivalve Kurtiella bidentata and polychaetes Magelona spp. and Eumida sanguinea made up 13 %, 4 %, and 3 % of total macrofauna density respectively. Though they only contributed 1 % of total macrobenthos density, the razor clam Ensis leei represented 45 % of total biomass collected followed by L. conchilega (22 %), Spisula subtruncata (bivalve, 14 %), and Echinocardium cordatum (heart urchins, 4 %).

Figure 5 Oxygen consumption (a) and macrobenthos densities (b) measured from treatment plots before (T0) and after (T1) fishing. The y axes are on a log scale.

Total macrofaunal densities (individuals per square meter) decreased from T0 (Tickler: 11 782±7360; Pulse-ON: 9783±14 072; Pulse-OFF: 6241±5120) to T1 (Tickler: 4949±3428; Pulse-ON: 3336±3529; Pulse-OFF: 5483±3623) but were only statistically significant after the Tickler treatment (p<0.05; Fig. 5b). Reference samples varied greatly but held the highest average densities from all the treatments (T0 = 19 174±17 493; T1 = 22 993±17 945; Fig. 5b).

L. conchilega maintained the highest average densities (individuals per square meter) in the reference stations (T0 = 10 703±9943; T1 = 14 707±11 336). L. conchilega densities declined between T0 and T1 for Pulse-ON (T0 = 5801 ± 9352; T1 = 1587 ± 3022), Pulse-OFF (T0 = 3081 ± 3711; T1 = 2251 ± 3547), and Tickler treatments (T0 = 5880 ± 3838; T1 = 2172 ± 1913) but was only significant for the Pulse-ON treatment (p<0.05). Amongst other abundant taxa, Abra alba (bivalve) density declined significantly after tickler trawling (T0 = 418 ± 632; T1 = 231 ± 420; p<0.05), and E. sanguinea was significantly lower after the Pulse-ON treatment (T0 = 555 ± 1315; T1 = 93 ± 195; p<0.001).

3.9 Effects of Lanice conchilega in trawled and untrawled sediments

L. conchilega was found in 68 out of 69 (98.5 %) benthic sediment samples collected in this study. In pre-trawled sediments, L. conchilega density explained 85 % of variation in the abundances of the other taxa (p<0.001), while biomass explained 38 % of variation in species richness (p<0.001; Fig. 6, Table 2). Tubeworm densities in pre-trawled samples were linked with smaller sediment grain sizes, accounting for 33 % of variability in the percentage of very fine sand (p<0.05), while biomass explained 18 % of variability in silt content (p>0.01) (Table 2). Densities of L. conchilega significantly predicted 87 % of SCOC in undisturbed sediments (p<0.001). Biomass of L. conchilega (compared to density) was a slightly better predictor of SCOC and explained 93 % of variation in O₂ fluxes (p<0.001), as well as 21 % of nitrate effluxes (p<0.01) in untrawled sediments (Fig. 6; Table 2). No significant relationships were found after trawling (Fig. 6, Table 2). Due to a lack of sufficient data, regressions between T0 and T1 within individual fishing treatments (Pulse-ON, Pulse-OFF, Tickler) did not yield consistent results.

Figure 6 Robust regression analysis of sediment community oxygen consumption (SCOC; mmol m⁻² d⁻¹) versus Lanice conchilega biomass (grams of wet weight per square meter) (a, b) and densities of other taxa versus L. conchilega densities (individuals per square meter) (c, d) in pre-trawled and trawled sediments. Note x-axis scales are different.

4 Discussion

The data gathered in this study exhibited a large amount of spatial variability, reflecting the dynamic nature of the Vlakte van de Raan ecosystem (Degraer and Hostens, 2016). Despite the challenge of detecting disturbances in such a variable system, we were able to provide a comprehensive assessment of acute bottom trawling with our 3×3 plot design. The dominance of L. conchilega in this habitat provided an additional opportunity to explore the effects of tubeworm abundance on environmental parameters in disturbed and undisturbed conditions.

This study follows up on similar research conducted in the Frisian Front region of the North Sea (Depestele et al., 2016; Tiano et al., 2019, 2020) with each study (including the current one) collecting information in the same season (early June) of its respective year. The offshore Frisian Front ecosystem (∼60 km from the Dutch Coast), however, experiences much lower levels of tidal disturbance, exhibits smaller grain sizes, and holds completely different macrofaunal assemblages and biogeochemical characteristics compared to the Vlakte van de Raan (Van Raaphorst et al., 1992; Dewicke et al., 2002; Degraer and Hostens, 2016; Depestele et al., 2018). By comparing the findings from these two contrasting areas, we can better understand how different ecosystems are affected by bottom disturbance. To our knowledge, this is the first coastal habitat study that simultaneously investigates acute beam trawl impacts on biological (Rabaut et al., 2008), physical (Depestele et al., 2016), and biogeochemical dynamics.

4.1 SPM during experimental trawling

The water column was 73 %–78 % more turbid during tickler trawling compared to both pulse trawl treatments. As tickler trawls tow at faster speeds compared to pulse trawls (Poos et al., 2020), they are predicted to mobilize more sediment (Rijnsdorp et al., 2021). However, due to the higher current speeds evident in the day prior to tickler trawling, we cannot be certain that this increased turbidity was the result of enhanced resuspension caused from bottom trawling alone. It is possible that trawling coupled with current-induced resuspension enhanced SPM measurements during the Tickler treatment. This would corroborate the findings in Depestele et al. (2016), in which no significant differences in sediment resuspension of pulse and tickler trawling were found in a similar coastal area. In the Frisian Front, where Tiano et al. (2019) reported conspicuous SPM peaks attributed to individual trawl passes, water current speeds were 62 %–88 % lower during fishing compared to the Vlakte van de Raan study. The more ambiguous link between trawling and SPM in our study (no obvious turbidity peak during trawl passes) suggests that the water column effects of trawling are more difficult to detect in areas with higher natural current speeds. These effects may be diluted through current-induced dispersion of resuspended particles and/or higher levels of natural turbidity.

4.2 Changes to the physical environment

A single trawl pass for all three treatments yielded a bathymetrical indentation of ∼1 cm in the sediment, similar to the 0.9 cm found after tickler trawling in coarse sediment off the Dutch coast (Depestele et al., 2016). Our study, however, did not replicate previous findings, which found greater bathymetrical effects for tickler beam trawling in finer North Sea sediments (Depestele et al., 2018). This is likely due to differences in grain size which are smaller in the Frisian Front (aforementioned study) compared to the Vlakte van de Raan (current study). Backscatter data, which produced similar values between T0 and T1 in our study, can approximate sediment characteristics and trawling disturbance (Roche et al., 2018) but do not purely reflect sediment grain size and can be distorted by “scatterers” such as shell debris and biogenic reefs, which may interfere with the acoustic signals (Degraer et al., 2008; Feldens et al., 2018).

Tubeworm reefs in untrawled sediments were positively associated with smaller grain size fractions (Fig. 3; Table 2). This pattern is commonly found with L. conchilega as high tube densities will slow down bottom current velocities leading to the settling of smaller particles (Van Hoey et al., 2008; De Smet et al., 2015; Foshtomi et al., 2018). Trawl-induced coarsening of sediments found in habitats dominated by L. conchilega may occur through both the resuspension and winnowing of finer-grained particles from direct trawling impacts (Palanques et al., 2014; Mengual et al., 2016; Depestele et al., 2018; Tiano et al., 2019) and the alteration of sediment characteristics caused by the rapid extraction of tube structures from the seabed.

The shift in particle sizes and changes in observed anoxic sediment after trawling suggest that the largest sediment alterations were attributed to the Tickler treatment followed by the Pulse-ON treatment, with no significant effects detected after Pulse-OFF trawling (Fig. 4; Table 1). As electricity from pulse trawls does not affect physical sediment characteristics (Tiano et al., 2021), differences between Pulse-ON and Pulse-OFF treatments may have been due to the variable mechanical impacts of pulse trawls (Depestele et al., 2018) or the environmental variability of the site (Degraer and Hostens, 2016).

4.3 Changes to biogeochemical characteristics

We detected significant declines in chl a from surface sediment (0–1 cm) of 63 % and 37 % after tickler and pulse trawling (Pulse-ON) respectively. These percent changes are lower relative to the reductions measured after tickler (83 %) and pulse (43 %) trawling in the siltier Frisian Front sediments (Tiano et al., 2019) even with the 3 to 6 times higher sediment chl a content in the more coastal Vlakte van de Raan system. Sandy dynamic habitats, like the Vlakte van de Raan, experience higher bed shear stress and are generally less vulnerable to trawling impacts compared to muddier more static environments like the Frisian Front (Kaiser et al., 2006; Queirós et al., 2006; Allen and Clarke, 2007; Sciberras et al., 2016; Rijnsdorp et al., 2020a).

Measured SCOC values in the present study were exceptionally high, indicating biogeochemically active sediments (Stratmann et al., 2019). There are several factors driving high SCOC in the Vlakte van de Raan including enhanced water temperatures during the experimental period, increased oxidation by bioirrigation from L. conchilega (evidenced in SPI images; Fig. 3), and the breakdown of fresh organic matter. High levels of the latter (implied by chl a concentrations; Table 1) coupled with elevated faunal and microbial activity due to increased temperatures are the likely mechanism for the high SCOC and mineralization found in this study.

Despite the lower trawl penetration into the sediment (∼1 cm Vlakte van de Raan versus ∼4 cm Frisian Front; Depestele et al., 2018), our study produced a stronger statistically significant reduction in SCOC (Tickler: 57 %) compared to Tiano et al. (2019; Tickler: 41 %; Pulse: 33 %). Significant reductions in total mineralization, denitrification, and nitrification after the Tickler treatment suggests that trawling reduced benthic nutrient removal and recycling functions in this study. Furthermore, higher percentages in the anoxic surface area from the SPI images taken in adjacent sites (significant for AD1) may suggest increased mineralization in undisturbed areas during the course of the experiments (Table S5). Mass budget estimates for mineralization were also higher on average at T1 for AD1 and AD2, though this was not statistically significant. If the benthic metabolism in untrawled areas did increase during the experimental period, the trawling effects detected in this study (lower metabolism after trawling) may be considered to be conservative estimates.

The large changes in SCOC may be explained by the decrease in fauna such as L. conchilega, which can be highly influential to SCOC and other biogeochemical fluxes (Braeckman et al., 2010). In addition to respiration, when L. conchilega (and other infauna) irrigate their burrows, they can drastically expand the surface area available for sediment–oxygen exchange (De Smet et al., 2016b; Kristensen and Kostka, 2013) with major impacts to biogeochemical dynamics (Olsgard et al., 2008; Braeckman et al., 2010). These results also suggest that areas with L. conchilega may be vulnerable to relatively shallow (∼1 cm) seabed disturbances as their tubes extend from over 10 cm within the sediment to above the seabed surface. Contrary to this, faunal-mediated biogeochemical functions in the Frisian Front are relegated to organisms such as the burrowing mud shrimp, Callianassa subterranea, which does not exhibit any protruding sediment surface structures and resides deeper within the seabed (Amaro et al., 2007; Tiano et al., 2020). These organisms are probably less affected by trawling despite their habitats displaying softer sediments which are prone to greater trawl penetration (Depestele et al., 2018; Tiano et al., 2019, 2020).

L. conchilega in untrawled sediments had a clear effect on O₂ and NO₃ fluxes, though these relationships disappeared after disturbance (Fig. 6, Table 2). The change in these relationships probably came from a combination of reduced average tubeworm densities (consistently lower after trawling but only significant after the Pulse-ON treatment) and direct trawl-induced habitat alterations. These results support the assertion by Hale et al. (2017) that trawling decouples faunal-mediated biogeochemical processes.

The Pulse-ON treatment, which led to the biggest decrease (72 %) and lowest L. conchilega densities after T1, was also the only experimental treatment that did not produce a significant increase in porewater nitrate concentrations. Previous research suggests a density-dependent enhancement by L. conchilega on denitrification and nitrification processes (Aller, 1988; Braeckman et al., 2010; Foshtomi et al., 2018). Trawling has also been shown to decrease denitrification by altering the sediment layers that optimize denitrification and reducing total mineralization rates (Ferguson et al., 2020; De Borger et al., 2021). It is possible that both mechanical trawl disturbance and a loss (Pulse-ON) or retention (Pulse-OFF, Tickler) of functional species roles after trawling led to the enhanced sediment nitrate concentrations, though it is difficult to pinpoint the exact cause of this increase.

4.4 Benthic community observations

Our study found a significant decrease in L. conchilega densities with the Pulse-ON treatment and not in the other trawl treatments; however, it is unlikely that electrical stimulation caused this reduction. Previous research has not found any direct mortality of fish or benthos associated with electrical stimulation used in the pulse trawl fishery (Polet et al., 2005; van Marlen et al., 2009; Soetaert et al., 2015b, 2016; ICES, 2018, 2019, 2020; Boute et al., 2020). Moreover, electric stimulation from flatfish pulse trawls cannot physically remove L. conchilega tubes from the sediment (Boute et al., 2020; ICES, 2020; Tiano et al., 2021). Our results, however, point to a higher physical disturbance occurring in the Pulse-ON areas as being the likely culprit for reduced L. conchilega densities (Fig. 5, Table 2).

The Pulse-ON treatment also caused a significant decrease in Eumida sanguinea, a species closely associated with L. conchilega (Rabaut et al., 2008), though other taxa, such as Abra alba, only showed notable decreases in response to tickler trawling. Nonetheless, when accounting for total macrobenthos densities, the only significant effect came from the Tickler treatment (58 % decrease), reflecting the more consistent effect associated with tickler beam trawl gears (Depestele et al., 2018; Tiano et al., 2019). Pre-trawled habitats in this study linked L. conchilega to the densities of other taxa and species richness, though this relationship was no longer detected after experimental trawling (Fig. 6, Table 2). This suggests that, although tubeworm densities were impacted to some extent, trawling disproportionally affected the corresponding taxa. Rabaut et al. (2008) reported significant changes to an intertidal community dominated by L. conchilega after beam trawling but not to the tubeworm densities themselves. Our results show a similar removal of species associated with L. conchilega (E. sanguinea in particular), though our findings of decreased tubeworm densities support the evidence that L. conchilega shows a certain level of sensitivity to physical sediment disturbances (Callaway et al., 2010). L. conchilega has been found to successfully recolonize previously disturbed habitats within 1–2 years, suggesting high resilience but not necessarily resistance to anthropogenic impacts (Callaway et al., 2010).

5 Conclusions

This study found various effects from conventional tickler-rigged beam trawls and electric pulse trawl techniques. Though pulse trawls (and any other bottom fishing technique) can undoubtedly cause significant habitat alterations, these effects tend to be quite variable (Tiano et al., 2019, and this study). The more consistent effects of tickler trawls are, on average, higher than that of pulse trawls as shown in the sediment, geochemical, and species community level changes seen in this study. However, we also show that pulse trawls may be at least equally as impactful on densities of certain benthos. It is likely that the lower average penetration depth of the lighter pulse trawl gears are less damaging to certain benthic communities compared to tickler trawls (Hiddink et al., 2017; Depestele et al., 2018); however, this reduced effect is dependent on species assemblages and sediment type. Compared to any difference in direct impacts between pulse and beam trawls (due to different sediment penetration depths), the potential for a reduced environmental footprint from pulse trawls (compared to beam trawls) stems more from its higher catch efficiency which leads to reduced time needed to catch their quotas (less time fishing) when exploiting North Sea sole (Solea solea; Rijnsdorp et al., 2020b; ICES, 2020).

The high SCOC values found at our experimental location imply that this community is an important site for carbon degradation (Stratmann et al., 2019). The trawl-induced declines in SCOC and total mineralization suggest that some disturbed locations in this community lost over half of their ability to mineralize organic matter. This highlights the propensity for bottom trawling to reduce carbon cycling in marine sediments and shift the mineralization processes to the overlying water (Tiano et al., 2019). With enhanced mineralization rates in the water column (due to higher oxic conditions compared to inside the sediment), trawling may lead to the removal of sedimentary carbon, increased nitrogen retention in the seabed (decreased denitrification), and a reallocation of nutrients away from trawled sediments (Tiano et al., 2019, 2021; Ferguson et al., 2020; De Borger et al., 2021).

Results in this study show that even highly dynamic coastal areas can be susceptible to significant trawl-induced changes depending on the type of benthic community present. Biogenic reef communities like the ones found in the Vlakte van de Raan can be highly productive habitats with enhanced biogeochemical characteristics. With these features, the potential for disturbance-induced changes becomes greater if keystone species such as L. conchilega are affected. When assessing the vulnerability of marine habitats, environmental managers and policy makers need to consider abundant ecosystem engineers as they perform important ecological and biogeochemical functions.