Quantifying habitat preference of bottom trawling gear

Quantifying habitat preference of bottom trawling gear. Continental shelves around the world are subject to intensive bottom trawling. Demersal ﬁsh assemblages inhabiting these shelves account for one-fourth of landed wild marine species. Increasing spatial claims for nature protection and wind farm energy suppresses, however, the area available to ﬁsheries. In this marine spatial planning discussion, it is essential to understand what deﬁnes suitable ﬁshing grounds for bottom trawlers. We developed a statistical methodology to study the habitat preference of a ﬁshery, accounting for spatial correlation naturally present in ﬁsheries data using high-resolution location data of ﬁshing vessels and environmental variables. We focused on two types of beam trawls to target sole using mechanical or electrical stimulation. Although results indicated only subtle differences in habitat preference between the two gear types, a clear difference in spatial distribution of the two gears was predicted. We argue that this change is driven by both changes in habitat preference as well as a change in target species distribution. We discuss modelling of ﬁsheries’ habitat preference in light of marine spatial planning and as support in benthic


Introduction
Continental shelves around the world are subject to intensive bottom trawling. Demersal fish assemblages inhabiting these shelves account for one-fourth of landed wild marine species (Amoroso et al., 2018). The North Sea is part of the European continental shelf and is extensively trawled by different fishing gears. Increasing spatial claims for nature protection and wind farm energy suppresses, however, the area available to fisheries that may hamper the ambition to increase food production from marine environments. As such, it is essential to understand what defines suitable fishing grounds for bottom trawlers to allow for informed decisions on the location and design of windfarms and marine protected areas (Stelzenmuller et al., 2008). This understanding is not only key for the spatial planning debate but also to illustrate that fishers are bound to certain hotspots in space and do not have the ability to move their activity without reducing the viability of their business. Beyond spatial planning, discussing the footprint of bottom fishing and comparing the impacts different types of fisheries have on seafloor integrity have increased in attention in recent years. This is likely driven by the Marine Strategy Framework Directive (MFSD) (EC, 2008) prescribing that member states in the EU need to ensure that seafloor integrity is at a level that ensures functioning of the ecosystem (Descriptor 6). Habitat characteristics, ecosystem functioning, and fishing impact are intertwined and hence all need to be appropriately addressed to evaluate the sixth MFSD descriptor. The societal debate also focusses on the ratio between seafloor impact (i.e. area impacted by a bottom trawl gear measured in km 2 ) and the amount of animal proteins obtained in the fishing activity, where generally fisheries with a lower ratio (i.e. footprint) are preferred. Studying habitat preference of bottom trawlers thus advances our understanding of benthic impact and ecosystem functioning as well as our ability to predict fishing impact at small spatial scales relevant for seafloor integrity, spatial planning, and fisheries footprint studies.
The spatial distribution of bottom trawlers differs among metiers (Eigaard et al., 2017;ICES, 2018a) and reflects the broadscale distribution patterns of the targeted marine resources (ICES, 2018b, c). At a fine spatial scale ($1 km scale), the distribution of a fishery is often patchy (Rijnsdorp et al., 1998;Murawski et al., 2005;Lee et al., 2010;Ellis et al., 2014), reflecting habitat heterogeneity (van der Reijden et al., 2018). Habitat heterogeneity will affect the local abundance of target species and determine the possibility to safely deploy a bottom trawl. As such, certain habitats are preferred over other habitats, i.e. fished with higher intensity, as they yield higher catch rates. This is referred to here as habitat preference. Because sensitivity of the seafloor and the benthic communities differs across habitats, knowledge on habitat preference is important for the assessment of fisheries impact (Kaiser et al., 2006;Lambert et al., 2014;Hiddink et al., 2017;Pitcher et al., 2017;Rijnsdorp et al., 2018;Hiddink et al., 2019). Disentangling whether the spatial distribution of the fishing fleet is defined by either fishing gear type or target species habitat preference is challenging. However, such information is vital for fisheries management and spatial planning because changes in gear type (e.g. due to innovation or policy changes) may result in changes in the distribution of fishing effort and may alter the interactions with other stakeholders using the marine environment. Furthermore, understanding how habitat preferences change with modifications made to fishing gears could lead to more tailored gear design that reduces the seafloor impact. In fishing gear technology, one needs to be able however to objectively evaluate how changes in fishing gear design result in changes in fishing footprint, an approach for this is presented in this study.
This study focusses on the beam trawl fishery in the North Sea targeting sole (Solea solea) and plaice (Pleuronectes platessa). Beam trawls have been used from the 1960s onwards when they replaced the otter trawl as the dominant gear to catch sole and plaice (Rijnsdorp et al., 2008). Although the large-scale spatial distribution of the beam trawl fisheries has shifted at decadal scales (Van der Pol et al., in prep.), the fine-scale distribution of fishing activity has been very stable since the 2000s . The fishery in the southern North Sea that primarily targets sole is known to be located in warmer, shallower, dynamic areas in the southern North Sea where sand ridges are common (van der Reijden et al., 2018). In between these ridges, fishers tend to achieve good catches and therefore return to such grounds year after year. Fishers tend to avoid areas with coarser substrate (van der Reijden et al., 2018;Hintzen et al., 2019). Furthermore, their distribution is affected by the availability of sole and plaice quota (Poos et al., 2010).
Traditionally, beam trawlers have fished with several tickler chains in front of their nets and a steel beam to keep the net open. Owing to increasing oil prices in the 2000s, the industry replaced the steel beam with a hydrodynamic foil (Sumwing) to reduce fuel consumption (Turenhout et al., 2016b, Depestele et al., 2019. A second innovation was the pulse trawl which replaced tickler chains with electrodes that emit electric pulses (van Marlen et al., 2014, Depestele et al., 2019. Although commercial electric fishing has been banned in EU waters since 1998 (EU, 1998), a study fleet received a temporary exemption (Haasnoot et al., 2016). Under this exemption, a large part of the Dutch beam trawl fleet switched to pulse fishing. Pulse trawling turned out to improve the economic profitability owing to a lower fuel consumption and improved catch efficiency for sole, although the catch efficiency for plaice and other species was reduced (van Marlen et al., 2014;Turenhout et al., 2016a;Poos et al., 2020). A large part of the Dutch beam trawl fleet switched to pulse fishing between 2009 and 2015 and by 2016 about 95% of the Dutch sole quota were caught with the pulse trawl (ICES, 2018d). This large-scale switch to pulse allows us to study differences in habitat preference affected primarily by the change in gear design. Stakeholder information suggests that pulse fishers started to use different habitats compared to their distribution while using tickler chains (ICES, 2018d).
In this article, we study the interactions between gear developments, habitat heterogeneity, and habitat preference. Habitat preferences can be studied making use of statistical models (Rushton et al., 2004;Bertrand et al., 2016) that relate spatial count data to environmental variables. In this study, fishing vessel GPS data [vessel monitoring by satellite (VMS)] provided a detailed view of the spatial distribution of the fishing fleet. The micro-scale (tens of metres) at which the VMS data are available allowed testing for subtle differences in habitat preference when comparing two gear types. A study by van der Reijden et al. (2018) indicated already that among other factors, depth profile, sediment type, and natural disturbance were key indicators to explain habitat hotspots for beam trawl fishers. Noting that bathymetric information is available at very fine spatial scales (van der Reijden et al., 2018), our ability to define habitat preference at micro-scale (tens of metres) is limited more by sediment and natural disturbance data, which are only available at lower resolutions (Wilson et al., 2018).
The results show that there is a substantial difference ($50%) in spatial distribution between the pulse and tickler chain fisheries, where the first prefers habitat with higher gravel content and more elevated areas (relative to its surroundings) and shows fishing activity in between sand ridges. This shift in spatial distribution has caused some habitats to be more intensively impacted than before the switch to pulse gears while other areas are less frequently trawled. The benthic impact associated with the distributional shift is discussed. We argue that there is a clear role for habitat preference modelling in spatial fisheries management, such as tailoring the design of marine protected areas and supporting benthic impact assessments.

Case study
We study the Dutch beam trawl fleet targeting mainly sole in the southern North Sea. In this mixed fishery, sole is caught in a mixture with plaice, turbot, brill, and dab. Beam trawling has been used in the North Sea since the 1960s and has been a dominant fishery ever since in the Dutch fishing sector. On each side of the vessel, a 12-m wide steel beam fitted to a shoe on each side of the beam is dragged over the seafloor (see Eigaard et al., 2016 for a graphical representation of the gear). The beam fixes the horizontal net opening and allows the fisher to deploy tickler chains perpendicular to the towing direction to chase flatfish from the seafloor into the net and increasing the catch efficiency of the gear (Daan, 1997;Rijnsdorp et al., 2008). The pulse trawl is similar in design as the traditional tickler chain beam trawl but uses longitudinal cables that emit electric pulses which invokes a cramp response that immobilise the fish (de Haan et al., 2016;Soetaert et al., 2019). The pulse trawl gear is lighter and penetrates less deep into the sediments than the tickler chain gear and hence may provide access to softer habitats or more coarse habitats as the gear can be more easily pulled over these habitats (Depestele et al., 2019. The fishery targets the same flatfish types, though catchability has increased fishing with the pulse trawl for sole and has been reduced for plaice compared to the tickler chain gear (Poos et al., 2020). Both the tickler chain and the pulse trawl fleet belong to the same fleet segment TBB (beam trawls). We refer to the beam trawl fleet when speaking of both the tickler chain and pulse trawl combined while the two gear types are singled out when we discuss the differences between the two gear types.

Spatial fisheries data
The analyses included VMS data and mandatory catch and effort logbook data from all Dutch flagged vessels that fished during the transition from the tickler chain gear to the pulse trawl gear. Only vessels with engine power >221 kW were selected. These larger vessels were not allowed to trawl inside the 12-nm zone and in the so-called Plaice Box (see Figure 1). The study area was delineated by the 51 latitude line in the south and the 56 latitude line in the north (55 latitude west of 5 longitude), excluding the 12-nm zone and Plaice Box, corresponding to the area where beam trawlers are allowed to fish with 80-mm codend mesh size. The years 2009-2017 were included, as the first vessels switched to pulse trawling in 2009. From 2015 onwards, pulse trawling represented $65-70% of the total area fished by the entire beam trawl fleet. During the study period, fishing effort declined from $14 000 fishing days in 2009 to $11 000 fishing days from 2014 onwards.
VMS observations (i.e. pings, a signal from a fishing vessel transmitted via satellites to a ground station) include information about vessel name, speed and heading over ground, a date-time stamp, and a GPS position. Fishing activity was defined based on speed profiles (Poos et al., 2013) and non-fishing pings were excluded from further analyses. Over time, the ping frequency of VMS has increased with more pings being submitted at 30-60min intervals rather than the common 2-h interval rate (from 80% 2-h interval in 2009 to 63% by 2017). Daily catch and effort logbook records provided information on vessel length, engine horse power, trip information such as gear and mesh size used, and the catch by species. Usually, a fishing trip lasts around one working week. Although the tickler chain and pulse trawl share one common gear code in the logbooks, i.e. "TBB", an independent database was included which contains more details on gear specifications and their introduction date. Each fishing trip gear usage was further validated by analysing mean fishing speed Quantifying habitat preference during a trip as this was found to be highly indicative of gear usage (Poos et al., 2020).

Creating count data
VMS and logbook data were carefully scrutinized for erroneous entries, following Hintzen et al. (2012). The final dataset included data from 70 vessels that were active during the transition from tickler chain to pulse trawl. This dataset was divided into two subsets: half of the vessels were randomly selected from the dataset and used for model fitting (training data), and the remaining vessels were used for cross-validation. Both subsets span all years and gears. The study area was divided into squares (i.e. grid cells) measuring one by one minute longitude-latitude ($2 km 2 ). For both the pulse trawl and tickler chain gear types separately, VMS pings within each square were summed and used as count data (i.e. counts) for the distribution model ( Figure 1). As such, the number of VMS pings within a grid cell was used as a response variable in the statistical model.

Environmental covariates
A priori seven covariates were selected to be included in the model, which were shown to be relevant in determining the distribution of fishing effort of bottom trawl fisheries (van der Reijden et al., 2018). Within the time-frame of this study the abiotic covariates are assumed to have remained constant. The following environmental variables were attributed to each of the grid cells in the study area: proportion gravel, proportion mud, proportion rock, depth, and mean tidal velocity as indictor of bed shear stress (bedstress). Note that proportion gravel, mud, and sand would sum to 100% and, for this reason, sand was excluded from the analyses to prevent having very high co-linearity among these. For each grid cell, distance to nearest Dutch harbour was calculated from the Euclidian distance between a grid cells' midpoint and the GPS midpoints of the Dutch harbours. For depth, Bathymetric Positioning Index (BPI) was used, being a measure of the depth at a specific location relative to the depth in the surrounding grid cells maximum r km away ( Figure 2). Two BPI values were used for each grid cell, with different r values: BPI 5 with r ¼ 5 km (small-scale features) and BPI 75 with r ¼ 75 km (large-scale features). These values were taken from (van der Reijden et al., 2018). Gravel, mud, rock, and bed shear stress estimates were obtained from (Wilson et al., 2018). Gear type (i.e. pulse trawl or tickler chain) was used as a covariate in the model (as a factorial covariate), which allows testing if there was a difference between the two gear types.
Furthermore, the inverse of average VMS interval time in each grid cell was used as a model offset. This offset was included to account for the change in a number of observations there are in the raw dataset owing to the decrease in interval rate of the VMS data from 2009 to 2017. There are more VMS observations when the interval time is low, resulting in higher VMS pings in a grid cell. This increased amount of pings should not be interpreted as an increase in fishing effort. The offset is a means to standardize the number of pings in each grid cell irrespective of the interval rate of the VMS data. Other available covariates such as sand proportion, depth, wave orbital velocity and intermediate BPI ranges were not considered owing to large co-linearity with the other covariates.

Model structure
We use a statistical framework to model the spatial distribution of fishing effort and hereby being able to objectively separate the response of fishing behaviour to different habitat characteristics. We used the Integrated Nested Laplace Approximation package in R, allowing for the inclusion of spatial latent fields to capture (residual) spatial autocorrelation in observations (Rue et al., 2017). Both spatial and temporal correlation are, by the way VMS data are collected, present in the dataset, i.e. the location of consecutive VMS pings depends on the previous position and the maximum speed of a fishing vessel that allows it to move to another area. A correction for this correlation is necessary to prevent drawing incorrect conclusions on the preference of a fishing gear to a certain habitat characteristic. The spatial autocorrelation is modelled using a Matérn correlation function that is commonly used to model the statistical covariance between observations of two data points that are x km away from each other. Estimating the parameters in the Matérn correlation function requires dividing the study area into a large number of non-overlapping triangles, called a mesh. We used a mesh with an average leg length of $14 km, $ 1 =4 of a degree longitude. The mesh, with the maximum edge of 25 (i.e. the maximum leg length on the edge of the mesh) and a cutoff of 1 (i.e. the minimum leg length between data points), covers the North Sea delineated by the 51 latitude line in the south and the 56 latitude line in the north (55 latitude west of 5 longitude), expanding over the edges of the study area. The mesh was used as input to the stochastic partial differential equations approach to estimate spatial correlation in continuous space (Lindgren et al., 2011).
For the distribution of the response variable, six options were explored: the Poisson, over-dispersed Poisson, Negative Binomial, Zero-inflated Poisson, and Zero-inflated Negative Binomial. Here, we started with a model including linear terms for all covariates and a separate spatial latent field for each gear type and selected one of the six statistical distributions leading to the lowest Watanabe-Akaike information criterion (WAIC, Watanabe, 2010) and Bayesian information criterion (BIC, Ghosh et al., 2007) score while checking observed vs. fitted values to be reasonable. Both WAIC and BIC are statistical representations of goodness of fit and lower values indicate better statistical fits to the observations.

Model selection
First, the most appropriate distribution for the response was selected fitting a full model including a spatial correlation factor for each gear type. Next, we extended the linear terms for the covariates with multi-order (first to seventh) polynomials. Also, the interaction between gear type and each polynomial function of the covariates was included through testing if confidence intervals of the respective covariate or covariate interaction of the more complex model were outside the bounds of the less complex model. In the third and final iteration, all covariates were evaluated adding or reducing the polynomial degree by one, as changes in step 2 could have resulted in small changes in the fit of one of the other covariates. This led to a final model generically formulated in (1). Best models were selected based on fit (visual), lowest WAIC, and ability to estimate the cross-validation counts (R 2 on observed $ fitted linear model).
The spatial correlation is estimated for both the pulse and the tickler chain gear types.
where Y represents the predicted total number of VMS pings (counts), I is the intercept term, and f p (Á) is a pth-order polynomial, where p can vary between 1 (linear term) and 7 (seventh order polynomial). Spatial correlation is defined by Aðs; s 0 Þ uðs 0 Þ where Aðs; s 0 Þ represents the projection matrix to project the process from the mesh nodes to the VMS locations. u s 0 ð Þ represents the random field at the mesh nodes. Note that Y represents either counts of pulse trawl or tickler chain owing to the Gear factorial covariate in the equation.
Finally, an analysis was undertaken adding covariates one at a time, keeping in each iteration the covariate that explains most of the remaining variance (and is associated with the lowest WAIC value). This analysis indicates which covariate is most important in explaining the differences between habitat preference of the two different gear types. Covariates that did not result in a significant reduction in WAIC value were omitted from the model used to make predictions on the standardized spatial distribution of both gears.
After the model was fitted to the data, 10 000 new sets of model parameters were simulated using the uncertainty estimates of these parameters and their joint posterior distributions that describes the correlation between all estimated model parameters. These parameter sets were used to obtain 10 000 predictions of VMS pings for each grid cell. VMS pings are equal to fishing effort here as the predictions are standardized for interval rate and as such each ping represents 1 h worth of fishing activity. To be able to quantify differences in the estimated relative distribution of the tickler chain and pulse trawl, each of the 10 000 samples were scaled by its maximum value. As predicted values depend on the total effort of each of the gears, which is different for the pulse and tickler chain fisheries, the scaling is necessary to account for the effort difference. The proportion of effort allocated was calculated for tickler chain and pulse trawl and compared grid cell by grid cell, assuming that if either pulse trawl or tickler chain trawl intensities were outside 95% of the predictions of the other gear, they would be considered statistically different. In those cases where 97.5% of the samples of pulse trawl had a lower value than the lowest 2.5% of tickler chain fishing, the area was marked as significantly favoured by the tickler chain. The same was applied to cases where 2.5% of the pulse trawl samples had a higher value than the lowest 97.5% of the tickler chain fishing. These areas were marked as being statistically significantly favoured by the pulse trawl gear.

Model fit
The observed and estimated counts are in good agreement and observations fit well within the uncertainty bounds as estimated in the model (Figure 3e and f). Although the spatial distribution of both the pulse and tickler chain trawling fleets in the model training data is markedly different from the data used for cross- Quantifying habitat preference validation (see Figure 1), the model fits well to the crossvalidation data. There is a slight underestimation of grid cells with 2-5 counts and overestimation for grid cells with 5-15 counts in the cross-validation data (panel f), potentially caused by temporal effects that were ignored in the model.
Each of the covariates were modelled with increasing flexibility in polynomial design. Model selection let to conclude that BPI 5 shows the observed log-frequency of 0, 1, 2, etc., counts in the dataset in red/blue bars, the black dots represent the estimated frequency of these counts by the model. Panels (c) and (d) show the 1:1 relationship between total observed and estimated counts. Panels (e) and (f) show the counts vs. log-frequency including a 95% confidence bound (dashed lines).
was modelled as a second order polynomial, BPI 75 as a fourth order, mud as a linear term, gravel as a third order, distance to harbour as a fourth order, rock as a third order and bedstress as a fourth order polynomial. For each of the covariates, the model fitting procedure yields a relationship between the covariate itself and the preference for it for each gear type (Figure 4). If values exceed y ¼ 0, including the confidence bounds, it is interpreted as a preference. This preference varies over the range of the covariate itself and can, e.g. illustrate negative preference at negative relative depths (BPI75) and positive preference at positive relative depths (BPI75).

Quantifying habitat preference
Visual inspection of the residuals, including fitted loess smoothers through the residuals, did not highlight any apparent pattern or deviation of balanced residuals. Spatial residuals showed randomness on the model fitted data and some minor patterns for the cross-validation fit.

Habitat preference
There is a clear preference for both gear types to fish at slightly elevated areas (BPI 75), in between sand ridges (BPI 5), and in areas with higher bedstress (Figure 4). Fishing in areas with higher gravel content and (to a lesser extent) more rocks is generally avoided. No preference for a specific range of mud fractions was found. For the distance to harbour variable, preference was rather similar over a broad range of distances between 50 and 300 km from harbour representing the fishing area for sole outside the 12 nm zone and below the northern border of the study area.
The interaction between habitat variable and gear type was significant for the habitat variables BPI 5, BPI 75, bedstress, rock and gravel, while no significant difference was observed for distance to harbour and mud ( Table 1). The relative depth, measured over a 75-km radius (BPI 75) is the most important explanatory variable, reducing the WAIC by 146 points. This represents 72% of the overall reduction in WAIC compared to the model without covariates. Furthermore, including distance to harbour and bedstress, up to $95% of the reduction in WAIC is explained. A final 3% is explained by adding percentage gravel. Adding rock, BPI 5 and mud did not improve the model and resulted in a minor increase in WAIC.
Comparison of the preference curves between gear types shows that pulse trawling has a slight but significant stronger preference for the intermediate depths of BPI 75 and a lower preference for deeper and slightly elevated areas than tickler chain gear. Deeper troughs (BPI 5) are preferred by both types of fishers, while pulse fishers also have a preference to fish at tops of sand ridges (BPI 5). The additional explained variance by the BPI 5 covariate is low however. The preference curve for bedstress of the pulse gear is shifted to a slightly higher bedstress. For gravel habitats, pulse gear has a slight but significant preference for low gravel fractions (i.e. preference is just above the y ¼ 0 line in Figure 4). Tickler chain gear tries to avoid gravel habitats under all circumstances.
A model consisting of fixed effects BPI 75, distance to harbour, bedstress, and gravel with a VMS interval rate set to 1 h was used to predict the pulse and tickler chain counts for each of the grid cells in our study area (hereby dropping rock, BPI 5 and mud given the minor contribution to overall gear differentiation and potential inflation of confidence intervals). Although the overall distribution of fishing activity is similar between the two gear types, there are areas where pulse trawls are more active, such as in the southwestern part of the North Sea, while tickler chains are more dominant in the south eastern part (i.e. German Bight) ( Figure 5).

Predicting spatial distribution
When determining the ratio in predicted pings between pulse and tickler chain gear, areas that show marked differences between the two gear types are given in darker colours ( Figure 5). The middle and right-hand panels show the areas that are significantly different for both gears (middle panel shows significantly higher pulse, right-hand panel shows significantly higher tickler chain). The area where tickler chain or pulse activity differed significantly amounted to 49.7% of the study area, with 16.8% of the area associated with higher pulse activity and 32.9% associated with higher tickler chain activity. A breakdown of the ratios is given in Figure 6, showing that 80% of the grid cells with significantly different intensity are associated with ratios between 3:1 and 1:3. Only 7% of the grid cells not significantly different in intensity are associated with ratios outside this intensity range. Table 2 shows the main characteristics of areas where tickler chain or pulse activity differed significantly compared to the average of the study area. This shows that pulse trawling is significantly more active in areas with higher gravel content, in more elevated areas compared to its wider surroundings (BPI 75) and in areas with higher bedstress (southern North Sea, which is also located, on average, closer to shore than areas further north). Tickler chain fishers fish in areas with lower gravel content, on less elevated grounds compared to its wider surroundings (BPI 75) and in areas with lower bedstress. The tickler chain fishers do show a preference for areas with higher bedstress (see Figure 4) and both groups prefer to fish in between sand ridges (BPI 5) rather than on the slopes or top.
The difference in aggregation of fishing activity between the pulse trawl and tickler chain fishers is best illustrated when we assume that an average grid cell is fished with one unit of effort. Pulse fishers deploy around three units of effort in each of the significant grid cells (i.e. cluster a large part of their effort in these areas: 50% of all effort units in 16.8% of the study area). Tickler chain fishers only deploy an additional 0.35 units of effort in the areas they have significant higher counts (i.e. show a more evenly distributed effort all over the fishable areas: 44% of all effort units in 32.9% of the study area).

Spatial distribution of bottom trawl fishery
The beam trawl fleet (both tickler chain and pulse trawl) preferentially selects elevated landscapes (i.e. higher BPI at large spatial scales, BPI 75), substrates with low gravel content and in-between sand ridges rather than on the top (lower BPI at small spatial scales, BPI 5). These results agree with the findings of van der Reijden et al. (2018) and suggest that habitat characteristics of fishing hotspots apply to areas with lower fishing intensity too. The beam trawl fleet seems to avoid either rocky or muddy substrates. Furthermore, there is no clear preference to fish closer or further from shore up till $300 km where after preference shows a clear dip, associated with the area north of 56 where larger mesh sizes are obliged, limiting the ability to catch sole.

Differences between pulse trawl and tickler chain trawl fisheries
Predictions showed that the spatial distribution was significantly different for the two gear types in almost 50% of the entire study area. The tickler chain fishery was most abundant both in the southern North Sea and in the German Bight in the eastern part of the North Sea north of the Netherlands and Germany. The pulse fishery was more concentrated in the southern North Sea, closer to the United Kingdom 12-mile zone. This southern area is characterized by higher bedstress and has more gravel patches and sand ridges, which are reflected by the higher variability in bathymetric position index at large spatial scales (BPI 75). Since the pulse fishery is more concentrated in the southern North Sea, the region closest to most Dutch harbours, the overall distance to harbour is lower for pulse fishers too. This spatial shift between tickler chain and pulse trawlers could be explained by the change in catchability of sole. Compared to the tickler chain trawlers, the pulse trawlers have a substantially higher catch efficiency for sole (and lower efficiency for plaice) (Poos et al., 2020). Sole abundance in the German Bight has declined since the 90s and increased in the southern North Sea (Vansteenbrugge et al., 2020); this change in the main distribution area of sole in the southern North Sea likely explains the observed spatial shift of effort. It is unlikely that the shift is a result of a change in habitat preference, as in general, habitat preference is similar for both gear types. The shape of the preference curves for abiotic variables such as BPI and gravel are very similar (Figure 4, e.g. for gravel the fitted preference curves are very close to each other over the entire range), though do occasionally differ significantly in absolute terms where, e.g. pulse trawl has a higher preference for BPI 75 in the range of 0-20 m.

Consequences of the transition from tickler chain to pulse trawl
The transition from tickler chain beam trawls to pulse trawls led to a higher catch efficiency, a lower towing speed, and a reduction in the impact on the benthic ecosystem (Poos et al., 2020;Rijnsdorp et al., 2020). Our study showed that pulse fishers spent around three times the effort per grid cell compared to an average grid cell in the study area, in areas where they have a significantly higher preference compared to the tickler chain fishers. Tickler chain fishers only spent an additional 35% of their effort per grid cell compared to an average grid cell in the study area. As such, they spatially aggregate their effort to a higher degree than the tickler chain fishers do, which implies that these areas are fished at higher fishing intensities (Ellis et al., 2014;Hintzen et al., 2019). The pulse preference areas are associated with higher gravel content. Coarser sediments have been shown to be more vulnerable to fishing Rijnsdorp et al., 2018) because they generally contain more sessile and longer-lived organisms. These organisms decline more rapidly in biomass under higher fishing pressures compared to communities with mobile and short-lived organisms (Hiddink et al., 2019). Furthermore, if pulse fishers moved to previously unfished areas, a substantial reduction in benthic biomass can be expected in those areas (Sciberras et al., 2018). However, the impact depends not only on the trawling intensity but also on the penetration depth Sciberras et al., 2018) and sensitivity of the benthic community, being related to the amount of natural disturbance (van Denderen et al., 2015b; Rijnsdorp et al., 2018;Hiddink et al., 2019). Bedstress caused by currents is higher in a large part of the pulse fishing area, mostly in the southern North Sea, compared to the tickler chain spatial distribution. The penetration depth of the pulse trawl is less than half the penetration depth of the tickler chain beam trawl, and depletion rates of epibenthos imposed by pulse trawls are $50% less than tickler chain beam trawls (Depestele et al., 2019). Indeed, direct mortality imposed by pulse trawling is less than by tickler chain beam trawling (Bergman and Meesters, 2020) and Rijnsdorp et al. (2020) demonstrated that overall benthic impact of pulse trawling was lower than that of tickler chains. In addition, further aggregation of fishing effort by the pulse fishers implies that a larger proportion of the seafloor outside of the preference areas remains unfished or is fished with lower intensity.

The importance of scale in habitat analysis
Here, we analysed habitat preference at a spatial resolution of 1 minute Â 1 minute. Detailed studies in the southern North Sea, however, revealed small-scale heterogeneity in bathymetry and sediment composition with alternating ridges and troughs at scales well within the 1 minute Â 1 minute grid cells (van Dijk et al., 2012;Koop et al., 2019;van der Reijden et al., 2019). This micro-scale heterogeneity may have been used by pulse trawlers. Anecdotal information from fishers indicates that pulse trawlers may have been able to fish areas that could not be fished by tickler chain beam trawlers due to the softness or higher gravel content of the sediment. Higher resolution data on the habitat covariables are required to investigate this hypothesis and assess the consequences on the benthic ecosystem. Furthermore, higher resolution data also help in pinpointing the habitat preference of target species and identify how much habitat is available and may be in need of protection when fish stocks are in decline. For example, the possible occurrence of untrawlable habitat fragments may be relevant to understand the effect of trawling on the population dynamics of sole because these habitat fragments could provide a network of refugia where sole may have been safe from exploitation. At the same time, being yet unable to identify the exact size and characteristics of fisheries hotspots hampers decision-makers to accurately value fishing grounds in a trade-off with other uses such as nature reserve or offshore energy farms.

General use of habitat modelling
The methodology used in this study shows the added value of fitting habitat preference models to VMS data, which is now routinely available for a large part of the global fisheries (Amoroso et al., 2018). Models like these allow testing for both spatial and temporal correlation in vessel abundance indicating the variability in time and space of specific fisheries. This information is valuable for bottom impact assessments where frequency of trawling and time for benthic communities to recover play an important role (van Denderen et al., 2015a). Furthermore, they allow for different management strategies to be developed, such as habitat credit systems (Kraak et al., 2012;Batsleer et al., 2018) for which there is a need to quantify the spatial overlap in habitat preference between different fishing gears. A habitat credit system requires information, on a fine spatial scale, on the likelihood of other fishers to use a specific fishing spot, this to set a cost to each individual fishing area. The statistical framework developed here can provide this estimate but is also able to show if nearby areas are accessible to the fishery that could reduce pressure on traditionally heavily fished grounds.
Having detailed overviews of spatial distributions of fishing activity in relation to benthic (micro)habitats allows for the evaluation of ecosystem functioning of these habitats. This is required under the MSDF (EC, 2008) where member states need to bring the seafloor in such a state that it ensures appropriate functioning of the ecosystem. The methodology developed here provides standardized estimates of fishing impact (contrary to raw VMWbased estimates that are spatially correlated) per habitat type, which has a direct link to ecosystem functioning . It furthermore shows the likelihood that fishers are willing to move to other areas, i.e. if habitat preference to a specific abiotic factor is high, the willingness of fishers to move from those grounds is likely to be small. As such, the estimated habitat preference provides a fisher-independent view on where fishers like to fish, regardless of their home port or quota share.
In the spatial planning debate where managers need to decide where to make room for fishing activity and where to locate windfarms and marine protected areas requires reliable information on ecosystem functioning, suitability of areas to serve as wind farms and habitat preference for fish species and their fishers (Stelzenmuller et al., 2008). The latter one is often inferred from recent observed distribution patterns while these could be biased due to other legal restrictions such as spatial closures and quota availability but also variability in the local productivity of fish stocks. For a long-term perspective, these variables need to be eliminated as can be achieved by making use of habitat preference models.
We demonstrated that changes in gear design had a marked impact on fishing distributions. However, these models could also be used to describe historic changes in distributions of fishing activity at small spatial scales in the absence of gear changes. In the marine spatial planning debate, optimizing the allocation of space for different uses such as energy production, nature conservation and fisheries is crucial. With help of habitat preference modelling, displacement of fisheries can be forecasted statistically, although longer-term forecasts would require inclusion of target species distribution. When predicting from the habitat model, a sum of total fishing effort is distributed over the grid cells according to each grid cells' preference, i.e. proportion of effort they will receive. If one wants to study how total effort would be distributed over space if certain grid cells would be closed, due to, e.g. windfarm development, the effort previously attributed to the windfarm area will be distributed over the remaining grid cells according to their preference. Such predictions are essential in decision-making by marine resource managers as decisions on, e.g. spatial closures or wind farm areas exclude other users. If there is a substantial change in target species distribution however, one needs to include this shift before making predictions based on habitat preference models. Even so, other factors such as temperature change or food availability for target species could become relevant covariates to include when studying the distribution of the fishing fleet over longer periods.
Beyond being able to predict spatial distributions from these statistical models is the ability to estimate the uncertainty of fishing intensity in relation to habitat use, i.e. for each grid cell, the uncertainty in habitat preference and hence uncertainty in predicted fishing counts are available. This is currently not possible with maps derived from raw VMS data. Given that VMS-based spatial analyses of fishing activity often include several assumptions (Hintzen et al., 2012), accounting for uncertainty not only reflects reality, it also provides a range in the footprint of bottom trawling fisheries considering the uncertainty.

Data availability statement
Primary VMS data and catch and effort data of the mandatory logbook are subject to confidential, privacy related, agreements. One should contact Sieto Verver, Head of the Centre for Fisheries Research (sieto.verver@wur.nl) for permission using these data.