Beyond connecting the dots: A multi-scale, multi-resolution approach to marine habitat mapping

of interests between and nature are in coastal seas, a growing need for evidence-based marine spatial planning. This requires accurate, high-resolution habitat maps showing the spatial distribution of benthic assemblages and enabling intersections of habitats and anthropogenic activities. However, such detailed maps are often not available because relevant biological data are scarce or poorly integrated. Instead, physiotope maps, solely based on abiotic variables, are now often used in marine spatial planning. Here, we investigated how pointwise, relatively sparse biological data can be integrated with gridded, high-resolution environmental data into informative habitat maps, using the intensively used southern North Sea as a case-study. We first conducted hierarchical clustering to identify discrete biological assemblages for three faunal groups: demersal fish, epifauna, and endobenthos. Using Random Forest models with high-resolution abiotic predictors, we then interpolated the distribution of these assemblages to high resolution grids. Finally, we quantified different anthropogenic pressures for each habitat. Habitat maps comprised a different number of habitats between faunal groups (6, 13, and 10 for demersal fish, epifauna, and endobenthos respectively) but showed similar spatial patterns for each group. Several of these ‘fauna-inclusive ’ habitats resembled physiotopes, but substantial differences were also observed, especially when few (6; demersal fish) or most (13; epifauna) physiotopes were delineated. Demersal fishing and offshore wind farms (OWFs) were clearly associated with specific habitats, resulting in unequal anthropogenic pressure between different habitats. datasets physiotope

Conflicts of interests between economic and nature conservation stakeholders are increasingly common in coastal seas, inducing a growing need for evidence-based marine spatial planning. This requires accurate, highresolution habitat maps showing the spatial distribution of benthic assemblages and enabling intersections of habitats and anthropogenic activities. However, such detailed maps are often not available because relevant biological data are scarce or poorly integrated. Instead, physiotope maps, solely based on abiotic variables, are now often used in marine spatial planning. Here, we investigated how pointwise, relatively sparse biological data can be integrated with gridded, high-resolution environmental data into informative habitat maps, using the intensively used southern North Sea as a case-study. We first conducted hierarchical clustering to identify discrete biological assemblages for three faunal groups: demersal fish, epifauna, and endobenthos. Using Random Forest models with high-resolution abiotic predictors, we then interpolated the distribution of these assemblages to high resolution grids. Finally, we quantified different anthropogenic pressures for each habitat. Habitat maps comprised a different number of habitats between faunal groups (6, 13, and 10 for demersal fish, epifauna, and endobenthos respectively) but showed similar spatial patterns for each group. Several of these 'fauna-inclusive' habitats resembled physiotopes, but substantial differences were also observed, especially when few (6; demersal fish) or most (13; epifauna) physiotopes were delineated. Demersal fishing and offshore wind farms (OWFs) were clearly associated with specific habitats, resulting in unequal anthropogenic pressure between different habitats. Natura-2000 areas were not specifically associated with demersal fishing, but OWFs were situated mostly inside these protected areas. We thus conclude that habitat maps derived from biological datasets that cover relevant faunal groups should be included more in ecology-inclusive marine spatial planning, instead of only using physiotope maps based on abiotic variables. This allows better balancing of nature conservation and socioeconomic interests in continental shelf seas.

Introduction
Continental shelf seas are subject to increasing anthropogenic activities, such as shipping (Sardain et al., 2019), sediment extraction (de Boer et al., 2011), offshore wind farms (OWFs; Grothe & Schnieders, 2011) and industrial fishing (Eigaard et al., 2017). Meanwhile, calls for effective marine conservation are increasing (Johnson et al., 2017), for example, in demanding the establishment of no-use marine protected areas (Costello and Ballantine, 2015). Ecology-inclusive marine spatial planning that balances the socio-economic and ecological interests can help to resolve current spatial conflicts of interest (Kaiser et al., 2016;White, Halpern, & Kappel, 2012). At present, anthropogenic activities are predominantly located in areas where they are most profitable (Grothe and Schnieders, 2011). However, to improve the balance between both socio-economic and conservation interests, ecological knowledge on, for instance, species abundance and diversity, community sensitivity, and ecosystem resilience should also be included in spatial zonation of anthropogenic activities. This requires accurate highresolution maps that capture the spatial heterogeneity, extent, and biological characteristics of different marine habitats (Kaiser et al., 2016;Reiss et al., 2015). These can be translated into maps representing the uniqueness, diversity, vulnerability and resilience of local demersal assemblages and their relation to anthropogenic pressures (Cooper et al., 2019;Kaiser et al., 2016).
Biological sampling (trough research-based trawls, cores, grabs, etc.) provides information on the spatial distribution of benthic assemblages. However, such sampling is very labour intensive and expensive. It furthermore is prone to methodological limitations in capture efficiency, site accessibility, and sampling of all biological strata (Beisiegel et al., 2017;Jørgensen et al., 2011), and to observer bias in taxonomic identification (Reiss et al., 2010). In addition, biological samples are by nature 'station-based' point samples and thus need interpolation to obtain full-coverage predictions of benthic assemblage distributions. Accuracy increases with better interpolation techniques (ICES, 2019a) and higher sampling resolution (Compton et al., 2013;Cooper et al., 2019).
Due to these difficulties of collecting and analysing biological data, high-resolution environmental ('abiotic') variables are often typically used instead. Variables as water temperature, bathymetric depth, salinity and sediment type are typically derived through combinations of remote sensing, modelling and point measurements, and often form the basis for habitat mapping processes. This includes the currently most commonly used habitat mapping approach: hierarchical classification schemes (Strong et al., 2019). These maps are structured in multiple (stacked) hierarchical 'levels'. Lower levels delineate broad habitat types based only on environmental factors, often combinations of bathymetry and sediment composition categories. The more detailed, higher levels subsequently divide these into more specific habitat types, based on additional environmental factors and/or station-based biological information (Galparsoro et al., 2012;Schiele et al., 2014;Strong et al., 2019). However, usage of these maps is often limited to the broad, low-level habitat types, so based on environmental variables only (Andersen et al., 2018;Galparsoro et al., 2012;ICES, 2020).
Here, we define maps that are exclusively based on environmental variables as 'physiotope maps', which are often seen as surrogates for benthic assemblage distributions (Huang et al., 2011;McArthur et al., 2010;Roland Pitcher et al., 2012;Vasquez et al., 2015). However, although benthic communities are strongly affected by environmental conditions, they typically are not defined by static environmental conditions alone (Galparsoro et al., 2012;Stevens and Connolly, 2004). Communities and environmental conditions can be heavily affected by biological processes like predation, competition, or ecosystem engineering effects (Jones et al., 1994;Menge and Sutherland, 1976), and can change under anthropogenic pressures, such as frequent bottom trawling (Kaiser, Ramsay, Richardson, Spence, & Brand, 2000). Moreover, the relevance of different environmental factors varies greatly between spatial scales (Damveld et al., 2018;Lecours et al., 2015), differs between faunal groups (e.g. demersal fish versus endobenthos; Hewitt et al., 2015), and not all relevant environmental variables may be included in physiotope maps. Therefore, physiotope maps at best represent the required biological reality only partially.
Using the North Sea as a case-study, we here explore a novel approach to derive high-resolution habitat maps that combine low resolution, point-based biological data of three faunal groups with higher resolution key environmental variables. The resulting habitat maps were compared with physiotope maps to visualize the contribution of biological input data. In this, the role of fishing intensity as an anthropogenic environmental variable was studied and spatial associations of anthropogenic uses (fisheries, OWFs, Natura-2000 areas) with the identified habitats were quantified. We conclude with a discussion on our findings can improve ecology-inclusive marine spatial planning in coastal seas.

Methods
In our study we compare habitat maps based on biological samples, spatially interpolated using high-resolution environmental gradients with physiotope maps that are solely based on the same environmental gradients (Fig. 1). Sampled species abundances of three faunal groupings were separately clustered into assemblages with a hierarchical clustering. Environmental gradients were determined by reducing 21 environmental variables (Table 1) to less dimensions (main environmental gradients) using Principal Component Analysis (PCA). The pointwise assemblages were interpolated based on these main environmental gradients into full-coverage habitat maps. The resulting habitat maps were then compared to physiotope maps, which resulted from a clustering of the main environmental gradients exclusively, to investigate the added value of biological information in the habitat mapping exercise. We additionally studied the potential of fishing intensity as an explanatory, anthropogenic variable by performing the whole analysis twice: first with the 21 environmental variables, and then with the 21 environmental variables in combination with fishing intensity.

Study area
In this study, we focused on the offshore Central and Southern North Sea, defined by the International Council for the Exploration of the Sea (ICES) as subdivisions IVb and IVc (Fig. 2). To avoid any near-coastal effects of modelled parameters, the national territorial waters (12 nautical miles from shoreline) are excluded from the analysis.

Demersal fish
We obtained demersal fish abundances from the annual North Sea Beam Trawl Survey (BTS). This survey samples demersal fishes with a beam trawl (8 m width, 4 cm mesh size), at relatively fixed stations throughout the North Sea (Fig. 3A) by performing 30-min hauls at 5 knots (covering ~ 0.037 km 2 ). All caught individuals were identified and counted (see ICES, 2009 for detailed methodology). We obtained all BTS-data from the DATRAS website for the years 2008-2015 (Millar et al., 2019). Only valid hauls with species recordings, a maximum distance of 10 km covered, and located within the study area were included. Species with missing abundances were removed, as were all non-fish species and species with recordings in < 5 hauls. For each haul, we determined species density (in number per m 2 ), using swept area (m 2 ) as the product of beam width and track distance.

Epifauna
For epifauna, we used the data collected in 2003 and 2004 by the internationally harmonized MAFCONS project (Managing Fisheries to Conserve Groundfish and Benthic Invertebrate Species Diversity) that sampled epifauna with a small beam trawl (2 m width, 5 mm mesh size). Five-min hauls were performed at 1 knot (covering ~ 300 m 2 ) at 283 stations (Fig. 3B). All caught organisms were determined to the lowest possible taxonomic level (see detailed methodology in Callaway et al., 2002). A checked and cleaned version of the dataset was used here, comprising a list of species abundances (in total number or weight per m 2 , depending on the species) per station (Robinson, Unpublished data). We removed stations outside our study area, species that were not epifauna, species that appeared in < 5 stations, and species without any weight or number recordings. Registrations with some unquantified presences were given the species-specific minimum observed weight/ number per m 2 .

Endobenthos
For infauna, we used the data sampled with a boxcore or grab (covering ~ 0.071 m 2 ) in the internationally harmonized North Sea Benthos Survey (NSBS) in 1986 (Kunitzer et al., 1992). The stations were regularly distributed across the entire North Sea (Fig. 3C). Samples were sieved over a 1 mm sieve and all organisms were subsequently identified to the lowest level possible, resulting in a dataset of species densities (see Heip et al., 1992 for detailed methodology). This dataset was obtained from the EMODnet Biology data portal (EMODnet Biology, 2018). Subsequent endobenthos surveys have been performed since the NSBS, but these could unfortunately not be merged in a single, sufficiently uniform dataset (results shown in Appendix A). To increase replicability of our study, we decided to only use the well-tested NSBS dataset, despite its age. Sampling stations outside our study area were removed, as were species that appeared in < 5 stations, and species without any recorded densities.

Biological clustering
Each biological dataset was analysed separately. We determined Bray-Curtis dissimilarities in community composition between sampling stations based on the fourth root of species densities (Reiss et al., 2010). Hierarchical clustering was then done using an average linkage sorting (Oksanen et al., 2019;Reiss et al., 2010), adopting a Flexible Dissimilarity Height (FDH-) approach. Based on visual inspection of the dendrograms, we performed a first cut-off at 0.50. Stations that were then in clusters with ≤ 5 stations were subsequently merged with larger clusters, as long as the dissimilarity height did not exceed 0.60. Remaining clusters with ≤ 5 stations were then removed as outliers. In Appendix D, we test the robustness of this approach by comparing it to a Fixed Number (FN-) approach, which identified a fixed number of clusters (here: 5, 10, 15) without considering dissimilarity height or cluster size.

Physiotopes
A total of 21 environmental variables were included in the physiotope clustering (Table 1; see Appendix B for detailed descriptions and figures of each individual variable). We additionally investigated the added value of fishing intensity as an explanatory, anthropogenic variable, which is further described in Appendix G.
All environmental variables were bilinearly interpolated to the highest resolution available (~180 × 180 m). Physiotope classification followed the methodology described in (Verfaillie et al., 2009). To avoid collinearity, all variables were summarized with a Principal Component Analysis (PCA) into environmental gradients. The first seven principal components (eigenvalue > 1) were kept, which together explained 77.5% of the observed variation. A hierarchical cluster analysis was performed on all grid cells, based on these seven main environmental gradients, using Wards method in the 'Rclusterpp' package (Linderman, 2013). Subsequent K-means partitioning used the hierarchical tree as starting points and were set to yield equal numbers of clusters as defined for biological assemblages, to facilitate the comparison between physiotopes and habitat maps.

Modelling habitat distributions
Habitat distributions were determined by applying Random Forest machine learning algorithm to the identified biological assemblages for the three faunal groups separately (Breiman, 2001). The Random Forest

Fig. 1.
Schematic overview of the methodology presented in this paper to create the final habitat (green boxes) and physiotope (yellow boxes) maps. Grey boxes depict input data. Purple letters refer to appendices where additional information can be found regarding the selection procedure of endobenthos data (A), the environmental variables used (B), the characteristics of resulting habitat maps (C), a comparison in clustering methodology (D), a comparison of interpolation methodology (E), the characteristics of resulting physiotopes (F), and the added value of fishing intensity as explanatory variable (G). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) classifier calculated the probability for a grid cell to belong to each assemblage, based on the main environmental gradients. This avoided collinearity of individual environmental variables and provided equal weight to the different gradients. Each model encompassed 500 internal runs. This resulted in three habitat maps with a resolution of ~ 180 × 180 m. Model accuracy was interpreted from 'Out-of-Bag' (OOB-)estimates and the median kappa value obtained from a cross-validation based on 100 iterations with half of the input data, using the 'rfUtilities' package (Evans and Murphy, 2018). We also determined spatial accuracy patterns. For this, a cross-validation was performed on five cluster-stratified folds after which the Shannon diversity index was determined for each grid cell. This index is 0 when grid cells were assigned to the same habitat (cluster) in all 5 predicted maps (so high accuracy), and increases to 1.6 when the cell was assigned to a different cluster in all 5 maps. The more different habitats a cell was assigned to in the 5 maps, and the more equal the proportional abundances of the five habitats were in the 5 runs, the more difficult it is to correctly predict which habitat the cell would be assigned to the next run. The Shannon entropy therefore quantifies the uncertainty (entropy or degree of surprise) associated with this prediction.
In Appendix E, the added value of environmental knowledge in the interpolation step is investigated. For this, we created alternative fullcoverage habitat maps, using the Nearest Neighbour method (NNapproach). This method applies Voronoi tessellation (Hijmans et al., 2017;Turner, 2019) on the point-wise biological clusters.

Anthropogenic pressures
We studied habitat associations of anthropogenic pressures by comparing the spatial distributions of habitats and anthropogenic activities. As anthropogenic variables we used international demersal fishing intensity, given as the average annual swept area ratio (SAR; year − 1 ) over 2010-2012 (Eigaard et al., 2017), and rasterized shapefiles of OWFs (EMODnet, 2019b) and Natura-2000 areas (EEA, 2018). Grid cells with a SAR > 1 were classified as 'fished', and average fishing intensity of the fished cells was calculated. In addition, habitat-specific extents were determined for the presence of OWFs and Natura-2000 areas. Finally, we compared fishing pressure and OWF presence inand outside of Natura-2000 areas.

Habitat maps
The datasets used for demersal fish, epifauna, and endobenthos included 1501, 192, and 149 stations and 73, 147, and 200 species, respectively. Hierarchical clustering resulted in 6, 13, and 10 assemblages for demersal fish, epifauna, and endobenthos, respectively, that were interpolated to full-scale habitat maps ( Fig. 4; Table 2; see Appendix C for detailed descriptions of the habitats and their characteristic species and environmental conditions). Herein, the FDH-approach was favoured over the FN-approach to determine the assemblages because it corrected for deviating stations that were unique clusters in the FN- The modelled fraction of mud (Stephens and Diesing, 2015;Stephens, 2015)  The modelled fraction of sand (Stephens and Diesing, 2015; Stephens, 2015) Gravel The modelled fraction of gravel (Stephens and Diesing, 2015; Stephens, 2015) Anthropogenic variable Fishing intensity* The average annual international fishing intensity over 2010-2012 (Eigaard et al., 2017) * The usage of this variable as an explanatory variable is further described in Appendix G.
approach (Appendix D). Secondly, the Random Forest models were favoured over the NN-approach (Appendix E) as their results matched better with our independent knowledge of the heterogeneity of the North Sea. Moreover, fishing intensity was not included as an explanatory variable as it did not alter the habitat distributions (Appendix G). The spatial extent of different habitat varied between faunal groups (DF: 0-53%; EF: 3-16%; EB: 3-23%; Table 2). However, general patterns could be observed among the three faunal groups (Fig. 4). Although the deeper waters northwest of the Dogger Bank were typified as one large demersal fish habitat (DF-1), they contained multiple epifauna and endobenthos habitats. These were located near the Norwegian trench (EF-2&3, EB-9), the Devil's Hole (EF-7, EB-3), the western coast (EF-1, EB-1&2) and the central North Sea (EF-4&8, EB-4) (Fig. 4). The shallower waters southeast of the Dogger Bank showed a similar pattern. Habitat DF-2 for demersal fish covered the entire area, whereas distinct epifauna and endobenthos habitats separated between the Brown Bank region (EF-11, EB-6), the Central Oyster Grounds (EF-10, EB-10), around the Dogger Bank (EF-5, 12&13, EB-5 & 7), and the   * Model accuracy details comprise the number of stations within an identified habitat (Nr. Stations), habitat-specific classification error (Class. Error) of the full model, and user (User acc) and prediction (Pred acc) accuracy estimates for the cross-validations. Overall Kappa and OOB-estimates are shown per model. ** Habitat details comprise the spatial coverage (%) and extent (in km 2 ). *** Anthropogenic pressures are displayed as the % of area that is fished (intensity > 1 year − 1 ), used for offshore wind farms, and designated as Natura-2000 area. Mean (±SD) fishing intensity (year − 1 ) is given for the fished (intensity > 1 year − 1 ) area. In addition, the habitat-specific percentage of Natura-2000 areas covered by the spatial extent of fishing activity and OWFs is given.
eastern coast (EF-9, EB-8). Interestingly, the two demersal fish habitats identified very locally near the outer Wash (DF-4, DF-5) did not have equivalent distinct epifauna or endobenthos habitats. The demersal fish Random Forest model had an overall Kappa of 0.82 and an OOB-estimate of 8%, showing a reliable prediction of the habitats (Table 2). Prediction accuracy was lowest for the smallest habitats (<0.01% coverage) DF-3 and DF-6, whereas the slightly larger sized (<1%) habitats DF-4 and DF-5 already showed accuracies of > 74% (Table 2). For epifauna and endobenthos, models showed reasonable, but lower overall Kappa values (0.59 and 0.58 respectively) and higher OOB-estimates (37% and 35% respectively; Table 2). Habitat extent, number of observations within the cluster, and habitat-specific accuracy were not clearly related (Table 2). Spatial patterns in accuracy differed between the three models (Fig. 5.) Uncertainty in demersal fish habitats was dominantly restricted to the Central Oyster Grounds and south of the Dogger Bank (Fig. 5A). For epifauna and endobenthos, higher Shannon diversity values were more common and widespread over the study area, representing lowered model consistency between the five cross-validations. Whereas the classifications of epifauna habitats were most consistent in the eastern and southern North Sea (Fig. 5B), endobenthos habitats were classified more consistently at the Dogger Bank, Central Oyster Grounds and the Brown Ridge Region (Fig. 5C).

Physiotopes
We applied hierarchical and K-means clustering on all grid cells, based on their score on the seven main environmental gradients resulting in high-resolution physiotope maps with 6, 13, and 10 clusters (P6, P10, P13; Fig. 4; see Appendix F for detailed environmental information of the physiotopes). To facilitate the comparison with the habitat maps, we visually matched the naming of individual physiotopes and habitats based on approximately shared locations. All three physiotope maps showed large ranges in spatial extent between physiotopes (P6: 1-41%; P10: 0-37%; P13: 0-23%; Table 3).

Habitat-specific anthropogenic pressures
All habitats were subjected to demersal fishing and/or offshore wind farms ( Table 2). The fished extent (in %) of each habitat varied between 9 and 65 % (Table 2), showing strong associations of demersal fishing with specific habitats. Preferred habitats (>40% fished extent) were dominantly located in the shallow waters southeast. Habitats with large OWF (>5%) and Natura-2000 area (>40%) extents were mostly located in the Brown Bank region and at the Dogger Bank. It should be noted, however, that especially the calculation of fished extent is probably less accurate for smaller habitats.

Table 3
Physiotope details as percentage spatial coverage, extent (in km 2 ), and anthropogenic pressures*. * Anthropogenic pressures are displayed as the % of area that is fished (intensity > 1 year − 1 ), used for offshore wind farms, and designated as Natura-2000 area. Mean (±SD) fishing intensity (year − 1 ) is given for the fished (intensity > 1 year − 1 ) area.

Discussion
Creating accurate, high-resolution marine habitat maps is challenging (ICES, 2019a), but essential for ecology-inclusive marine spatial planning (White et al., 2012). Here, we developed a novel method that combines point-wise biological data with grid-based environmental variables to create full coverage, high resolution habitat maps. We expand earlier work in the spatial modelling of species assemblages (see e.g. Ferrier & Guisan, 2006;Rubidge, Gale, & Curtis, 2016;Smoliński et al., 2017) by (i) the modelling of both physiotopes and observed assemblages for different faunal groups, (ii) our explorations in the usage of environmental gradients based on environmental and anthropogenic variables (Appendix G), and (iii) the direct comparison between habitat types and anthropogenic activities. As a result, our final habitat maps represent benthic assemblage distributions, congruent to spatial patterns in benthic assemblages as previously described by station-based studies (Kröncke et al., 2011;Reiss et al., 2010). Although several habitats showed close resemblance to physiotopes delineated from environmental variables, some essential differences were observed. Several physiotopes did not distinguish between multiple identified habitats, and one physiotope had no equivalent habitat, suggesting that important environmental variables or ecological processes might not have been captured in the delineation of physiotopes and habitats. Furthermore, we showed that separating maps for three faunal groupsdemersal fish, epifauna, and endobenthosyielded substantially different habitat maps, as expected based on their differences in ecology and mobility. By spatially comparing the habitats identified herein with anthropogenic pressures such as demersal fishing and OWFs, we confirmed that these activities have spatial associations with specific habitats (Grothe and Schnieders, 2011;. Hence, to sustainably manage anthropogenic pressures, habitat-specific impact assessments should be undertaken for the main faunal groups of benthic assemblages. Our findings provide an important new approach to marine habitat mapping, which will enable improvement of habitatspecific impact assessments for ecology-inclusive marine spatial planning. The presented methodology enabled the inclusion of distinct environmental variables which were relevant to specific communities and faunal groups (Hewitt et al., 2015;Lecours et al., 2015) and avoided subjective choices regarding included variables and exact classification boundaries. As a consequence, the resulting habitat maps predicted a substantially different spatial distribution of benthic assemblages compared to physiotopes maps, including the ones currently used for marine management and spatial planning (Ferrari et al., 2018;ICES, 2020;Rubidge et al., 2016). The Broad Habitat Types under the Marine Strategy Framework Directive (MSFD) are an example of such a currently used habitat map for management of the North Sea. The MSFD aims to achieve Good Environmental Status (GES) in all marine EU waters, for instance by the conservation of seafloor integrity (European Commission, 2008). Seafloor integrity is assessed by habitat-specific analyses of fishing pressures and sensitivity (ICES, 2019b;2020). This assessment uses the MSFD Broad Habitat Types (MSFD-BHT 1 ), which are dominantly based on water depth and sediment type (EMODnet, 2018;ICES, 2020). However, these MSFD-BHT show different habitat distributions compared to our maps, with higher habitat complexity than the demersal fish level, but lower complexity for both epifauna and endobenthos. Most remarkably, the clear distinction between the deeper northwest and shallower waters southeast of the Dogger Bank, which was dominant in all three presented habitat maps, is not captured by the MSFD-BHT map. Some resulting discrepancies between expected benthic assemblages and MSFD-BHT were already recognized during the assessment (ICES, 2019b), and more of such mismatches are likely given our results. These mismatches probably affect the required estimate of habitat sensitivity, which is parameterized by both habitat-and gearspecific depletion and recovery rates after a bottom trawl pass (Hiddink et al., 2017;2019;;Rijnsdorp et al., 2018). Especially habitatspecific recovery rates, based on longevity estimates of local endobenthos abundances will probably change if the boundaries of endobenthos habitats are adjusted (ICES, 2020;Rijnsdorp et al., 2018). In addition, a similar or integrated analysis for epifauna and demersal fish assemblages would allow for an assessment of the overall demersal ecosystem. Hence, our case-study showed that the current assessments of fishing impact on seafloor integrity probably fail to capture the true ecological status of benthic and demersal assemblages. Similar discrepancies between benthic assemblages and physiotopes have been observed in other areas than the North Sea as well (Ferrari et al., 2018;Rubidge et al., 2016).
Our methodology relies on the interpolation of biological assemblages with environmental variables, and its reliability is thus related to the quality of input data and the ability to apply a correct interpolation. Our input data, especially the endobenthos dataset, showed a spatiotemporal mismatch with the environmental data, increasing the uncertainty of the habitat maps (ICES, 2019a). Environmental data is often surveyed at or modelled to higher spatial and temporal resolution than what is available for biological data, which are generally snapshots in time. This can result in temporal mismatches, both in actual timing and represented time period. As a consequence, the correlations between benthic clusters and environmental gradients used for the interpolation have lower certainty. This can be observed in Fig. 5A, where the demersal fish Random Forest model shows high consistency between all cross-validations, probably because the input data has no spatiotemporal mismatch. The presented endobenthos habitat map has much lower certainty, caused by a large spatiotemporal mismatch in the input data in combination with (anthropogenic) changes in the North Sea since the initial sampling. Recent efforts to merge available endobenthos data for EMODnet demonstrated that the NSBS86 survey is still the leading survey in several regions of the North Sea (P.M.J. Herman, pers. comm.), with limited sampling stations in the northern part of our study area (Rees et al., 2007). And although a comparison of endobenthos composition in 1986 and 2000 showed relatively little difference in the large-scale community distributions (Kröncke et al., 2011), we strongly recommend for a new North Sea Benthos Survey that is internationally harmonized and covers the entire North Sea. Nevertheless, data availability will always dictate what can be used as input data, and therefore spatiotemporal problems will remain. International harmonization of national sampling surveys is therefore required to reduce spatiotemporal mismatches in marine habitat mapping (Appendix A).

Conclusions
We demonstrate a new approach to derive high-resolution, fullcoverage habitat maps representing benthic assemblages for different faunal groups. The detailed spatial information on benthic assemblages captured by our habitat maps are a first step to true ecology-inclusive marine spatial planning (Ferrari et al., 2018;Kaiser et al., 2016;Rubidge et al., 2016;White et al., 2012). Nonetheless, we find that the quality of the input data is key for accurate output, which demands international harmonization of existing sampling surveys. Our detailed overview in spatial distribution of both habitats and anthropogenic pressures identifies potential areas of conflicting interests and facilitate discussions on a fair balance between economic and ecological values. Assembly-specific knowledge based on species traits and speciesenvironment feedbacks, leading to assessments of resilience, sensitivity, and longevity, can now further improve habitat-specific impact assessments and subsequent management.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.