Predicting the Distribution of Indicator Taxa of Vulnerable Marine Ecosystems in the Arctic and Sub-arctic Waters of the Nordic Seas

In the deep waters of the Nordic Seas and adjacent areas, several benthic habitats such as cold-water coral reefs, coral gardens, and deep-sea sponge aggregations have been classified as vulnerable marine ecosystems (VMEs), due to their uniqueness, limited spatial extent, physical fragility, and slow recovery rate. In the last decade observations carried out by habitat mapping programmes in Norway, Iceland, and more recently in the Faroe Islands have substantially increased knowledge on the distribution of VMEs in the Nordic Seas. Nevertheless, large areas have not been explored due to the cost and logistics of obtaining observations in the deep-sea. Species distribution models can be used to predict the distribution of VMEs and their indicator species. Here we present the predicted distribution of 44 VME indicator taxa including 20 sponges, 17 cold-water corals, and 7 seapens in the Nordic Seas based on data compiled and models developed by the NovasArc project (2016–2018). Models for 44 VME indicator species were developed using the maximum entropy algorithm MaxEnt, using an extensive database compiled from habitat mapping surveys, by-catch data from bottom fish surveys, and records from reports and peer reviewed publications. Modeled distributions showed good agreement with observations. Niche overlap measures were used to identify seven groups and four subgroups of VME indicator taxa that co-occur. These were consistent with the species composition of known biotopes in the study area. A VME Index that combine the predictions for all VME indicator species was computed to identify particularly valuable and vulnerable ecosystems that should be targets of further exploration and conservation efforts. Such areas were identified at shelf break and slope off Iceland, the Faroe Islands, and central Norway, and the continental shelf off southern Greenland. The predicted distribution of VMEs in Arctic and sub-Arctic waters allows for the evaluation of interactions with fisheries and other anthropogenic activities and provides an important input for managers.


INTRODUCTION
In the deep sea, bottom trawling is the main source of anthropogenic impacts (Benn et al., 2010;Ramirez-Llodra et al., 2011). Of particular concern are the effects of bottom trawling on Vulnerable Marine Ecosytems (VMEs), ecosystems dominated by large epibenthic organisms (e.g., corals or sponges) which are likely to experience substantial alterations, and where recovery occurs very slowly (Wheeler et al., 2005;Clark et al., 2010;Williams et al., 2010;Buhl-Mortensen et al., 2013). In addition, VMEs are increasingly being threatened by pollution (Fisher et al., 2014) and by the effects of climate change including increasing water temperatures and ocean acidification (Guinotte et al., 2006;Levin and Le Bris, 2015).
The recognition that some deep-sea ecosystems are particularly susceptible to bottom trawling led the United Nations General Assembly (UNGA) to adopt resolutions 59/25, 61/105, and 64/72, calling for member states and regional fisheries management organisations (RFMOs) to identify areas beyond national jurisdiction (ABNJ) where VMEs occur, or are likely to occur, and to prevent significant negative impacts from damaging fishing practices. The Food and Agriculture Organisation of the United Nations (FAO) defined a set of criteria to identify VMEs, including their uniqueness or rarity, the functional significance of the habitat they form, structural complexity, fragility, and life history traits that make recovery difficult (e.g., slow growth, limited mobility; FAO, 2009). In general VMEs are identified by the presence of indicator taxa (e.g., stony corals, sponges), although merely detecting the presence of a VME element (species, habitats or features) is not sufficient to identify a VME (FAO, 2009). In the north-east Atlantic several benthic habitats have been classified by the Convention for the Protection of the Marine Environment of the North East Atlantic (OSPAR) and by the North East Atlantic Fisheries Commission (NEAFC) as vulnerable to human impacts (VMEs). In the sub-Arctic waters of the Nordic seas these habitats include seapen fields, cold-water coral reefs, coral gardens, and deep-sea sponge aggregations.
Vulnerable marine ecosystems are by definition susceptible to low levels of fishing pressure, and some types of spatial management have been shown to be effective to protect VMEs in the deep-sea (Clark and Dunn, 2012;Schlacher et al., 2014). These include confining the bottom trawling effort to already established footprints and establishing spatial closures to protect vulnerable species and habitats (Hourigan, 2014;Clark et al., 2015). In the North East Atlantic several closures have been established by NEAFC and OSPAR in areas beyond national jurisdiction (NEAFC, 2014). In addition, spatial closures to protect VMEs from fishery impacts have been established within the territorial waters of Norway (Fossåand Skjoldal, 2010), Iceland (Ólafsdóttir and Burgos, 2012), and the Faroe Islands (Anonymous, 2014).
A fundamental requirement for spatial management is detailed information on the distribution of VMEs. The optimal and non-destructive methodology to identify the occurrence of VMEs is through the use of underwater imagery, which allows an accurate description of the species composition and the abundance or density of organisms (e.g., Fabri et al., 2014;Anderson et al., 2016a;Beisiegel et al., 2017;Buhl-Mortensen et al., 2017). Nevertheless, the collection of underwater images in deep-sea environments is an expensive and complex endeavor, and therefore in most areas only a limited proportion of the seabed has been mapped visually. Often, the only available information is from the occurrence of VME indicator species in by-catch from fishery surveys and commercial trawls (Murillo et al., 2011;Durán Muñoz et al., 2012;Jørgensen et al., 2014). The lack of information on the distribution of VMEs in the deep sea is hampering the development and application of measures to protect these habitats from impacts of anthropogenic activities (Weaver et al., 2011).
Given the lack of extensive biological data on VMEs in most offshore environments and their presumed wide distribution, species distribution models (SDMs), also known as habitat suitability models or environmental niche models, are increasingly recognized as an effective way to obtain knowledge on the likely distribution of VMEs (Hourigan, 2014;Vierod et al., 2014;Clark et al., 2015). SDMs are models that predict the potential distribution of a species or a group of species in a given area using environmental variables as suitability predictors. Several studies have used SDMs to predict the distribution of VMEs (Howell et al., 2011(Howell et al., , 2016 and of VME indicator species (e.g., Rengstorf et al., 2013;Ross and Howell, 2013;Guinotte and Davies, 2014;Anderson et al., 2016b). The use of these models has been recommended for designing management plans to protect VMEs from fishing impacts (Ardron et al., 2014;Vierod et al., 2014). This includes the evaluation of the risk of fishing impacts (Penney and Guinotte, 2013) and the selection of areas for spatial closures (Lagasse et al., 2015;Rowden et al., 2019).
Within the Nordic Seas, the presence of VMEs has been documented by visual habitat mapping programmes in Norway and Iceland (Ólafsdóttir and Burgos, 2012;Buhl-Mortensen et al., 2015b), and records of VME indicator taxa have been obtained from commercial fisheries and scientific surveys. But predictive models of the distribution of VMEs and VME indicator taxa have not been developed in this area, except for Howell et al. (2016), who modeled the distribution of deep-sea sponge aggregations, and the predicted distribution of biotopes produced by the MAREANO programme in several regions in Norwegian waters (Elvenes et al., 2014;Buhl-Mortensen et al., 2015b;Gonzalez-Mirelis and Buhl-Mortensen, 2015). Here we present predictive models for a suite of indicator taxa of the most important VMEs in the Nordic Seas and adjacent areas. Fisheries Zone of the Faroe Islands. The study area also includes the entire NEAFC's Regulatory Area 2 (known as the "Banana Hole") and part of Regulatory Area 1. The study area is encompassed mostly in region I of OSPAR (e.g., the Arctic), but includes also small proportions of regions II, III, and V.
This area can be divided into three main basins separated by the northern extension of the mid Atlantic Ridge and the Greenland-Iceland-Scotland Ridge (GISR). The oceanography of the area is characterized by relatively warm surface water supplied from the south by the North Atlantic Drift, an extension of the Gulf Stream, overlying colder water masses (Norwegian Sea Deep Water, Arctic Intermediate Water) supplied from deepwater formation in Arctic areas (Buhl- Mortensen et al., 2015c). In coastal areas, the water is influenced by run-off from land. The seasonal variation is much less in the deeper waters than in the upper layers. Current velocities are controlled by the flow of the water masses and the tide, modified by the seabed topography. The GISR has a major impact on the distribution of water masses. The main pathway of water crossing this ridge is through the Wyville-Thomson Ridge between the Faroe Islands and Scotland. Here, the warm NAD passes into the Norwegian Sea above a sill of approximately 500 m. South of the Wyville-Thomson Ridge, the NAD water extends deeper and overlies a watermass characterized by water from the Mediterranean Ocean (the Mediterranean Outflow Water). The ridge system from Greenland to Scotland represents a major geographic barrier with great implications for distribution of marine species (Brix and Svavarsson, 2010;Dauvin et al., 2012;Omarsdottir et al., 2013).

Biological Data
Records of VME indicator species in the Nordic Seas were compiled from an extensive set of sources. Data were extracted from several databases in Norway, Iceland, and the Faroe Islands, including data from the Benthic Invertebrates of Icelandic waters (BIOICE) and the Marine Benthic Fauna of the Faroe Islands (BIOFAR) projects, and the Institute of Marine Research (IMR) coral database. In addition we used unpublished data from habitat mapping surveys by the Marine and Freshwater Research Institute (MFRI) in Iceland (Ólafsdóttir and Burgos, 2012) and the MAREANO project in Norway (Buhl- Mortensen et al., 2015b), and recent video observations carried out in the Faroe Islands (Buhl- Mortensen et al., 2019). We included by-catch data from the Joint Annual Norwegian-Russian Ecosystem Surveys in the Barents Sea (Jørgensen et al., 2015), and from the MFRI autumn surveys. We also extracted data from the ICES VME database (Morato et al., 2018), and the Ocean Biogeographic Information System (OBIS, Grassle, 2000). Finally, records were extracted from the literature including published data from the early expeditions, and the more recent work by Copley et al. (1996), Klitgaard and Tendal (2004), Mortensen et al. (1995Mortensen et al. ( , 2001, Cárdenas and Rapp (2015), and Hestetun et al. (2017). A complete list of the publications used to obtain records of VME indicator species can be found in Buhl- Mortensen et al. (2015b) and as an appendix in Buhl- Mortensen et al. (2019).

Selection of VMEs and Indicator Species
For the selection of the relevant vulnerable marine ecosystems and their indicator species for the Arctic and sub-Arctic area of this study, we considered: the VME classifications made by the Convention for the Protection of the Marine Environment of the North East Atlantic (OSPAR, 2010b); the list of VMEs compiled by the North East Atlantic Fisheries Commission (NEAFC, 2014) and the Northwest Atlantic Fisheries Organization (NAFO, Fuller et al., 2008); and the revised list of deep-water VMEs with characteristic taxa for ICES/NAFO waters by the ICES Workshop on Vulnerable Marine Ecosystem Database (WKVME, ICES, 2016). In addition we used recent experience gathered by national mapping projects in the study region. We selected 44 indicator taxa of 11 VME types and sub-types in the study area (Table 1). In total, 21 models were at the species level (e.g., Acanella arbusculla). The remaining 23 models were fitted at the genus level, either because most of the compiled record originated from video observations and the species could not be identified (e.g., Stryphnus sp.), or because there were several species of the same genus and the number of records was too low to allow modeling of individual species (e.g., Cladorhiza sp.). The VME types included in this study are the following:

Soft Bottom Sponge Aggregationse
In the Nordic Seas, demosponges of the order Tetractinellida form dense aggregations commonly known as "ostur" or "cheese bottom." These species can occur at depths between 150 and 1,700 m, on gravel and coarse-sand bottoms (Klitgaard and Tendal, 2004;Murillo et al., 2012;Maldonado et al., 2015). Two main types of ostur assemblages were recognized by Klitgaard and Tendal (2004) : the boreal "ostur" and the cold water "ostur." The boreal "ostur" is, according to Klitgaard and Tendal (2004), characterized by Geodia barretti, G. macandrewii, G. phlegraei, G. atlantica, Stelletta normaini, and Strypnhus ponderosus, although more recently it has been suggested that the latter species correspond to Strypnhus fortis (Cárdenas and Rapp, 2015;Maldonado et al., 2015). Boreal "ostur" assemblages were observed on some areas of the western Barents Sea, the Norwegian shelf (Kutti et al., 2013;Gonzalez-Mirelis and Buhl-Mortensen, 2015) and Faroese shelf (Klitgaard et al., 1997;Davison et al., 2019), and off southern Iceland (Klitgaard and Tendal, 2004). The cold water "ostur" is characterized by G. hentscheli, G. parva, and Stelletta raphidiophora, and it is found off northern Iceland, the Denmark Strait, off East Greenland, and north of Spitzbergen (Klitgaard and Tendal, 2004). The models of the distribution of species of the genus Geodia in the North Atlantic made by Howell et al. (2016) agreed with the observed distribution patterns. We fitted eight models of sponges considered indicators of soft bottom sponge aggregations: six models based on the Geodia species, and two models for sponges of the genera Stryphnus and Stelletta.

Hard Bottom Sponge Aggregations
A range of medium-to large-sized sponges occur on hard substrates including bedrock, lithified crust, and rocks. In the study area these include various axinellid sponges from the genera Axinella and Phakellia, and the demosponges Antho dichotoma and Mycale lingua. Off northern Norway, hard bottom demospongiae represents a single community (Gonzalez-Mirelis and Buhl- Mortensen, 2015). In addition to these four taxa, the ICES WGDEC (ICES, 2016) considered the family Tetillidae (genera Crainella and Tetilla), as well as sponges of the genera Polymastia and Tethya also to be indicators of hard bottom sponge aggregations, and these are frequently recorded in the Nordic seas (Buhl- Mortensen et al., 2012Mortensen et al., , 2015b. Models were fitted to seven indicator taxa of this VME.

Deep Arctic Sponge Aggregations
Several species of hexactinellid sponges are found in relatively high densities in deep cold (<0 • C) waters. One of the most common species in the Norwegian Sea is Caulophacus arcticus, which is generally found on hard bottoms at the lower part of the continental slope (Tendal and Barthel, 1993;Buhl-Mortensen et al., 2015b), and has been observed on the base of the Schultz Massif Seamount, at depths below 1,400 m (Roberts et al., 2018). Hexactinellid sponges of the genus Asconema can also constitute sponge grounds, although in restricted geographical settings (ICES, 2008). Asconema foliata has been observed on seamounts (Roberts et al., 2018), and is considered as a main habitat builder associated to cold water "ostur" habitats (Maldonado et al., 2015). In addition, poecilosclerid demosponges of the family Cladorhizidae become numerous at depths below 400 m, and at greater depths they constitute a large fraction of the sponge fauna (Hestetun et al., 2017). Several species of these carnivorous sponges of the genera Chrondocladia, Cladorhiza, and Lycopodina have been reported in the Nordic Seas (Hestetun et al., 2017). They are usually found in low densities, although aggregations of Chondrocladia grandis and Cladorhiza sp. have been observed off northern Iceland. Models were fitted to five indicator taxa of this VME.

Soft Bottom Coral Gardens
The term "coral garden" refers to relatively dense aggregations of colonies or individuals of one or several coral species (OSPAR, 2010a). They can be classified by substrate type (soft and hard bottoms) and the main representative taxa (ICES, 2016). Soft bottom coral gardens can be comprised by gorgonians of the families Isididae and Chrysogorgiidae, which can form dense aggregations on sandy mud (Buhl- Mortensen et al., 2015d). Among these, Isidella lofotensis is found almost exclusively off Norway (Buhl- Mortensen et al., 2015c), although it has been reported off east Greenland (Mayer and Piepenburg, 1996). Radicipes sp. aggregations have been observed off Norway only on the area known as the Bjørnøya slide, but it seems to be more widely distributed south of Iceland (Buhl- Mortensen et al., 2015c). In the warmer waters off southern Iceland the bamboo coral Acanella arbuscula is also relatively common. Soft-bottom coral gardens can also be comprised of solitary scleractinean corals of the genus Caryophyllia and Flabellum aggregated in relatively high densities forming what is known as "cup coral fields" (Baker et al., 2012;Buhl-Mortensen et al., 2015d). Models were produced for six indicator taxa of this VME.

Hard Bottom Coral Gardens
Hard-bottom coral gardens often occur in locations with strong currents. In the study area three of the subtypes from the ICES VME classification (ICES, 2016)

Reef-Forming Scleractineans
In Nordic waters only three species of scleractinean corals are reef building: Lophelia pertussa, Madrepora oculata, and Solenosmilia variabilis. Among them, L. pertusa is the most common and has been recorded frequently on the Norwegian shelf, around the Faroe Islands and off southern Iceland. M. oculata is less abundant, has a more limited framework-building capacity, and it often co-occurs with L. pertusa . In our study area, S. variabilis has been observed deep on the Reykjanes Ridge south of Iceland (Copley et al., 1996). Reef forming scleractineans do not always form reefs. For example, on vertical solid substrates coral debris cannot aggregate and reefs do not develop (Buhl- Mortensen et al., 2015d). In the North Atlantic reef-forming scleractineans can also form densely-packed "thickets, " as part of hard bottom coral gardens, or as isolated colonies Buhl-Mortensen, 2004, 2005;Davies et al., 2017). These growth forms are usually considered to represent different habitats. For example, the ICES VME classification distinguishes between cold-water coral reefs, non-reefal scleractineans, and colonial scleractineans on rocky outcrops (ICES, 2016). Because the growth form is seldom known or reported all observations of each taxon were grouped under a single VME type. Three models were fitted using presence data of each of the three reef-forming scleractineans.

Shallow Sea Pen Communities
Sea pen communities are usually defined as areas of bioturbated fine sediments with relatively high densities of sea pens. In OSPAR's list of threatened and/or declining habitats, this biotope is termed "sea-pens and burrowing megafauna communities" (Curd, 2010). This biotope is found in the relatively warm Atlantic water shallower than 700 m. The most common sea pen species are Funiculina quadrangularis, Virgularia mirabilis, Pennatula phosphorea, and Kophobelemnon stelliferum. Here we fitted four models based on records of these four species.

Deep-Sea Sea Pen Communities
The sea pen species Umbellula spp. and Anthoptilum spp. occur in deep waters (below 700 m) in an environment with colder temperatures and less anthropogenic activities than shallow water sea pens, and therefore should be regarded as a separate sea pen VME or at least a distinct sub-type (Buhl- Mortensen et al., 2019). High densities of Umbellula encrinus are found in deep waters north of Iceland and on the Norwegian slope, at depths below 800 m in the Norwegian Sea-Arctic Intermediate water. This large sea pen can reach a height of three meters. Off southern Iceland, sea pens of the genus Anthoptilum are also found in deep, albeit warmer waters. Models were fitted for both species.

Environmental Predictors
Bathymetry data for the study area was obtained from the General Bathymetric Chart of the Oceans (GEBCO version 20150318, Weatherall et al., 2015), a global relief model with a resolution of 30 arc-seconds. The data was projected using a Lamberts Equal Area projection centered at 69 • N and 4 • W and bilinearly interpolated to obtain a grid with a resolution of 500 m. All other environmental data sets were adjusted to match the same projection and resolution using bilinear interpolation. The seabed morphology was characterized following Lecours et al. (2017), using the following parameters derived from the 500 m bathymetry grid: local mean depth, slope, aspect (divided into northness and eastness), bathymetric position index (BPI), and vector roughness. The BPI indicates if a particular pixel forms part of a positive or negative feature of the surrounding terrain (Wilson et al., 2007). Vector roughness on the other hand measures the topographic surface roughness by quantifying the local variability in slope and aspect (McKean and Roering, 2004). Terrain analysis variables were calculated using a moving window of 3 cells, corresponding to a scale of 1,500 m, by fitting a bivariate quadratic polynomial to each window size using least squares. In addition BPI and vector roughness were also calculated using a moving window of 21 cells, corresponding to a scale of 10,500 m.
Temperature and salinity profiles for the study area were obtained from the NISE project (Norwegian Iceland Seas Experiment; Nilsen, 2008). Near-bottom temperature and salinity gridded fields were estimated following the methodology described in Jochumsen et al. (2016). The data going into the gridding corresponded to the deepest observation point. To avoid including shallow profiles, we only used observations obtained deeper than 80 % and within 80 m of the bottom depth. The gridding was then performed on a spatial resolution of 0.2 • longitudinal by 0.1 • latitudinal using an objective analysis and with an influence radius of 50 km. Minimum and maximum temperature, temperature difference (the difference between the lowest and highest temperature values), and mean salinity were calculated by interpolating along topography following Davis (1998) using a topography length scale of 300 km (Voet et al., 2010;Skagseth and Mork, 2012).
The aragonite saturation state for the study area was obtained from data provided by Jiang et al. (2015), and interpolated using a similar methodology to the temperature and salinity. Since the aragonite data is much more limited, the criteria for data to be included in the gridding was relaxed. We only required the observation depth to be within 80% of the bottom depth, and further we increased the influence radius to 200 km.
Monthly averages of mean net primary productivity (mg C m −2 day −1 , NPP) estimated from MODIS data using the carbon-based Production Model (CbPM) (Behrenfeld et al., 2005;Westberry et al., 2008) were obtained from the Ocean Productivity site 1 . Data was downloaded for the period 2006-2015 with a resolution of 5 arc min. Particulate organic carbon flux to the sea bottom (POC flux; g C m −2 year −1 ) was estimated from the bottom depth and the seasonal variation in NPP, defined as the ratio between the standard deviation and the mean of monthly NPP values (Lutz et al., 2002(Lutz et al., , 2007. Data on near-bottom average current speed and concentrations of nitrate, phosphate, and silicate were obtained from the Bio-ORACLE v2.0 database using the package "sdmpredictors" (Assis et al., 2018) in the R statistical environment (R Core Team, 2019), which provides layers of near-bottom physical and chemical parameters. Current velocity data (m.s −1 ) was produced by the Global Ocean Physics Reanalysis (ECMWF) using the OCEAN5 system at a native resolution of 0.25 • in the horizontal and 75 levels in the vertical, with the separation between vertical levels increasing with depth. Nutrient concentrations (mmol.m −3 ) were derived from the Global Ocean Biogeochemistry Non-assimilative Hindcast (PISCES). In both cases, data was obtained from the E.U. Copernicus Marine Service Information 2 , and statistically downscaled to a resolution of 5 arcmin using a kriging model (Assis et al., 2018).
Collinearity among environmental layers was explored by computing the Variance Inflation Factor (VIF) (Dormann et al., 2013). Variables with high collinearity were eliminated through a stepwise procedure in which the VIF was calculated for all variables, the variable with highest VIF was removed, and VIFs were recalculated until all variables had a VIF value lower than 10 (Naimi et al., 2014). This procedure selected maximum temperature, nitrate and phosphate as variables causing collinearity, therefore these parameters were not included as predictors. The correlation among remaining variables is shown in Figure 2. In the case of MaxEnt, it has been suggested that variables should be selected based on previous knowledge on the biology and ecology of the modeled species, but that stricter selection of variables is unlikely to improve models Elith et al., 2011). Therefore we did not attempted to select particular sets of predictor variables for individual taxa, and all models were constructed using all available variables.

Modeling Approach
In this study, rather than modeling the distribution of a VME by using the combined records of all indicator taxa (as done e.g., by Buhl- Mortensen et al., 2019), we opted to model each taxon individually. The rationale for this approach is that some of the VMEs include indicator taxa with different environmental requirements. For example, the corals Isidella lofotensis and Radicipes sp. are both considered indicators of soft-bottom coral gardens, but I. lofotensis is almost exclusively found on the Norwegian shelf, while Raicipes sp. is much more common off southern Iceland (Buhl- Mortensen et al., 2015c). A single model with records of both taxa will overestimate their distribution because a wider range of environmental settings would be FIGURE 2 | Correlation matrix among the environmental variables used as predictor for the species distribution models: Aragonite saturation state (Arag), broad-scale bathimetric position index (BPI BS), small-scale bathimetric position index (BPI SS), bottom depth (Depth), eastness (East), northness (North), mean net primary productivity (NPP), particulate organic carbon (POC), mean salinity (Sal), silica concentration (Si), bottom slope (Slope), near-bottom current speed (Speed), temperature difference (Tdiff), minimum temperature (Tmin), broad-scale vector roughness, and small-scale vector roughness. Colors and size of circles indicate correlation values, from 1 (blue) to -1 (red). considered suitable. To avoid this, we choose to (a) model individually each of the 44 taxon, (b) measure the similarity among the predicted distributions using the "I" similarity statistic (Warren et al., 2008), (c) use cluster analysis to identify groups of taxa with similar predicted distributions, and (d) combine the predictions of the taxa of each group using a stacked species distribution approach in order to obtain a predictive map of the distribution of each group. Finally, a VME index was computed to map the distribution and relative vulnerability of VMEs in the Nordic Seas.

Models of Individual Taxa
The distribution of VME indicator taxa was predicted using species distribution models (SDMs), which predict the geographic distribution of a species by identifying the combinations of environmental variables where they are observed to be more prevalent and then mapping that combination of variables into geographic space. We predicted the distribution of suitable habitats for VME indicator taxa using the maximum entropy algorithm MaxEnt (version 3.4.1). MaxEnt is a machine learning model that uses a presence-only approach to quantify the relationship between environmental variables at locations where a species has been observed vs. background locations in the study region (Phillips et al., 2006). MaxEnt uses transformations of the original environmental variables named "feature classes" (FC), namely the linear, product, quadratic, hinge, threshold, and categorical feature classes (Elith et al., 2011). Different combinations of feature classes allow the construction of very flexible models. By default, MaxEnt selects the number of feature classes based on the number of presence observations, increasing the number of feature classes with the number of records. To avoid overfitting, MaxEnt uses regularization, which penalizes the inclusion of parameters that produce small improvements in the model (Merow et al., 2013). Regularization is controlled by a parameter termed regularization multiplier (RM, default value = 1). Higher RM values reduce the flexibility in the relationships between species presence and environmental predictor variables. The performance of SDM models is sensible to model specifications (Elith et al., 2011;Merow et al., 2013;Warren et al., 2014). Recent studies have shown that the default MaxEnt options (i.e., the RM value and feature classes used) can produce models that perform poorly (Radosavljevic and Anderson, 2014). To select model settings that approximate optimal levels of model complexity, models were fitted to each VME indicator taxon using different combinations of feature classes and regularization multiplier values using the "ENMeval" (Muscarella et al., 2014) and "dismo" (Hijmans et al., 2017) packages in the R statistical environment (R Core Team, 2019). Each model was fitted using k-fold validation with five bins. To predict the distribution of each VME indicator taxon we selected the model with the highest average test AUC (area under the curve of the receiver operating characteristic plot), averaged across the 5-folds. This was the model with the best capability to successfully discriminate occurrence from background localities (Muscarella et al., 2014). In addition, we evaluated the Symmetric Extremal Dependence Index (SEDI, Wunderlich et al., 2019). SEDI is analogous to the widely used True Skill Statistic (TSS, Allouche et al., 2006) but better behaved in presence-background models because its error weighting reflects the low confidence in the pseudo-absence data, in particular in models with low prevalence and a high number of background points as the models in this study (Wunderlich et al., 2019). The SEDI for each model was estimated using a confusion matrix obtained after converting the model prediction into a binary presence-absence raster using the threshold that maximizes the sum of sensitivity and specificity (maxSSS), minimizing commission and omission errors (Liu et al., 2016). This threshold is commonly used to transform the output of SDMs into a binary output (Liu et al., 2005;Elith et al., 2006).
To evaluate if the models were overfitting we examined two metrics. The first is the difference between training and testing AUC, averaged across the 5 random folds (Warren and Seifert, 2011). This metric is expected to be high if models are overfitting the data (Muscarella et al., 2014). The second is the values of the 10% training omission rate (OR10). This is a threshold-dependent metric equivalent to the proportion of test localities with predicted suitability values lower than excluding the 10% of training localities with the lowest predicted suitability. Omission rates higher that the expected value of 10% typically indicate model overfitting (Muscarella et al., 2014;Radosavljevic and Anderson, 2014). Higher values are sometimes used to distinguish between degrees of overfitting. Here we followed Kivlin et al. (2017) by considering values below 0.2 as indicators of models with relatively low degrees of overfitting.
The model selected was used to indicate the distribution of VME indicator taxa in the study area based on the predicted occurrence of suitable habitat. Predictions were obtained for each cell in the same 500 m grid used for the environmental predictors. Model predictions were exported in the cloglog scale, which under specific conditions can be approximated to a probability of presence . In each model we computed the permutation importance of each predictor variable, which is the drop in AUC resulting from randomly permuting the values of the predictor variable on the training and background data sets and reevaluating the model (Phillips et al., 2006). Finally, models were examined visually to evaluate that high suitability areas corresponded to the locations of the majority of observations (Radosavljevic and Anderson, 2014).

Target Group Background
The spatial distribution of the available records of VME indicator species showed a strong sampling bias within the study area. In some areas the sampling intensity was high, in particular on areas of the Norwegian shelf and in the Barents Sea that have been mapped by the MAREANO programme (Buhl- Mortensen et al., 2015b), but also to some degree the Icelandic and Faroese shelves. In other areas like the eastern Greenland shelf and the deep basins the sampling effort has been very low or non-existent, or existing data is unavailable. Sampling bias can strongly influence the reliability of predictions of presence-background modeling. One way to reduce the effect of sampling bias in presence-only models is to use a set of background points with the same bias as the sampling effort (Phillips et al., 2009;Fourcade et al., 2014). To do this we modeled the sampling effort in the study area by fitting a kernel density estimate (KDE) to the locations of all indicator taxa compiled in the database, an approach known as target-group sampling (Elith et al., 2011;Merow et al., 2013). The KDE produced an estimation of the density of samples in each cell of the 500 m grid. These estimates were normalized so the sum of all cells was equal to one. These values were used to select 50,000 background points using the normalized kernel density values as a probability grid.

Niche Similarity
Niche overlap among all VME indicator taxa was estimated using the "I" similarity statistic (Warren et al., 2008), which ranges between 0 (no overlap) and 1 (niches are identical). To verify if VME classes currently used (e.g., ICES, 2016) consisted of indicator taxa with similar predicted distribution we carried out a cluster analysis using the Ward method, which defines groups by minimizing the within-group sum of squares (Legendre and Legendre, 1998). As a measure of dissimilarity between predicted distributions we used the complement of the "I" similarity statistic (1 -I). Groups of VME indicator taxa with similar distributions were identified from a dendrogram using a dissimilarity cutoff value of 0.4. This value was selected to produce groups roughly similar to known VME classifications.

Stacked Species Distribution Models
To map patterns of environmental suitability for VMEs, we merged the predicted distribution of VME indicator taxa to form stacked species distribution models (S-SDMs, Ferrier and Guisan, 2006;D'Amen et al., 2017;Wiltshire et al., 2018). S-SDMs is an approach that follows the "predict first, assemble latter" strategy in which the distribution of each individual taxon is modeled first, and then the predictions are combined (or "stacked") to produce a community prediction (D'Amen et al., 2017). This approach allows for the use of presence-only data where the presence records originated from different sources and where records of species of the same community or habitat are not usually collected at the same location. This is important when modeling deep-sea benthic megafauna in areas where the majority of records originated from fisheries by-catch and from scientific surveys using gear with relatively low sampling efficiency like dredges or bottom trawls (as opposed to underwater video surveys, which provide a more complete description of the benthic megafauna in a particular location). The S-SDMs approach also provides the flexibility of letting different environmental variables influence the distribution of individual species with distinct speciesenvironment relationships .
Here we produced S-SDMs for groups of VME indicator taxa with similar predicted habitat suitability, as defined by the cluster analysis of the "I" similarity statistic. SDMs were produced by averaging the predictions of individual VME indicator taxa (Calabrese et al., 2014;D'Amen et al., 2015;Wiltshire et al., 2018), and scaling the resulting average to a range of values between 0 and 1.
In the absence of knowledge on species prevalence, the output of presence-only models like MaxEnt is monotonically, but not proportionally, related to the relative probability of presence (Elith et al., 2011;Wiltshire et al., 2018), and cannot be interpreted as a measure of abundance or compared between species (Elith et al., 2011;Merow et al., 2013). As a result, the average suitability in each cell is likely to be related to the relative species richness in that location, although the relationship is not necessarily directly proportional (Aranda and Lobo, 2011;Guillera-Arroita et al., 2015;Wiltshire et al., 2018). Therefore, rather than attempting to predict species richness or community composition, we considered the average suitability as a tool to examine patterns of habitat suitability in each VME and to highlight areas of high average suitability that should be examined more closely and targeted for conservation measurements.

VME Index
To obtain a general overview of the distribution and relative vulnerability of the VMEs in the Nordic Seas, we computed a SDM-based VME Index analogous to the VME index developed for the ICES VME database (Morato et al., 2018). This index is a combination of a set of indicator scores, which quantifies in very broad terms the vulnerability of a taxa or group of taxa to anthropogenic impacts, and abundance scores based on abundance data of each taxon usually obtained from by-catch.
The indicator scores in Morato et al. (2018) were based on the vulnerability criteria defined by FAO (2009): uniqueness or rareness, functional significance, fragility, life-history that makes recovery difficult, and structural complexity. The degree to which each group fit each of the five criteria were scored using a scale between 1 (low) and 5 (high) by a group of experts (Morato et al., 2018). As the five indicators are considered to be approximately orthogonal, an indicator score was computed for each group using the quadratic mean (Morato et al., 2018). We used the indicator scores of the seven VME indicator groups present in our study area ( Table 2). Each of the 44 VME indicator taxon in our study were assigned to one of the seven VME indicator groups. Following Morato et al. (2018), sponges of the genera Asconema, Craniella, Geodia, Polymastia, Strypnhus, Tetilla, and Thenea were included in the "large sponges" group, while the remaining genera were considered as "generic sponges." The ICES VME index combines the indicator scores with an abundance score where by-catch weights are classified into a 1-5 scale based on current encounter threshold for live corals and live sponges established by NEAFC and the European Union (Morato et al., 2018). In an analogous manner, we scaled the predicted habitat suitability of each VME indicator taxon to a scale of 1-5 by computing the maxSSS threshold (Liu et al., 2016). Cells with values below this threshold received an abundance score of zero. Suitability values above the threshold were linearly scaled into a 1-5 scale, where 1 corresponds to the presence threshold and 5 corresponds to 1 (the maximum possible suitability value). The final SMD-based VME index was calculated by thresholding and scaling each of the predicted suitability for the 44 SDMs into a 1-5 scale, and multiplying by the corresponding VME indicator score. As opposed to Morato et al. (2018) we opted to average the VME index among the 44 indicator taxa, in order to highlight areas suitable for multiple VME indicator taxa.

Model Performance
All models had average test AUC values above 0.85, with 30 models having AUC values above 0.9, which indicates good model performances and capacity to distinguish between observation and background points ( Table 3). In addition, all models had relatively high SEDI values which indicate good model performance and a balance of omission and comission errors. In most cases (n = 33) the model with the highest AUC included all four feature classes (linear, quadratic, hinge, and product). This is the default behavior of MaxEnt, which selects the number of feature classes depending on the number of presence records. A total of 17 models had a regularization parameter (RM) of more than 1, which is the default in MaxEnt.
In general, models had a relatively low degree of overfitting, with an average OR10 (10% training omission rate) of 0.17. Only six models had OR10 values equal or above 0.25, which we took as indicator of high overfitting. The models with higher levels of overfitting were Anthothela grandiflora (OR10=0.50, Figure S6 Figure S10). In addition to these, the models of Geodia atlantica (Figure S2), Asconema sp. (Figure S9), and Duva florida (Figure S9) produced relatively high values of AUC diff (the mean difference between training and testing AUC) also suggesting overfitting.

Importance of Environmental Variables
There was high variability in the explanatory power of the 16 predictor variables, as measured by their permutation importance, on the 44 MaxEnt models ( Table 4). Minimum temperature was the most important predictor across all models, with an average permutation importance of 27.4% and explaining more than 10% of the variability in 40 of the 44 models. In individual models the importance of temperature reached up to 78.7% and it was particularly high for cold-water sponges (Caulophacus, Cladorrhiza, Chrondocladia Geodia phlegraei), but also for some cold-water corals (Lophelia pertusa, Madrepora oculata, Flabellum Gersemia), and gorgonians (Paragorgia arborea, Paramuricea, Primnoa resedaeformis, and Anthomastus). In addition, temperature difference had an average importance of 4.8%.
Depth was the second most important predictor, with an average importance of 20.2% and reaching up to 67.4%. Variables describing the morphology of the seafloor (small and large scale BPI and roughness, slope, northness, eastness) did not have high average importance by themselves (0.13-3.6%), but their averaged combined contribution was relatively large (11.86%) and reached up to 26.3%.
The combined effect of variables related to seawater chemistry (salinity, aragonite saturation state and silica concentration) was also relatively large (average importance 21.9%). Aragonite saturation state had an average importance of 6.5%, but had higher importance on some taxa including Caryophyllia (52.0%), Geodia phlegraei (20.21%), Stylasterids (32.9%), Drifa glomerata (14.6%), and Madrepora oculata (12.5%). In general salinity had a low explanatory power, with an average importance of 4.4%. For two Geodia species, G. parva and G. hentschelli salinity explained a large proportion of the variance.

Predicted Distributions and Niche Similarity
Predicted distributions of individual taxa are shown in Figures S1-S10. Pairwise niche similarity, measured by the "I" similarity statistic, ranged between 0.03 and 0.95, indicating a wide range of similarities among the predicted distribution of the VME indicator taxa. The cluster dendrogram shows seven groups (1-7) of VME indicator taxa at a dissimilarity level of 0.4 (Figure 3). Three of the groups can be divided into two subgroups each based on their dissimilarity value and the similarity of the spatial patterns of the predicted distribution. These groups represents VME indicator taxa with similar habitat suitability as predicted by the models. The following groups were identified:

Group 1
The first group includes the reef-forming corals Lophelia pertusa and Madrepora oculata, the gorgonians Parmuricea sp., Primnoa resedaeformis, and Pargorgia arborea, as well as corals of the family Stylasteridae. These species are often found in close proximity (Buhl- Mortensen et al., 2015c). Depending on local conditions, L. pertusa and M. oculata may be the dominant species and form coral reefs or coral thickets, or they can be found as isolated colonies forming part of hard-bottom coral gardens together with the other species in this group. This group of VME taxa is predicted to be distributed in narrow areas on the southern and western Icelandic shelf, around the Faroe Islands, off southern Greenland, and broadly on the central Norwegian shelf (Figure 4A).

Group 2
The second group includes three Geodia species (G. atlantica, G. macandrewi, and G. phlegraei) considered characteristic of the boreal "ostur" community (Klitgaard and Tendal, 2004), as well as the sponges Stelletta sp. and Strypnhus sp. Klitgaard and Tendal (2004) included Stelletta normani and Strypnhus ponderosus as part of the boreal "ostur, " while Stelletta raphidiophora was considered associated to the cold "ostur" community. Although we modeled Stelletta sp. and Strypnhus sp. at the genus level, given the geographic distribution of our samples it is likely that the majority of our records correspond to the species associated boreal "ostur" (Figure S2). High suitability for this group was observed off western Iceland, the Denmark strait, and the southern Greenlandic shelf, off the Faroe Islands, and in broad areas of the central and northern Norway shelf ( Figure 4B).

Group 3
The next group includes a number of sponge taxa usually associated to hard bottoms. Two subgroups can be recognized here. Subgroup 3A included Axinella sp., Phakellia sp., Antho (Antho) dichotoma, and Mycale (Mycale) lingua. Areas of high suitability for this group includes the central Norwegian shelf, the western and southern Icelandic shelf, and off the Faroe Islands ( Figure 4C). Subgroup 3B included Thetya sp., Polymastia sp. and Tetillidae, and also Geodia baretti, which is one of the Geodia species usually considered as part of the boreal "ostur" community (Klitgaard and Tendal, 2004). This group was associated to the colder waters of the Barents Sea, the Greenlandic shelf and off north Iceland, and to the shelf break off Norway ( Figure 4D).

Group 4
This group includes four indicator taxa of the shallow sea pen VME: Funiculina sp., Pennatula sp., Kophobelemnon sp., and Halipteris sp. Individual models for these taxa suggest that the four taxa have areas of high suitability south of the GISR, on the Skagerrak and in the North Sea, while only Funiculina sp. and Kophobelemnon sp. have relatively high suitability on the Norwegian shelf. The SSDM of these four taxa indicate some areas of high combined suitability on the southern Icelandic shelf and on the Skagerrak strait ( Figure 5A). Only predictors with permutation importance >5 in at least one model were included.

Group 5
This group incorporates a mixture of taxa with predicted suitability mostly south of the GISR including the gorgonians Radicipes sp., Anthotella grandiflora, and Anthomastus sp., the deep-water sea pen Anthoptilum sp., the cup coral Flabellum sp., and the reef-forming coral Solenosmilia variabilis (Figure 5B).

Group 6
This group consists of two subgroups of sponge taxa associated to cold waters. Subgroup 6A included Geodia parva and G. hentschelli, two species characteristic of the cold "ostur" assemblage (Klitgaard and Tendal, 2004). This group showed high suitability on the Greenlandic shelf ( Figure 5C). Subgroup 6B included three taxa of carnivorous sponges Caulophacus (C.) arcticus, Chladorhiza sp., and Lycopodina sp., and is predicted to be present mostly on the deep basins on the Norway and Greenland seas (Figure 5D).

Group 7
This group incorporates VME indicator species mostly associated to the continental slopes north of the GISR. Two subgroups can be identified in this group. Subgroup 7A included the deep-water sea pen Umbellula sp. and the carnivorous sponge Chrondocladia (C.) grandis, with predicted distribution mainly restricted to the deep continental slopes (Figure 6A). Subgroup 7B consisted of the VME cauliflower coral fields indicator species Duva florida, Drifa glomerata, and Gersemia sp., together with the sea pen Virgularia sp. and the carnivorous sponge Asconema sp. Its predicted distribution included the continental slope but also on broader areas of the shelves off Greenland, northern Iceland and the Faroe Islands, and on the Barents Sea ( Figure 6B). The two VME indicator taxa Isidella lofotensis and Caryophylla smithii had predicted distributions distinctly different from the others and were not part of any cluster. The gorgonian I. lofotensis is mainly restricted to Norwegian waters, but has been reported off northern Greenland (Mayer and Piepenburg, 1996;Buhl-Mortensen et al., 2015c), with a predicted distribution also including some areas around the Jan Mayen archipelago (Figure S10). The observed and predicted distributions of the cup coral C. (C.) smithii was restricted to some areas between the Shetland Islands and Norway, and within the North Sea.

VME Index
Values for the computed VME index ranged between 0 and 8.4 (Figure 7). Areas with high VME index values indicate locations where multiple VMEs are expected to be present, and/or where the VMEs present had high values in the indicator scores quantifying the vulnerability criteria defined by FAO (2009). The VME index suggest that even though some VME indicator taxa are predicted to have broad distributions, areas with high VME index values are more restricted. These areas include most of the Norwegian continental slope between 62 • N and 71 • N, coastal areas in the Barents Sea, the shelf break off the Faroe Islands and on the Faroe Bank, the southern and western Icelandic shelves, areas in the Reykjanes Ridge and the Kolbeinseyn Ridge, and some areas in the southern Greenlandic shelf and slope.

DISCUSSION
Here we presented the first comprehensive broad-scale modeling effort for VME indicator taxa on the Nordic Seas, including Icelandic and Faroese waters. Models of VME indicator taxa have been produced at more local scales within this area, in particular off Norway (Gonzalez-Mirelis and Buhl- Mortensen, 2015;Sundahl, 2017). Previous efforts in broad-scale models for cold-water corals (e.g., Davies and Guinotte, 2011;Yesson et al., FIGURE 5 | Predicted distribution of VMEs based on stacked species distribution model (SSMD) of (A) the sea pens Funiculina sp., Pennatula sp., Kophobelemnon sp., and Halipteris sp., (B) the gorgonians Anthotella grandiflora,/Radicipes/ sp. and Anthomastus sp., the deep-water sea pen Anthoptilum sp., the cup coral Flabellum sp., and the reef-forming coral Solenosmilia variabilis, (C) the sponges Geodia parva and G. hentschelli, and (D) the sponges Caulophacus (C.) arcticus, Chladorhiza sp., and Lycopodina sp. 2012) did not include our study area, and with the exception of Howell et al. (2016), who focused on Geodia sp., no broad scale models have been produced for VME indicator taxa for this area. The models presented in this study significantly expand the knowledge on the potential distribution of VME indicator taxa in Arctic and sub-Arctic waters, and provide a baseline for the evaluation of the presence of VMEs in the Nordic Seas.

Limitations of the Modeling Approach
SDMs are subjected to an array of sources of uncertainty (Vierod et al., 2014). Some degree of uncertainty is introduced by the observations of the taxa modeled. Positional uncertainty (Moudrý and Šímová, 2012;Naimi et al., 2014) and imperfect detection (Monk, 2014) are issues when using historical records, and data from fisheries by-catch and bottom trawl surveys. Sampling bias (Beck et al., 2014;Fourcade et al., 2014) is an important factor in our study, given that in some areas the number of observations is high (i.e., the Norwegian shelf, and to a lesser degree the Icelandic and Faroese shelves), while much lower in other locations (the Greenlandic shelf and the deep basins). Methods like the target group background approach used in this study can reduce the effects of sampling bias in the model prediction but cannot eliminate it completely. Additional sources of uncertainty arise from the environmental predictors. Oceanographic parameters are derived from databases like the World Ocean Atlas (Locarnini et al., 2013) or from physical ocean models (e.g., Logemann et al., 2013), often at relatively coarse resolutions. Estimates from both sources have their own uncertainty, but modeling approaches like MaxEnt do not incorporate this uncertainty in the predictions. Finally, the selection of the modeling approach itself can introduce biases in the prediction (Piechaud et al., 2015). In our study we utilized a single modeling framework, although some authors have suggested the use of ensemble models, averaging predictions from different models (Georgian et al., 2019). Given all the potential sources of the uncertainty, we echo Piechaud et al.  (2015) and urge caution when using the output of these models for management purposes. Modeled distributions should be used not as evidence of the presence of VMEs, but rather as a datadriven approach to identify areas where the presence of these habitats is likely (Hourigan, 2014).
In this study we have not produced uncertainty estimates for the predicted distribution of VMEs. Before predicted VME distributions can be used for management applications, it is necessary to quantify their uncertainty and to develop methods to incorporate the uncertainty in management decisions (Guisan et al., 2013). Although the internal uncertainty of MaxEnt models is difficult to quantify, bootstrap methods have been used to quantify some aspects of the uncertainty of MaxEnt predictions (Anderson et al., 2016b). For example, if planning tools like Zonation (Lehtomäki and Moilanen, 2013) or Marxan (Watts et al., 2017) would be used to prioritize areas for protection, it would be possible to prioritize locations with high conservation value (i.e., with high VME suitability and low uncertainty, Anderson et al., 2016b). Locations with confirmed VME presence, for example from underwater video surveys, would have no uncertainty and receive the highest priority for conservation. Uncertainty maps could also inform which areas should be targeted by future surveys, by highlighting locations where VMEs are predicted to occur and where the predictions are uncertain because of the lack of samples. An analysis of this type should follow this study.
We need to perform independent validation of the models (Elith et al., 2006;Davies and Guinotte, 2011) which can confirm model predictions (e.g., Rooper et al., 2018) or can highlight limitations of the predicted distributions. One form of validation is to compare model predictions with new observations. For example, recent observations in the Schultz Massif Seamount indicated the presence of Geodia parva, G. hentscheli, Stelleta sp., Caulophacus articus (Roberts et al., 2018), in agreement with our predicted suitability for these taxa. A comparison between predicted distributions and independent observations can provide useful information about the performance of the models. In an illustrative example, Anderson et al. (2016a) validated models for four reef-forming corals in the South Pacific Ocean using data from photographic surveys collected independently from the data used to fit the model. They found that the observed frequency of corals was much lower than predicted and that the correlation between observed and predicted coral distribution was not particularly high. The poor performance of the models was attributed to the low precision of the global bathymetry data, and to the lack of data on geomorphology and substrate data at the scale appropriate to the taxa modeled (Anderson et al., 2016a). These factors may be also relevant for the models in our study. An inspection of highresolution bathymetry derived from multibeam data available for the Norwegian shelf and some regions on the Icelandic shelf indicates that the GEBCO global bathymetry models are much less detailed and do not resolve small geomorphic features that may be important for the distribution of VMEs (Davies et al., 2009;Henry et al., 2010;Rengstorf et al., 2012). The lack of information describing substrates is also likely to affect the results of our models, as sediment composition is highly variable and is known to influence the distribution of epibenthic sessile organisms (Davies and Guinotte, 2011;Tracey et al., 2011). The effect of the lack of substrate data in our models can be illustrated by the fact that the cold-water coral model predicts high suitability in regions of the Skagerrak known to be dominated by soft sediments and where cold-water corals are usually not observed. The effect of the lack of sediment data is accentuated by the low resolution of the bathymetry model, because terrain variables derived from high-resolution bathymetry can play a better role serving as proxy variables for sediment composition (Dunn and Halpin, 2009). Given these factors, there is a need to produce SDMs at finer scales, incorporating high resolution bathymetry and sediment distribution data, if available.
When estimating the present distribution of VMEs and VME indicator taxa, depending on the goal of a study, it can be argued that the effect of historical fishing should be included (Ross et al., 2012). Some areas may have high predicted VME suitability, but if these areas are continuously being trawled, they may not have high concentrations of VME indicator species because of the cumulative effect of fishing-induced mortality. Penney and Guinotte (2013) suggested computing a "discounted suitability", where the suitability of each cell is reduced proportionally to the swept-area ratio. This method assumes that VME indicator species do not survive the impact of a single trawling event, and therefore in cells that are fished more than once per year the suitability is reduced to zero. This assumption can be adjusted to incorporate differences in the vulnerability of each VME to bottom trawling. An analysis of historical fishing patterns may highlight relatively pristine areas with high suitability for VME indicator taxa. These areas should be targeted for exploration and conservation.

Environmental Factors Influencing the Distribution of VME Indicator Species
Temperature is an important factor determining the distribution of cold-water corals (Davies and Guinotte, 2011;Yesson et al., 2012;Buhl-Mortensen et al., 2015c) and sponges (Klitgaard and Tendal, 2004;Howell et al., 2016). Given the strong bottom temperature gradients in the study area, it is no surprise that this parameter explained a large proportion of the predicted distribution patterns of the VME indicator taxa. In our study, minimum temperature was the most important factor with a permutation importance higher than 10% in 42 of the 44 taxa modeled. It was particularly high for species associated with cold waters including the sponges Geodia parva, Caulophacus (C.) arcticus, Cladorrhiza sp., Chrondocladia (C.) grandis, and the Neptheidae coral Gersemia sp. Temperature was also an important predictor for the scleractinean corals Lophelia pertusa, Madrepora oculata, and for taxa with predominantly southern distributions like Flabellum sp., Anthomastus sp., and Anthoptylum sp.
Next to temperature, depth was the next most predominant factor predicting the distribution of the VME indicator taxa in this study. Depth had a permutation importance higher than 10% in 33 of the 44 models. Similar to temperature, our study area included a wide depth range (Buhl- Mortensen et al., 2015c). Depth is not considered a direct explanatory variable, as hydrostatic pressure does not limit the distribution of VME indicator species. Instead it acts as a surrogate for other environmental parameters that are usually correlated with depth (Thresher et al., 2014;Howell et al., 2016). In our study there was some correlation between bottom depth and aragonite saturation state, particulate organic carbon and silica concentration. In addition, the environmental data sets used in this study had relatively low spatial resolution. In these conditions depth can act as a proxy explaining some proportion of the spatial patterns controlled by more direct explanatory variables.
Although the permutation importance of individual variables describing the morphology (i.e., terrain variables) of the seafloor was relatively low, their combined contribution was considerable in most of the models. This is consistent with previous modeling efforts for cold-water corals and sponges in which terrain variables have been important predictors (e.g., Rengstorf et al., 2013;Gullage et al., 2017;Rowden et al., 2017). It is well-known that the topography of the seafloor has a strong influence on the distribution patterns of filter feeders like cold-water corals and sponges, which are often associated with complex and elevated topographic features where locally accelerated currents increase the provision of food particles (Thiem et al., 2006;Duineveld et al., 2007;Navas et al., 2014) and influence the transport of larvae (Piepenburg and Müller, 2004;Dullo et al., 2008). The permutation importance of near-bottom current speed as a predictor variable was higher than 5% in 17 of the 44 models. These include four of the six models of sponges of the genus Geodia, and in the models of Stryphnus sp. and Steletta sp, taxa considered indicators of soft bottom sponge aggregations. Current speed had also an importance >5% in five of the seven models of sea pens (Pennatulacea). On the other hand, the importance of current speed was below 5% for most other taxa including reef-forming scleractinean corals, gorgonians, and sponges associated to hard bottoms. This is likely the result of the relatively low resolution of the near-bottom current speed data, which helped to predict the distribution of species with relatively broad distributions associated to soft sediments, but could not resolve more localized effects to predict the distribution of species associated with more complex terrain features. In the later cases, and similar to previous models for VME indicator taxa, the terrain descriptor variables act as proxies of nearbottom current speed. When the output of high-resolution oceanographic models is used as a predictor of the distribution of scleractinean corals like Lophelia pertussa or Madrepora oculata, near-bottom current speed do explain much of the observed patterns (Mohn et al., 2014;Bargain et al., 2018). In addition to acting as proxies of near-bottom current speed, terrain variables can also explain some of the spatial patterns associated to substrate type, which is another factor that strongly influence the distribution of benthic megafauna including corals and sponges (Gass and Roberts, 2006;Greathead et al., 2014;Baker et al., 2019). For example, hard bottom habitats are associated with complex topographies (Dunn and Halpin, 2009) where sediments are less likely to accumulate.
As expected, silicate was an important variable explaining the distribution of several sponge taxa. The permutation importance of silicate was higher than 10% in 12 of the 20 models for sponges. It was particularly high in the models of sponges associated to hard bottoms, including Axinellida sp., Phakellia sp., Antho (A.) dichotoma, and Mycale (M.) lingua. Silicate is required by sponges that form siliceous spicules, and its concentration may be a factor limiting the spatial distribution of sponge habitats (Leys et al., 2004;Howell et al., 2016). Silicate was also an important predictor for the distribution of some corals and sea pens, i.e., Isidella lophotensis, Anthomastus sp., Umbellula sp., and Pennatula sp. It is likely that for these taxa the concentration of silicate is not a direct limiting factor but rather is acting as an indicator of water masses. In the case of corals and sea pens, aragonite saturation had a lower importance than expected with only seven of 24 models having a permutation importance above 10%.
Geodia parva and G. hentscheli are the only two species for which salinity was an important predictor (permutation importances of 30.6 and 52.8% respectively). Both species are considered indicators of the cold "ostur" assemblage, with their distribution restricted to the Greenland shelf and Denmark strait.
Here it is likely that salinity is a proxy for water masses  rather than being a direct physiological constraint, acting as an indicator of the waters on the Greenland shelf characterized by the low temperatures and low salinity due to the input of glacial meltwater. This is supported by the fact that for G. hentscheli the permutation importance of temperature is lower, and the permutation importance of salinity is higher than for G. parva, but the combined importance of both parameters is similar for both species (70.1 and 60.9%, respectively).

Spatial Distribution of Predicted Taxa and Taxa Groups
As expected, there was a good agreement between the observed locations of VME indicator species and the areas of predicted high suitability, in particular for taxa distributed along a narrow range of environmental variables like Lophelia pertusa and the other taxa in cluster group 1 (Figure S1), Umbellula sp., and for taxa with rather limited geographic distributions like Isidella sp. and Caryophila. For other taxa the prediction appears less precise, the models tended to predict broad areas of high suitability. Some examples are the models for Geodia phlegraei, Anthothela gradiflora, Anthomastus sp., Flabellum sp., and Lycopodina sp. This could be because the taxa modeled are generalist species with a broad distribution, or could be an artifact due to low sample sizes or the lack of environmental predictors that limit their distribution.
Our predicted distributions for Lophelia pertusa, Paragorgia arborea, and Primnoa resedaeformis are very similar to those obtained by Sundahl (2017) in Norwegian waters using mostly data from the MAREANO programme (Buhl- Mortensen et al., 2015a,b). This is reassuring, given that Sundahl (2017) produced models at a higher spatial resolution (176 m) and included environmental predictors that were not available for our models, including sediment type and the output of a high-resolution ocean model. In addition, our areas of predicted high suitability for Lophelia pertusa in the territorial waters of the United Kingdom and Ireland are similar to the areas identified by Ross and Howell (2013) as having medium and high probability of reef presence.
In general, the predicted distributions for sponges of the genus Geodia were comparable to those predicted by Howell et al. (2016) in our study area. In particular, the predicted presence of G. mandrewi were very similar to the areas of high suitability predicted by our model. Our predictions for G. parva and G. hentschelli were also similar,  Elvenes et al. (2014) and Gonzalez-Mirelis and Buhl-Mortensen (2015) with associated VME indicator taxa, and the cluster dendrogram groups in which those taxa were assigned in our study.  Howell et al. (2016), whose predicted areas with presence included the deep basins of the Norwegian Sea and off southern Iceland where our models predicted low suitability. The cluster analysis of the "I" similarity statistic computed among the predicted distributions of VME indicator taxa produced 7 groups and 6 subgroups. These groups are not biotopes or communities in the strict sense, as they are not the result of direct observations of species living in close proximity, but rather are groups of taxa with similar predicted distributions. Predictions from the broad-scale models in our study are not expected to agree with the species composition and distribution of biotopes identified from high-resolution local predictions, but reveal similarities at large scales relevant for the management of large marine ecosystems. Nevertheless, in several instances the groups identified were analogous to known VMEs in the study area, and in some cases to biotopes identified on the Norwegian shelf by the MAREANO programme (Table 5).
Group 1 included the gorgonians Paramuricea sp., Paragorgia arborea, and Primnoa resadaeformis, which are indicators of hard-bottom coral gardens (ICES, 2016) and were identified as an homogeneous community in the Norwegian shelf (Gonzalez-Mirelis and Buhl- Mortensen, 2015). Group 1 also included the scleractinean corals Lophelia pertusa and Madrepora oculata, which are the most important reef-forming corals in the study area. These two species are indicator taxa of cold-water coral reefs, but also of two types of hard-bottom coral gardens: colonial scleractineans on rocky outcrops, and non-reefal scleractinean aggregations (ICES, 2016). In our study we could not distinguish between L. pertusa and M. oculata records originating from different VMEs, and therefore we consider them indicators of a single VME termed reefforming scleractineans. In addition, in the study area the taxa in this group are often found in close proximity (Buhl- Mortensen et al., 2015c), and therefore their predicted broadscale distribution was similar. At local scales, it is likely that there are differences in the distribution patterns among the taxa in the groups identified in this study. For example, Elvenes et al. (2014) identified ten biotopes on the Norwegian shelf off the Lofoten and Vesterålen archipelagos (Elvenes et al., 2014), seven of which included VME indicator taxa. Lophelia and Primnoa were assigned to the same biotope (biotope class 10, Table 5), while Paragorgia was included in biotope class 8. These differences are expected when comparing predicted distributions at relatively large scales with biotopes identified at more local scales.
Five of the groups identified in our analysis were dominated by sponges and corresponded well to sponge VMEs within the study area. Our groups 2 and 6A include indicator species of softbottom sponge aggregations (ICES, 2016) and are analogous to the boreal and cold "ostur" assemblages described by Klitgaard and Tendal (2004), respectively. The stacked distribution models for groups 2 and 6a suggest that "ostur" habitats are widely distributed on the continental shelves and slopes of eastern Greenland, Denmark Strait, Iceland, the Faroe Islands, the Norwegian shelf, and the Bartents Sea. This agrees with previous studies by Klitgaard and Tendal (2004), Christiansen (2010), Howell et al. (2016). Geodia aggregations have been observed at depths between 150 and 1,700 m (Maldonado et al., 2015), which correspond well to the depth ranges of predicted high suitability for group 2 (boreal "ostur"). Areas of predicted high suitability for group 6a (cold "ostur") includes areas in the Norwegian and Greenland seas at depths below 1,500 m. The sampling effort at these depths is very low and no Geodia aggregations have been previously reported. Group 3 included indicator taxa for hard-bottom sponge aggregations (ICES, 2016). The fact that these taxa form a distinct group from sponges associated to soft-bottoms supports the notion that deep-sea sponge aggregations form distinct habitats depending on their bottom type preferences as proposed by Gonzalez-Mirelis and Buhl- Mortensen (2015).
The list of indicator taxa for ICES/NAFO waters (ICES, 2016) also includes Caulophacus articus as an indicator for soft-sponge aggregations. In an analysis of community structure, Gonzalez-Mirelis and Buhl- Mortensen (2015) concluded that Hexactinellid and other sponges associated to deep, cold waters constitute a distinct habitat separated from soft-bottom and hard-bottom sponge aggregations. Our analysis placed C. articus, together with Cladorrhiza sp. and Lycopodina sp. into group 6B, which can be considered analogous to the deep sea sponges habitat defined by Gonzalez-Mirelis and Buhl- Mortensen (2015). This supports the need to define a distinct VME type for taxa in deep, cold waters, as suggested by Buhl- Mortensen et al. (2019). Caulophacus was associated to a distinct biotope (biotope class 2) by Elvenes et al. (2014), and our predicted distribution for Caulophacus sp. shows high suitability at depths below 1,500 m, which agrees with the predicted distribution of biotope 2 (Elvenes et al., 2014).
Group 4 included four sea pens considered indicators of the VME type "sea-pen fields" in the list of indicator taxa for ICES/NAFO waters (ICES, 2016): Funiculina sp., Pennatula sp., Kophobelemnon sp. and Haliperis sp., and it is analogous to the biotope "sea pen and burrowing megafauna communities" identified by Gonzalez-Mirelis and Buhl-Mortensen (2015) on the Norwegian shelf. The ICES/NAFO list also includes deepsea sea pens Umbellula sp. and Anthoptilum sp., which our analysis did not included in group 4 given the differences in their predicted distributions. Similarly, Gonzalez-Mirelis and Buhl- Mortensen (2015) concluded that the sea pens Funiculina cuadrangularis, Kophobelemnon stelliferum, Pennatula sp., and Virgularia spp. formed an homogeneous community on the Norwegian shelf distinct from the Umbellula stands (e.g., areas with high densities of Umbellula sp). This supports the definition of two separate VME types which include shallow and deepsea sea pen species (Buhl- Mortensen et al., 2019). Gonzalez-Mirelis and Buhl- Mortensen (2015) predicted the distribution of Umbellula sp. stands in two regions in the continental slope at depths below 500m roughly off the Lofoten and Vesterålen archipelagos, and in the Eggakanten north of 71 • N. Our model for this taxa predicts very high suitability values in similar locations ( Figure S8). In our analysis Umbellula spp. clustered together with the sponge Chrondocladia (C.) grandis, as both species have a distinct predicted distribution in the slopes north of the GISR. A biotope comprised by both species was also identified on the slopes of the Norwegian shelf (biotope class 3, Elvenes et al., 2014). Our models for both species and our stacked model for group 7a predicted high suitability along the continental slope between 800 and 1,500 m, which agrees well with the predicted distribution biotope class 3.
As opposed to most other groups, groups 5 is formed by taxa indicators of different VMEs characterized by a predicted distribution mostly in the deep basins south of the GISR. Among the taxa in this group are indicators of soft bottom coral gardens, including the gorgonians Radicipes sp., Acanella arbuscula, and the cup coral Flabellum, and indicators of hard bottom coral gardens including the gorgonians Anthomastus sp., Anthotella grandiflora. A. arbuscula has been reported to co-occur with Flabellum alabastrum (Buhl- Mortensen and Buhl-Mortensen, 2018). Our model for Radicipes sp. correctly predicted high suitability in the Bjørnøya slide area, which is the only area in Norwegian waters where Radicipes meadows have been observed (Buhl- Mortensen et al., 2015c;Gonzalez-Mirelis and Buhl-Mortensen, 2015, Figure S6). The group also includes Solenosmilia variabilis, a reef-forming scleractinean coral than in our study area was recorded south of the GISR, and the deep-sea sea pen Anthoptilum sp. This group includes taxa reported in the Mid-Atlantic ridge at a relatively wide depth range (800-2,400 m, Mortensen et al., 2008).
Group 7B includes mostly corals of the family Nephtheidae family that are indicator taxa of the VME subtype cauliflower coral fields (ICES, 2016;Buhl-Mortensen et al., 2019). These corals are found over a wide range of substrates including semiconsolidated mudstone (Buhl-Mortensen and Buhl- Mortensen, 2018), sometimes in relatively high densities as observed off the Westfjords in Iceland. On the Norwegian shelf, Elvenes et al. (2014) identified a biotope (biotope 8) that included Gersemia and Drifa, but also Paragorgia ( Table 5).
Our VME Index provides a summary of the distribution of all 44 taxa, giving more weight to taxa that are considered to more closely fulfill the FAO criteria for VME identification (FAO, 2009). Areas with high VME Index include much of the shelf break and slope off Norway and the Barent Sea, Iceland, and the Faroe Islands, the shelf off southern Greenland, and the areas in the Reykjanes Ridge and the Kolbeinseyn Ridge. Several of these areas are being targeted by the MAREANO programme in Norway, and by the habitat mapping efforts by the Marine and Freshwater Research Institute (MFRI) in Iceland. Recent video observations carried out in the Faroe Islands are providing new information on the distribution of VME indicator species (Buhl- Mortensen et al., 2019). Nevertheless, there are vast areas where the number of observations is very low and should be the target of new research efforts. The continental shelf off southern Greenland is of particular interest because the predicted suitability of several VME indicator taxa is high in this region, and the area has not been subjected to intense fishing effort which increases the probability of finding pristine VMEs. Recent mapping efforts have been producing information off the west coast off Greenland (Yesson et al., 2016), but little is know off the eastern coast. Also there is a lack of data on the fauna on the deep basins of the Norwegian and Greenland sea.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
LB-M, PB-M, JB, SÓ, SR, and PS conceived the study. LB-M secured funding. LB-M, PB-M, SÓ, PS, and SR compiled the database with records of VME indicator taxa. ØR provided data on temperature and salinity and performed spatial interpolations. JB compiled the remaining environmental predictors and carried out species distribution models and other analysis. The manuscript was written by JB, and was contributed to and edited by LB-M, PB-M, SÓ, SR, and ØR.

FUNDING
This article is a product of the Nordic Project On Vulnerable Marine Ecosystems And Anthropogenic Activities In Arctic And Sub-Arctic Waters (NovasArc, 2016(NovasArc, -2018, a collaboration between the Institute of Marine Research (Norway), the Marine and Freshwater Research Institute (Iceland) and the Faroe Marine Research Institute (Faroe Islands). NovasArc was supported by the Marine Group (HAV) and the Working Group for Fisheries (AG-Fisk) of the Nordic Counsel of Ministers.