Pan-Arctic suitable habitat model for Greenland halibut


 Deep-sea marine fishes support important fisheries but estimates of their distributions are often incomplete as the data behind them may reflect fishing practices, access rights, or political boundaries, rather than actual geographic distributions. We use a simple suitable habitat model based on bottom depth, temperature, and salinity to estimate the potential distribution of Greenland halibut (Reinhardtius hippoglossoides). A large presence-only dataset is examined using multivariate kernel densities to define environmental envelopes, which we link to spatial distribution using a pan-Arctic oceanographic model. Occurrences generally fit the model well, although there were gaps in the predicted circum-Arctic distribution likely due to limited survey activity in many of the ice-covered seas around the Arctic Ocean. Bottom temperature and depth were major factors defining model fit to observations, but other factors, such as ecosystem interactions and larval drift could also influence distribution. Model predictions can be tested by increasing sampling effort in poorly explored regions and by studying the connectivity of putative populations. While abundances of Greenland halibut in the High Arctic are currently low, some areas are predicted to be suitable habitat for this species, suggesting that on-going sea-ice melt may lead to fisheries expansion into new areas.


Introduction
Understanding the distribution of fishes is important for the efficient and sustainable management of fisheries. Targeted fish species may populate habitats outside known fishing grounds that are difficult to access, and consequently, estimates of their population structure and spatial distribution can be incomplete.
The data to generate distribution estimates are often collected directly from fisheries and national research surveys. Such results may reflect fishing practices, access rights, and political boundaries rather than the actual geographic distribution of a species (Kerr et al., 2017). Compiled datasets can help, but these data are limited to areas with frequent fishing. Modelling the V C International Council for the Exploration of the Sea 2021. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
distribution of a species based on constraints set by oceanography and ecology can be useful to identify regions where little or no data are available but where fish might be present (Planque et al., 2011).
Identifying the factors controlling the distribution of a species has been a central question in ecology and the focus of a century of research (Grinnell, 1917;MacArthur, 1972). It has been suggested that distributions of marine fishes are controlled by a combination of geographical attachment, environmental conditions, density-dependent habitat selection, spatial dependency, demographic population structure, species interactions, and population memory (Planque et al., 2011). Numerical models of species distributions have been the object of intensive ecological research in the last few decades (Guisan and Thuiller, 2005;Elith and Leathwick, 2009;Zurell et al., 2020). Most modelling efforts have focused on climate envelope models that relate to environmental niche theory (Hutchinson, 1957;Soberó n and Nakamura, 2009). Bioclimatic envelope models, also known as niche-based, habitat suitability, and species distribution models, use associations between aspects of climate and species' occurrences to estimate the conditions that are suitable to maintain viable populations (Araú jo and Peterson, 2012). Typically, a multidimensional set of environmental variables limiting the distribution of a species is first defined from empirical data. A measure of distribution in the environmental space (e.g., boundaries, probability of occurrence, or density) is then calculated using models (Lauria et al., 2015), kernel-density (Broennimann et al., 2012), or principal components (Soberó n and Nakamura, 2009) to form a multivariate envelope, also called environmental or habitat space. Finally, the environmental envelope is projected to geographic space typically using a physical model, an oceanographic model in the case of marine fishes, providing the advantage that the empirically defined environmental space and the modelled geographic space are independent of each other. Many such approaches, however, are limited to variables derived from the physical model with sufficient empirical data and rarely account for species dispersal or ecosystem interactions, leading to speculative models for which "suitable habitat model" (SHM) is arguably a more accurate term than "species distribution model".
For marine fishes, SHMs have been used to investigate factors controlling species distributions at seasonal and interannual time scales (Planque et al., 2007;Loots et al., 2010), investigate possible climate-induced shifts in species distributions (Wisz et al., 2015;Beaumont et al., 2016;Morato et al., 2020) as well as to project shifts in communities and biodiversity (Jones and Cheung, 2015). The nature of the data used to construct SHMs can vary greatly, from occurrence only, to presence/absence, to biomass or abundance, and so can the degree of complexity of the models. The challenge lies in finding a model that balances fit and predictive power (pp. 37 in Hastie et al., 2017) while adequately capturing the underlying hypothesis (Araú jo and Peterson, 2012) and ensuring model transferability (Yates et al., 2018) given the constraints set by data availability. Failure to meet these challenges can lead to poor model performances (see, e.g. the case of Barents Sea cod projections discussed in Ingvaldsen et al., 2015). Despite their shortcomings, SHMs have proven valuable tools to conceptualize our knowledge of real-world species distributions. Such models provide a useful framework to predict where and when habitats may be suitable, even if the species has not been observed at those specific locations.
Greenland halibut (Reinhardtius hippoglossoides; also known as Greenland turbot or black halibut) is a commercially important flatfish supporting demersal fisheries in northern regions of both the Atlantic and Pacific Oceans (Smidt, 1969;Godø and Haug, 1989;Bowering and Brodie, 1995) with some records along the continental slope of the Arctic Ocean ( Figure 1; Hedges et al., 2017;Majewski et al., 2017;Orlova et al., 2017;Mecklenburg et al., 2018). Although considered a demersal species, Greenland halibut largely feed on pelagic prey (Hovde et al., 2002;Vollen and Albert, 2008;Dwyer et al., 2010;Giraldo et al., 2018) and individuals may be in the pelagic realm for 10-20% of the time (Albert et al., 2011). Outside the feeding trips, the adults are demersal, inhabiting a depth range between 200 and2000 m (Bullough et al., 1998;Bowering and Nedreaas, 2000;Morgan et al., 2013), while early life stages are exclusively pelagic (Domínguez-Petit et al., 2013). After spawning, which takes place in the mesopelagic layer, the eggs drift with ocean currents developing to larvae within 1-2 months (Domínguez-Petit et al., 2013). A few days before the eggs hatch, their buoyancy changes, and the larvae likely hatch in the euphotic zone where they continue to drift with ocean currents for up to four months (Sohn et al., 2010). The total drift time is dependent on the temperature and can vary between 3 and 8 months allowing dispersion with ocean currents over large distances (Sohn et al., 2010;Domínguez-Petit et al., 2013). After the pelagic larval stage, juveniles become demersal, moving deeper as they grow (Godø and Haug, 1989;Sohn et al., 2010). Greenland halibut typically inhabits a temperature range between 0 and 4 but have been recorded from temperatures between À1.9 and 10 (Peklova et al., 2012;Morgan et al., 2013;Wheeland and Morgan, 2020). Greenland halibut observations are generally rare from regions covered in sea-ice throughout much of the year. This may be due to the ecology of the species but might also reflect low fishing effort in these high Arctic waters (Christiansen and Reist, 2013;Chernova, 2017;Majewski et al., 2017). The distribution may also reflect spawning locations and their relations to prevailing ocean currents (Å dlandsvik et al., 2004;Sohn et al., 2010). Currently, three major spawning regions have been reported in the literature ( Figure 1): (i) the Davis Strait between West Greenland and Arctic Canada (reviewed in Bowering and Brodie, 1995), (ii) north of the Aleutian Islands in the Bering Sea (Sohn et al., 2010;Bryan et al., 2018), and (iii) the continental slope along the Norwegian coast towards Svalbard in the Barents Sea (Godø and Haug, 1989;Albert et al., 2001b). Greenland halibut also appear to spawn outside the major spawning regions and winter spawning season (Bowering and Brodie, 1995;Albert et al., 2001b). Although these localized spawning events may produce fewer recruits than the major spawning regions, multiple small spawning areas may also influence the distribution and population structure of the species and can remain undocumented in rarely surveyed regions. Localized spawning has been conjectured to occur southwest of Iceland (Magnusson, 1977), along the southeast coast of Greenland (Gundersen et al., 2013), in the West Greenland fjords (Riget and Boje, 1989), along with the continental slope of Labrador and Newfoundland (Bowering, 1982;Junquera and Zamarro, 1992;Bowering and Brodie, 1995 and the references to grey literature therein), in the Gulf of St. Lawrence (Ouellet et al., 2010;Domínguez-Petit et al., 2013), and the Beaufort Sea (Chiperzak et al., 1995).
Here, we use a simple SHM approach together with an occurrence dataset of Greenland halibut compiled from the databases of institutions working with the assessment of the species in the Arctic. We use bottom depth, bottom temperature, and salinity in our multivariate environmental envelopes to study how these factors may influence Greenland halibut distribution and to examine missing limiting factors. Further, we use the model to identify potential locations of under-sampling and to discuss the potential connectivity of populations suggested by the model in a pan-Arctic perspective. The approach is applied to Greenland halibut, but may be used for other species for which relevant environmental data and ecological knowledge are available.

Material and methods
We used a non-parametric, non-linear multivariate approach to define suitable habitat following the environmental envelope or space-based concept (Soberó n and Nakamura, 2009; Araú jo and Peterson, 2012), which we call "habitat spaces". The concept relied on the assumption that bottom depth, temperature, and salinity are together limiting the distribution of Greenland halibut (Morgan et al., 2013;Wheeland and Morgan, 2020). The modelling approach consisted of distinct steps, each of which is explained in Figure 2. First, a cloud of observations of individual fish was plotted in the bottom temperature-logarithm transformed depth space (TD-space). The data for suitable bottom condition estimation was acquired from five sources (Data for habitat spaces). Next, the suitable bottom conditions were limited Figure 2. The suitable habitat model. First, available temperature-depth (TD) observations were plotted on a scatter plot (1). Next, the TDspace was estimated from the kernel-density of binned TD data (2). Salinity limits were added to the TD-space resulting in a cylindrical TDSspace (3). NEMO temperature, depth, and salinity estimates for the Arctic were constrained using the TDS-spaces (4-5) resulting in a binary NEMO coordinate grid representing the suitable habitat estimate (6). The estimate was rasterized to polar-stereographic projection and small disconnected habitat fragments removed leading to a suitable habitat model (SHM,7). Finally, the model was compared to observational data by binning both to a similar grid (8).
using two-dimensional kernel density estimation (KDE; Habitat spaces). Salinity limits (29-37 g/kg) were added to the TD-spaces, due to missing observational data on salinity, leading to cylindrical three-dimensional TDS-spaces (temperature, depth, salinity). This three-dimensional shape formed the habitat space of Greenland halibut. Decadal averages of bottom temperature, depth, and salinity estimates were acquired from the Nucleus for European Modeling of the Ocean (NEMO) physical oceanography model for the Arctic, north of 40-50 latitude (The oceanographic model, Figure 1). The habitat spaces were then used to limit the NEMO bottom temperature, depth, and salinity grid to estimate suitable habitat for Greenland halibut (The suitable habitat model). Finally, the acquired SHMs were compared to distribution datasets assembled from multiple sources (Distribution data and Model fit).
Our approach has the advantage that the geographic location of samples and the projected suitable habitat become disconnected from each other allowing examination of model fit using the same data as used to estimate the habitat spaces. The concept, however, introduced further assumptions in addition to the above-mentioned limiting factors: (i) the habitat selection of Greenland halibut remained unchanged throughout the period covered by the temperature-depth data , and (ii) our temperature-depth data covered the extremes of these variables representatively. The effects of possible deviations from these assumptions are explained in Discussion and Supplementary Text S1.
It was also assumed that environmental preferences vary with Greenland halibut body size, which serves as a proxy for ontogenetic changes. The fish were allocated to three size classes: "small" fish with total length (TL) [9, 30] cm (Albert et al., 2001a), "medium" with TL [30, 60] cm, and "large" with TL > 60 cm. Greenland halibut smaller than 9 cm are typically pelagic larvae and were excluded from the analysis (Sohn et al., 2010). Greenland halibut have marked sexual dimorphism when it comes to lifespan and size-at-maturity. Male and female length at 50% maturity is around TL 40-45 and 55-70 cm, respectively (Albert, 2003;Gundersen et al., 2010;Hallfredson et al., 2011;Núñez et al., 2015). Consequently, the "small" category consisted of juveniles, the "medium" category contained a mixture of immature individuals and mature males, and the "large" category comprised predominantly maturing or mature females. The Norwegian/Russian trawl survey data were compiled from the IMRs database and covered the Norwegian Continental Slope and the Barents Sea from 2004 to 2019 (Table 1). Bottom temperature from two sources was included: (i) Mean bottom temperatures from Alfredo bottom trawl surveys on the eastern continental slope of the Norwegian Sea 62-80 N. These data were collected using trawl mounted temperature and pressure loggers  (Eriksen et al., 2018). Data collected from Icelandic waters are from a variety of sources including samples from commercial operations, annual scientific surveys, and occasional research cruises. The bulk of the measurements were conducted in the Icelandic spring and autumn trawl groundfish surveys, which started in 1985 and 1996, respectively, and cover the Icelandic continental shelf and surrounding waters (Marine and Freshwater Research Institute, 2010). The scientific cruises, dating back to 1969, had multiple fisheries-related purposes (Magnú sson et al., 1998). The data were compiled from the MFRIs databases that store scientific data on fisheries in Icelandic waters. Bottom temperature was measured with a temperature and depth logger attached to the trawls.

Data for habitat spaces
Greenlandic data came from trawl, longline, and gillnet surveys. Trawl surveys covered the West Greenland shelf and continental slope from 59.5 to 72.5 N and East Greenland slope, from  Habitat model for Greenland halibut 1991 to 2019. Gillnet and long-line surveys covered the inshore fiords of north-western Greenland (Disko Bay, Uummannaq, and Upernavik) from 1993 to 2018. Bottom temperature was estimated with a temperature logger (Starmon Mini) attached to the trawl doors (trawl surveys) or the gear (long-line and gillnet surveys). The Alaskan data were collected during the annual eastern Bering Sea shelf bottom trawl survey (1982,1984, and the biennial Bering Sea slope bottom trawl survey , 2004, 2008, 2012, 2016Stauffer, 2004). Both surveys were conducted from June-August and used the same trawl gear (83-122 eastern trawl). Depth and bottom temperature were measured using a self-contained depth and temperature logger mounted near the centre of the head rope.
Survey data from Canadian waters were available from 1999 to 2017 and covered shelf and slope areas within Baffin Bay, Davis Strait, Hudson Strait, and Ungava Bay. Bottom temperature was recorded using a Seamon temperature logger attached to the trawl door from 1999 to 2004. Since 2004, the bottom temperature has been collected using a CTD attached to the trawl.

Habitat spaces
The SHMs were limited using combined bottom temperature and depth occurrence data for each size group of Greenland halibut (Data for habitat spaces). Bottom temperature (x) and 10-based logarithm transformed depth (y) values for each fish were first binned to 30 by 30 grid, and the number of fish in each grid-cell was calculated. Next, a two-dimensional kernel density estimate (KDE; Chacó n and Duong, 2018) was calculated using the ks package (Duong, 2019) for R (R Core Team, 2020), leading to an estimate of selected conditions in the TD-space. Multivariate kernel density (mKDE) estimation is a non-parametric way of estimating the probability density of a set of variables. The method can be used for many purposes in ecology and statistics due to its flexibility but has a disadvantage of bandwidth parameters influencing the result, making it difficult to exactly reproduce KDE runs. We used the default bandwidth parameters in the ks::kde v 1.11.7 function letting the algorithm optimize the parameters as described in Duong and Hazelton (2003) (refer to Supplementary Text S1 for details). Since the dataset was collected from multiple sources over many years, it likely contained errors. To address this possibility, the 0.999 probability contour was extracted from the kernel density estimates assuming 0.1% of grid-cell probabilities lying outside the main density area were outliers. The threshold was chosen iteratively by adjusting the TD-space borders to encompass most points that were close to each other. The concave shapes of the resulting TD-spaces likely reflected fishing practices, variable correlations, and geographic habitat availability rather than the actual habitat preferred by Greenland halibut (refer to Supplementary Text S1 and Discussion). Therefore, a convex hull was applied to remove any concave shapes. Greenland halibut demonstrated different bottom temperature-depth selection between the Bering Sea (the NOAA dataset) and the North Atlantic (the IMR, MFRI, GINR, and DFO datasets). Consequently, the TD-spaces were calculated separately for the Atlantic and Pacific sides of the Arctic.
Robustness of TD-spaces to single values was examined using bootstrapping: a TD-space was defined for 10 000 temperaturedepth values randomly sampled (with replacement) from the cloud of observations (Steps 1 and 2 in Figure 2; n small ¼ 171 863, n medium ¼ 354 240, n large ¼ 63 053, and n general ¼ 589 156). The procedure was repeated 100 times and the percentage overlap of TD-spaces used as a measure of model sensitivity.
Because much of our survey dataset lacked salinity measurements, we could not use the same approach as detailed above for this parameter. Consequently, the TD-spaces were constrained by the minimum salinity in the available data and the maximum bottom salinity estimated by the NEMO model (S A 2 ½29; 37 g/kg) leading to orthogonal salinity limits compared to temperature and depth in the temperature-depthsalinity (TDS) space ( Figure 2).

The oceanographic model
We used the NEMO-NAA10 km configuration newly developed by the IMR to estimate bottom temperatures and salinities. NEMO-NAA10km is a regional ocean modelling configuration based on the NEMO ocean engine mainly designed to study the Nordic Seas and the Arctic Ocean dynamics (Madec and NEMO System Team, 2015). It covers the North Atlantic Ocean, the Arctic Ocean, and part of the North Pacific Ocean, and utilizes a curvilinear rotated 781 by 888 grid to avoid the North Pole singularity (Figure 1). The horizontal resolution of the model is approximately 10 km, depending on the location. The model uses the ETOPO2 (National Geophysical Data Center, 2001) bathymetry database and has a vertical resolution of 10 to 200 m at the bottom depending on depth. Refer to Supplementary Text S2 for further details about the NEMO-NAA10km model, which will be referred to as "the NEMO model" from here on.
Daily simulations of the model were run from 1 January 1960 to 31 December 2015. Monthly means were calculated for the spatial model grids. The simulations between 1 January 2000 and 31 December 2009 were aggregated over each grid cell using averages to acquire bottom temperature and salinity estimates for the general and size-based Greenland halibut SHMs. Averages of monthly simulations were calculated for each decade from the 1960s until 2010s to study potential changes in Greenland halibut habitat.

The suitable habitat model
The SHMs were generated from mean bottom temperature, depth, and salinity estimates for each cell in the NEMO model. The TDS-space for each size group was used to limit the bottom conditions giving a binary value ("suitable" or "not suitable") for each NEMO grid cell. The curvilinear four-dimensional (two longitude and two latitude dimensions) coordinates for model cells were converted to 350 Â 350 spatial raster grid using the Arctic Polar Stereographic projection (EPSG:3995 with latitude true scale at 71 N, Figure 2) covering the extent of the NEMO model. Disconnected habitat patches smaller than 30 000 km 2 were removed from the SHMs during the conversion. Bering Sea TDspaces were applied to the Pacific side of the Bering Strait and Atlantic TD-spaces were used for the rest of the Arctic with latitude cut-point at 2:5 Â 10 6 m.
The SHMs did not cover the entire range of Greenland halibut in the Pacific, specifically along the west coast of North America, Japan, and Eastern Russia because the NEMO model did not extend over these regions (reviewed in Mecklenburg et al., 2018). Further, the SHMs were unable to estimate Greenland halibut habitat close to land due to the resolution of the underlying oceanographic model. Refer to Supplementary Text S1 for a more detailed explanation of the modelling approach.

Distribution data
Available Greenland halibut data, including those used to develop the habitat spaces, were used to evaluate the fit of SHMs to observations. A size-based dataset was compiled from the sources detailed in the "Data for habitat spaces" section. However, all available data, including fish without bottom temperature or depth records, were used in the model evaluation (Table 1, Figure 3). This size-based dataset consisted of individual fish as data points. Since the size-based dataset was geographically limited, we also compiled an additional general dataset without size, bottom temperature, or depth information (Table 1, Figure 3). The dataset consisted of point observations of Greenland halibut regardless of abundance, size, sampling gear, or sampling time. Both distribution datasets contained presence-only data. The general station-based dataset was compiled from the same sources as the size-based dataset plus: (i) a subset of Norwegian commercial fisheries electronic logbooks for >15 m fishing vessels (2011)(2012)(2013)(2014)(2015) and >13 (Christiansen, 2012) and finally, (vi) reported Greenland halibut occurrences from the East Siberian Sea (Chernova, 2017). The combined dataset spanned across the Arctic but was missing data from sea-ice covered regions around the Arctic Ocean. Since the distribution dataset contained hundreds of thousands of observations over a century of research (Table 1), positions of single fish may be erroneous and should not be over-interpreted.

Model fit
The fit of SHMs to observations was compared by gridding the model area and observational data to 100 Â 100 evenly sized hexagons ($6776 km 2 each). Hexagons even partly covered by a given model with at least one observation were assigned to "habitat with observations", marked using green colour in figures and called "In" in the tables. Hexagons containing the modelled suitable habitat, but no observation were named as "habitat without observations" or "No-obs" in the tables and marked using blue colour. Hexagons that did not contain suitable habitat but contained at least one observation were named as "observations outside habitat" or "Out" in the tables, marked using red colour in figures. Hexagons that did not contain modelled habitat or observations were ignored in the analysis.

Practical implementation
The NEMO model was run using the Fram supercomputer (https://www.sigma2.no/Fram). The TD-spaces and suitable habitats were calculated using R (R Core Team, 2020). The scripts used in the suitable habitat modelling were compiled to an R package and are openly available as supplementary material and online (Vihtakari, 2021a). Maps were made using the ggOceanMaps package (Vihtakari, 2021b) with land polygons from Natural Earth Data and bathymetry from Amante and Eakins (2009).

Temperature-depth occurrence
Mean bottom temperature and depth selection generally increased throughout the ontogeny of Greenland halibut. Small individuals were caught from colder and shallower bottom conditions than the medium and large fish (Figure 3), this was particularly true in the Pacific. Nevertheless, the TD-selection was relatively similar for the size classes within each region, as indicated by the TD-spaces (Figure 4). Fish on the Atlantic side of the Arctic occurred deeper (150-1500 m) than on the Pacific side (50-1000 m; dashed lines in Figure 3).
The separation of realized and available TD-spaces for the North Atlantic at temperatures between À2 C and 0 C was caused by fish caught north and east of Iceland ( Figure 3). Bottom temperatures measured during fishing events were often <0 while modelled NEMO temperatures were >0 for this region. This separation did not lead to a mismatch in modelled suitable habitat since >0 temperatures were already identified as suitable habitat. Further, there was a mismatch at depths >750 m and temperatures between 1 and 4 C for both the Western Atlantic Arctic and North Atlantic. The realized and available TD-spaces overlapped relatively well for the East Atlantic Arctic and the Bering Sea.

The suitable habitat model
The SHM predicted a circumpolar distribution for Greenland halibut with potential habitat in parts of the shallow Bering Strait and the deep Lomonosov Ridge close to the North Pole ( Figure 5). The model estimated shelf seas and continental slopes around the Arctic Ocean as well as the high Arctic Canadian archipelago as a potential habitat for Greenland halibut. Hudson Bay was estimated as a suitable habitat for small and mediumsized Greenland halibut, whereas only a minor area in the central part of the bay qualified as suitable habitat for the large fish.

The model fit
The SHMs were robust to single temperature-depth values throughout the geographic range of SHMs apart from a small region along the coast of Southern Norway ( Figure 5). Size-based SHMs generally explained the available Greenland halibut occurrences with length measurements (Table 2, Figure 6). Suitable habitat was underestimated for medium and large fish west of Iceland and at Rockall bank northwest of Scotland ( Figure 6). The underestimation was caused by the top right corner of TDspace (high bottom temperatures and depths) possibly stemming from the offset between temperature measurements during surveys and the 2000-2009 means of daily NEMO model simulations (Figures 4, 7, and 8). Further, the size-based observational data indicated occurrence of medium and large Greenland halibut in the North Sea along with the Norwegian Trench (n ¼ 1212 and 81, respectively) which the SHMs failed to simulate ( Figure 6). Despite these occurrences, the numbers of Greenland halibut encountered in the North Sea were negligible compared to the Habitat model for Greenland halibut number of observations in the size-based dataset (Table 2). Modelled bottom temperatures in the Norwegian Trench were generally higher than the TD-space limits for Greenland halibut causing the difference between modelled habitat and the rare observations ( Figure 8).
Large parts of eastern Barents Sea contained hexagons without observations of large fish (Figure 6), resulting in 35% of hexagons with no observations. This indicated that the estimated habitat was too wide for the Barents Sea since the observational dataset was representative for this region. A similar effect was evident for medium and large fish in the Bering Sea.
The non-size based general occurrence data demonstrated patterns similar to the size-based data (Figure 6), and also highlighted the model underestimation west of Iceland, but also around Nova Scotia in Canada. Regions such as Hudson Bay, High Arctic Canada, Arctic Ocean basins and ridges, Beaufort, Chukchi, East Siberian, and Kara Seas as well as the Pacific Coast of Alaska contained more modelled Greenland halibut cells without observations than cells with observations.

Limiting factors
Depth limited the modelled suitable habitat in large parts of the Arctic (Figure 8, Table 3). High temperature was a limiting factor along the coast of Europe, the west coast of North America, Rockall bank, and around the Faroe Islands. Further, temperature limited the modelled habitat of large fish in the northern Bering Sea (too low) and around Nova Scotia (too high). Salinity only limited the modelled habitats together with other factors. Both temperature (too high) and depth (too shallow) explained the lack of suitable habitat around the British Isles and in the North Sea. The corners of TD-spaces (i.e. simultaneously limiting temperature and depth) removed suitable habitat from southwest of Iceland, Rockall bank, around Nova Scotia, along with the Norwegian coast, and parts of Hudson Bay, Barents, Kara as well as Bering Seas.

Projected changes in the suitable habitat
The model indicated generally no substantial changes in the suitable habitat of large Greenland halibut from the 1960s to the 2010s (Figure 9). Nevertheless, there was a simulated habitat loss for Hudson Bay, west of Iceland, the Norwegian Trench, and the Rockall bank. Habitat gains were estimated off the Norwegian coast and northeast of Iceland.

Discussion
The SHMs were based on the assumption that only abiotic conditions such as bottom depth, temperature, and salinity would restrict the distribution of Greenland halibut. This assumption is, obviously, a simplification but may nevertheless provide a framework to compare the suitable habitats to observed distributions. The SHMs help to identify aspects we still do not understand in the ecology of the species: they can be used to examine whether other factors influence the distribution of the species and to identify regions that might not have been sampled adequately. This, in turn, allows us to conjecture about the pan-Arctic connectivity of Greenland halibut populations and provide testable hypotheses for future research.
Physiological tolerance curves for a single variable, such as temperature, and a given (ectotherm) species have been established through decades of ecological and physiological research (see e.g. Fry and Hart, 1948;Brett and Groves, 1979;Pörtner, 2002). Here, we applied a similar, but binary, concept using multiple variables. Instead of a measurable physiological response, we used the density of observations as a response variable assuming that the individuals would select their habitats according to the limiting factors. While this assumption has to be based on  Figure 4. Bottom temperature-depth (TD) spaces used in suitable habitat models (blue area), by Greenland halibut size and area. The coloured polygons indicate the regional TD-spaces shown in Figure 3. Shading indicates the model robustness to outliers. The y-axis is on a logarithmic scale. ecological knowledge of the species, the approach has a benefit that the limiting factors do not necessarily need to be based on the physiological limits of Greenland halibut but are rather a consequence of physiological and ecological interactions experienced by individuals. For example, bottom depth where Greenland halibut occurs (50-2000 m) is unlikely to set physiological constraints for a deep-sea adapted species without a swim-bladder (Sebert, 2002). Notwithstanding, of the environmental variables examined in our model, bottom depth limited the habitat of Greenland halibut across the largest area, and is an important factor also in the Barents Sea (Husson et al., 2020). The highest prey densities for medium and large individuals of this opportunistic predator are available in the pelagic layer and demersal habitats on and along the slopes of continental shelves (Michalsen et al., 1998;Dwyer et al., 2010). Consequently, the vertical and horizontal distances to available prey may influence the depth selection of the species. Further, competition between species and predation may influence the depth of occurrence as demonstrated by Nogueira et al. (2017): Greenland halibut occurred shallower than previously after the collapse of the Atlantic cod population in Flemish Cap in the 1990s and shifted back to deeper waters once the cod population began to recover, potentially in response to cod predation. Greenland halibut eggs and larvae are pelagic for months and can be carried large distances by prevailing currents along the continental slopes (Sohn et al., 2010;Domínguez-Petit et al., 2013). The currents may generate upwelling leading to high productivity and attracting juvenile and adult fish to increased prey abundance (Cushing, 1971). Further, hydrographic conditions and distance to spawning grounds can influence the migration of adults (Chumakov and Savvatimsky, 1990) as well as the location of nursery areas. These complex interactions are all interwoven and contribute to the temperature-depth (TD)-spaces studied in this paper. The TD-spaces were based on data from scientific surveys. While these surveys attempt to cover the main distribution of Greenland halibut in their respective locations, the surveys Figure 6. The fit of observations to the general (large map) and size-based (small maps) suitable habitat models. Red hexagons indicate observations outside the modelled habitat, green hexagons observations inside the predicted suitable habitat, and blue hexagons modelled habitat without observations. Dashed lines in size-based maps indicate the coverage of observational data for regions in Table 2. The "Region" column refers to the regional comparisons in Figure 6. The first "Observations" column gives the total number of fish (size-based models) or stations (general model) used in the comparison. The following columns show the percentage of these observations that were within and outside the suitable habitat hexagons. The first "Hexagons" column represents the total number of hexagons used in the model fit comparison. The following two columns indicate the percentage of these hexagons inside and outside the suitable habitat. The last column indicates the percentage of total hexagons that did not contain overlapping observations of Greenland halibut.
certainly have not covered all habitats where the species occurs. For example, maximum bottom depth limits are set somewhat arbitrarily, and the species is likely to occur in greater depths than estimated by the approach. However, surveys at continental slopes, like the one off Norway, extend to 1500 m and are designed to cover the depths of high abundance, typically down to around 1000 m. Thus occurrences at depths >1500 m are likely rare and represent low densities. Our analysis suggests that while temperature limits the southern boundary of Greenland halibut distribution (European coast, North Sea, and the West coast of North America in Figure 8), the northern boundary seems to be limited by other factors. Greenland halibut inhabits shallower waters on the Pacific side of the Arctic compared to the Atlantic side. This difference could be explained by average bottom conditions selected by Greenland halibut on the Atlantic side not being available on the Pacific side ( Figure 3). Notably, the selected average conditions on the Pacific side were available for the West and East Atlantic Arctic regions. Consequently, while high temperatures limit the distribution, other unidentified factors appear to cause the difference in habitat selection between the Atlantic and Pacific sides. The SHMs suggest a connected circumpolar distribution for Greenland halibut as also proposed by Mecklenburg et al. (2018).
While the classic understanding describes the distribution of the species as two separated groups of either Pacific or Atlantic origin with no connection (Smidt, 1969;Godø and Haug, 1989;Bowering and Brodie, 1995), recent literature challenges this view. Orlova et al. (2019) suggested that the Greenland halibut found in the Laptev Sea belongs to the Atlantic strain indicating that juveniles from the Barents Sea spawning region drift to the Russian Arctic (Å dlandsvik et al., 2004). Reports of Greenland halibut from the East Siberian, Chukchi, and Beaufort Seas (Rand and Logerwell, 2011;Logerwell et al., 2015;Chernova, 2017;Majewski et al., 2017) further indicate that connected circumpolar Greenland halibut populations might exist; yet likely in low abundances as indicated by low catches of the species from these ice-covered Arctic shelf seas.
Studying the genetic connectivity of Greenland halibut around the Arctic could test the hypothesis of connected pan-Arctic populations. Genetic similarity between specimens from the Pacific side of the Arctic Ocean (mainly Beaufort, Chukchi, East Siberian, and Laptev Seas) and the Atlantic side would give support for the model, while lack of such connectivity could be used to falsify the hypothesis in favour of the classical understanding of Greenland halibut distribution. The connection between the Pacific Bering Sea population and the putative Arctic Ocean population is particularly interesting. While our model suggests a relatively large geographic distance between the eastern Bering Sea spawning area (Sohn et al., 2010;Duffy-Anderson et al., 2013) and the suitable habitat in the Arctic Ocean, it is still possible that individuals move between these locations as has been shown for the Atlantic side where Greenland halibut migrate from Svalbard to Iceland (Albert and Vollen, 2015). So far, the only studies comparing the Bering Sea population with Atlantic populations using microsatellites indicated that these populations might be genetically separated (Orlova et al., 2017(Orlova et al., , 2019. Further, studies have suggested that Greenland halibut is divided into two populations in the North Atlantic: the Northwest containing Newfoundland, Davis Strait and Baffin Bay, and the Northeast containing Iceland, Norway, and Russia (Knutsen et al., 2007;Westgaard et al., 2017), while another study in the Northwest Atlantic found evidence for panmixia (Roy et al., 2014). Low gene-flow between the hypothesized Atlantic populations could indicate that there are several relatively isolated Greenland halibut populations across the Arctic due to the ecology of the species and potential migration barriers. The abovementioned studies, however, have been far from conclusive and a pan-Arctic study on the genetic structure of Greenland halibut is required to better understand the connectivity of populations and migrations of individuals among them.
Localized spawning outside the three major spawning regions has been suggested to contribute recruits, especially to the Greenland halibut fishery around Iceland and East Greenland (Gundersen et al., 2013). If localized spawning was an important part of Greenland halibut life history strategy, the gene flow between populations, especially within the hypothesized lowabundance regions around the Arctic Ocean, would show high inter-annual variability through variable recruitment and possibly lead to low genetic differentiation over time. Further, the connectivity would not be dictated by larval drift from the three major spawning regions, but multiple minor spawning locations could contribute to the genetic structure of populations. Such a scenario would make our SHM a more realistic distribution model of Greenland halibut as the model lacks larval drift simulations  larvae would influence the distribution of Greenland halibut as a whole and therefore the lack of larval drift simulation within the model likely leads to an overestimation of the potential distribution of the species. Incorporating larval drift into the model would be a useful exercise for future studies. Such studies could examine the influence of different spawning location scenarios across the Arctic on the observed distribution. Future field surveys should aim to capture pelagic Greenland halibut larvae outside the known spawning regions. Our model suggests that regions, which have been or are typically covered by perennial ice (Lomonosov Ridge, High Arctic Canada, and continental slopes along the Central Arctic Ocean) are potential habitat for Greenland halibut. While only limited observations of the species from ice-covered regions exist, the demersal fish communities have generally not been adequately studied in locations where bottom trawling is difficult. For example, the model suggests Lomonosov Ridge as a suitable habitat for adult Greenland halibut. So far, <15 fish species have been recorded from the Central Arctic Ocean and Greenland halibut is not among them (Christiansen and Reist, 2013). Another interesting suggestion by the model is that Hudson Bay could be, at least partly, suitable habitat for Greenland halibut. Despite many observations of the species in Hudson Strait and records of other fish species in central Hudson Bay (Fisheries and Oceans Canada, 2019), Greenland halibut has not been reported from Hudson Bay (Mecklenburg et al., 2018). One explanation might be the prevailing current from Hudson Bay through Hudson Strait that runs west to east and is fresher and colder than the Atlantic water current that flows east to west along the north coast of Hudson Strait. Another explanation for these offsets could be that the species is present in the locations, but in very low abundances, since the model does not make any predictions of the abundance of Greenland halibut. Consequently, the presence of the species could have gone unnoticed given the low sampling effort in most of the Arctic. The model may, and very likely does, overestimate the range of Greenland halibut, because the species is likely not only constrained by bottom temperature, depth, and salinity as discussed previously. For example, the above-mentioned lack of larval drift to these locations would restrict the distribution of the species. The underlying oceanographic model and the decadal medians for NEMO model grid cells are also likely to bias model predictions. Further, the parameters selected during the definition of the environmental TD-spaces from kernel densities could lead to slightly differing SHMs.
The SHM indicates that Greenland halibut may have lost parts of its habitat from the 1960s to the 2010s in the region west of Iceland, which currently supplies 5-25% of the Icelandic survey biomass estimate (MFRI, 2019). This projected habitat reduction is due to ocean warming. Therefore, Greenland halibut may shift their range northward within this region, and future catches by Iceland may concentrate on the region between Iceland and Greenland. A change in fishing areas would affect the economy in these fisheries. On the other hand, the opening of the Arctic Ocean (Screen and Deser, 2019) and estimated suitable habitat along the continental shelf margins of the Arctic Ocean may lead to new fishing opportunities and increased catches of this species globally (Christiansen, 2017, and references therein). Nations such as Canada, the United States, Russia, Norway, and Greenland/Denmark may benefit from the projected changes to Greenland halibut distribution.
Species distribution models have been subject to intensive research for the past decades, and a variety of approaches have been described in the literature (see e.g. Elith and Leathwick, 2009). Due to the abundance of literature, wildly varying terminology, and simplicity of our approach, we cannot exclude the possibility that the method has already been described, although we are not aware of exactly similar approaches. Unlike many common niche-or environmental space-based approaches (Elith and Leathwick, 2009), our model design is close to empirical. The values of environmental variables are acquired directly from observations, without the spatial model stage (Elith and Leathwick, 2009), making the spatial observations and geographic projections independent of each other. The kernel-density estimation is  designed to avoid assumptions regarding the distribution of points in the habitat space, setting it apart from many alternatives using parametric statistics (Elith and Leathwick, 2009 and references therein) or complex combinations of such statistics (Elith et al., 2011) to estimate the probability density. The use of convex hulls to correct for correlations and biases in the underlying data sets our method apart from the closest analogy (Broennimann et al., 2012) where the authors used multivariate methods to enable comparison of niches. Our approach is simple, flexible, non-linear, and can be run using multiple limiting variables (up to 6 at the time of writing). It works well as an exploratory model to guide the examination of unknown species distributions, their relation to limiting factors, and to generate testable hypotheses.
The flexibility of the model comes at the cost of utility, however. Instead of making assumptions on the distribution of data, the approach requires a large observational dataset and an assumption that the used environmental variables limit the distribution. Further, it is assumed that the extremes of observational data represent all possible combinations of environmental variables where the species is encountered across the domain of the projection. These assumptions make sense when applied to the convex hull boundaries of habitat spaces but may not apply to the probabilities inside the habitat spaces as presence-only observations tend to be heavily biased by sampling. The probabilities cannot be estimated reliably without assumptions on the distribution of data and their multidimensional relation to each other. Consequently, the model only provides binary suitable habitat estimates without any predictions of the probability of occurrence. This is a limitation that could be overcome by using a smaller domain of projection and better data (e.g. frequencies of occurrence). There are, however, other methods that may be better suited for such cases (Elith and Leathwick, 2009).

Conclusions
Despite shortcomings, the model does indicate that bottom depth and temperature are major constraining parameters of Greenland halibut distribution and that even such a simple modelling approach does succeed in simulating the known distribution of the species. In this article, we used a suitable habitat model to pinpoint poorly surveyed regions where Greenland halibut may occur, and to raise questions for future research (e.g. genetic studies and larval drift modelling). We did not attempt to estimate the abundance or probability of occurrence of Greenland halibut across its range as such a model would have required further assumptions, absence data, and knowledge about the ecology of the species. Our model indicates the presence of a suitable habitat for a potential connected circumpolar distribution of Greenland halibut. The on-going melt of sea-ice and associated changes in marine productivity may result in a northward shift in Greenland halibut distribution that could negatively impact some fisheries, but allow for the opening of others in previously unfished areas of the Arctic. Demersal fish surveys in poorly explored regions and genetic analyses to study the connectivity of putative Greenland halibut populations are required to test our predictions.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.