Hierarchical Bayesian model reveals the distributional shifts of Arctic marine mammals

Our aim involved developing a method to analyse spatiotemporal distributions of Arctic marine mammals (AMMs) using heterogeneous open source data, such as scientific papers and open repositories. Another aim was to quantitatively estimate the effects of environmental covariates on AMMs’ distributions and to analyse whether their distributions have shifted along with environmental changes.


| INTRODUC TI ON
The decline of the sea ice has changed the Arctic landscape and the habitats of the marine species in the Arctic marginal seas (Durner et al., 2009;Gaston, Gilchrist, & Hipfner, 2005;Laidre et al., 2015).
Longer ice-free periods and less ice have facilitated marine traffic and extraction of natural resources (Smith & Stephenson, 2013), creating environmental risks such as oil spills (Nevalainen, Helle, & Vanhatalo, 2017). Accurate information on species' distributions helps to assess species' vulnerability to changes in their habitats (Laidre et al., 2008) and to prevent their exposure to human caused hazards (Helle, Jolma, & Venesjärvi, 2016). Arctic marine mammals (AMMs) spend most or all of the year in sea areas (Laidre et al., 2008) and their specific activities (foraging, moving, breeding and denning) depend on ice cover and type (Bluhm & Gradinger, 2008;Kovacs, Lydersen, Overland, & Moore, 2011). By inferring how habitat characteristics (sea ice, depth, distance to the coast and hydrography) correlate with the abundances of AMMs, we can predict the distributions of AMMs. The predictions can be utilized in planning conservation actions and in assessing the risks of different species-human interactions Wilson, Trukhanov et al., 2017), whereas habitat utilization functions can underpin spatially explicit demographic analysis for better population assessments (Kearney & Porter, 2009;Lunn et al., 2016). In areas with high survey costs and little designed survey data, flawed population assessments may lead to improper conservation or harvest actions (Regehr, Wilson, Rode, Runge, & Stern, 2017).

Species distribution modelling (SDM) is a cost-efficient method
for studying how species respond to environmental covariates and where they occur (Elith & Leathwick, 2009). The scarcity of in situ observations on AMMs has so far restricted the attempts to assess their densities in large areas (as in Matishov, Chelintsev, Goryaev, Makarevich, & Ishkulov, 2014). Utilizing species observations from complementary sources poses a high potential to overcome the data shortage (Pacifici et al., 2017) but increases the uncertainty about the spatiotemporal accuracy of observations and about the survey effort, and thus evokes the need for methodological development (Guillera-Arroita, 2017;Warton, Renner, & Ramp, 2013).
We developed distribution models for three AMMs, polar bears (Ursus maritimus), Atlantic walruses (Odobenus rosmarus rosmarus) and ringed seals (Phoca hispida), and studied their seasonally varying areal densities in the Kara Sea ( Figure 1). We analysed previously F I G U R E 1 The study transect lines (red) and species observations (green) for polar bears (a-c), walruses (d-f) and seals (g). A grid cell through which a transect goes or at which there is a species observation is treated as a cell with observation. The observation is the number of reported individual species members and zero otherwise. Other grid cells are treated as missing data and do not contribute to the analysis. The s 1 and s 2 axis (h) denote the axes of the coordinate system used in spatiotemporal model [Colour figure can be viewed at wileyonlinelibrary.com]

(a) (b) (c) (d) (e) (f) (g) (h) (i)
published data on species observations that had been collected by heterogeneous sampling methods and were incompletely reported, for example in some cases containing presence-only observations.
Along with the hierarchical modelling, the Poisson point process (PPP) has become an important modelling methodology for data with varying and potentially uncontrolled sampling effort (Dorazio, 2014;Fithian et al., 2015;Warton & Shepherd, 2010). In a spatially discretized (grid based) model setup PPP is an efficient tool for modelling species' density per area, which is a more intuitive and accurate measure of density than the scale dependent density (Maxent;  or occurrence probability (logistic regression ;Fithian & Hastie, 2013) per grid cell. The novelty of our work is in tailoring a recently developed SDM methodology to answer a question that has been out of reach with past modelling tools.
We quantified the utilized habitat characteristics of the species, created a hindcast of density predictions and based on predictions, evaluated potential changes in species' relative densities in the Kara Sea during the study period. This is essential information in the Siberian shelf area, which has been little studied compared to other marginal regions of the Arctic (Wassmann, Duarte, Agustí, & Sejr, 2011). The examined species are high in the marine food web and represent the full diversity of marine mammals in the Kara Sea excluding beluga whales (Delphinapterus leucas). The predator species are essential for the functioning of the marine ecosystem as they control populations of species below them in the food web and the balance in utilization of resources (Baum & Worm, 2009;Myers, Baum, Shepherd, Powers, & Peterson, 2007). Despite earlier estimates of the population sizes of polar bears (3,200 ± 1,100 individuals [Matishov et al., 2014]), ringed seals (90,000-150,000 [Kelly et al., 2010]) and walruses (<500 [Born, Gjertz, & Reeves, 1995]  , the lack of reasonable covariates or a mismatch of spatiotemporal scales between covariates and observation records still constrain the utilization of the data (Rocchini et al., 2011). There are also other poorly documented Arctic areas alongside the Kara Sea, where biodiversity estimates would benefit from a more efficient use of available data (Laidre et al., 2015).

| Data
Our data consist of species observations from the years 1996 to 2013 and of environmental data (the covariates). We use bathymetry, distance to the coast, ice concentration and sea surface salinity (SSS) as environmental covariates to explain species observations. The covariates were selected based on their accessibility and hypothesized impact on species occurrence (see Tables 1 and 2). The study area was defined as a region between the Eurasian continent, Novaya Zemlya and Severnaya Zemlya (see Figure 1). We used a 5 × 5 km lattice grid over the Kara Sea as our modelling layer (see Section 2.2) and a temporal resolution of 1 month. The spatial and temporal resolutions compromise between computational effort and a temporal range of ecosystem responses to environmental variation (Mannocci et al., 2017). More information about the spatiotemporal scaling and spatiotemporal variation of covariates is provided in the supplementary material (Appendix S1).
We collected the species data from scientific articles and books (see Table 3).  There were some exceptions in the data quality between source studies (see Table 3). The observations of Matishov et al. (2014) Pilfold et al. (2014), Reeves (1998), Smith, Hammill, andTaugbol (1991) and Tynan and DeMaster (1997) Sea surface salinity (SSS) SSS is a surrogate for the productivity of the marine ecosystem in the Kara Sea. SSS is a regionally important covariate for the ecosystem functioning and its' spatial trend may be one factor determining the distribution of Arctic marine mammals Ref. Bluhm and Gradinger (2008) and Miquel (2001) Distance to the coast Polar bears inhabit coastal or pelagic areas according to the subpopulation. In many parts of the Arctic and in the Kara Sea polar bears inhabit coastal land fast ice and avoid open sea In summer walruses stay on the coast and dive in the coastal zone for benthic prey. Mostly they do not follow the marginal ice zone. This means that they do not leave too far from the coast Ringed seals are mostly coastal, but their density varies according to the distance to the coast. They are not dependent on depth but more on the ice type varying along with the distance to the coast Ref. Born and Knutsen (1997), Ferguson et al. (2000), Mauritzen et al. (2002Mauritzen et al. ( , 2003 and Pilfold et al. (2014) Freitas et al. (2009) Freitas, Kovacs, Ims, and Lydersen (2008), Krafft et al. (2007) and Reeves (1998) Relative seal density Ringed seal is the most important prey for polar bears and foraging habitat of polar bear follows the distribution of seal denning --Ref. Born and Knutsen (1997), Derocher, Wiig, and Andersen, (2002) and Stirling and Oritsland (1995) --Bathymetry Not included. Bathymetry is an important covariate in Arctic wide studies reflecting shelf edges and basins but our study area is located mostly in the shelf area where the effect of bathymetry on the density of polar bears is less important than in other shelf seas Walruses feed on benthic species, which makes them dependent on the bathymetry. Walruses are able to dive deeper than 250 m, which is enough for them to access the bottom in the most of the shelf area in the Kara Sea observations. We linked this information to our study grid by manually setting presence observations to the grid cells corresponding to the reported locations at the reported time.

| Data analysis
We used the hierarchical Bayesian PPP model to analyse the data.
A rationale for using a Bayesian approach is that it provides tools to combine heterogeneous data through hierarchical model structures that account for variations in data originating from spatiotemporal dynamics in species' density and differences between data collection methods (Gelfand et al., 2005;Latimer, Banerjee, Sang, Mosher, & Silander, 2009). A PPP is a point process that is a widely used building block for many spatial abundance data (Banerjee, Gelfand, & Carlin, 2015;Gelfand, 2010) and in recent years its potential for analysing opportunistic, non-design based data has been demonstrated by several authors (Dorazio, 2014;Warton et al., 2013;Yuan et al., 2017). Here we give an overview about the hierarchical model structure and how covariates and random effects were included in the model. More information on how our methodology arises from the PPP framework is given in the supplementary material (Appendix S2).
The Bayesian hierarchical framework allows us to model the observations, density processes and process parameters at separate levels (Wikle, 2003). In the observation model level, we describe the conditional distribution of the number of individuals, y(s,t,j), observed in a grid cell with coordinates s (kilometres) and at time t (months) during survey j (in total three surveys of polar bears and walruses each, and one survey of seals are summarized in Table 3) with a negative Binomial distribution function, where the latent function, f(s,t,x s,t ), denotes the logarithm of the relative density of a species, x s,t the vector of environmental covariates at grid cell s at time t, ε j the effect of sampling effort of survey j and r the overdispersion parameter. The negative Binomial distribution is an overdispersed version of the Poisson distribution.
We parameterized it as in Vanhatalo et al. (2013)  corresponds to decreasing variance, and at the limit, as r approaches infinity, the negative binomial approaches a Poisson distribution.
The overdispersion parameter r accounts for spatially and temporally uncorrelated variation that is not explained by covariates or spatiotemporal random effect or sampling effect (to be described below).
The original publications do not contain information on the probability of observing an individual nor on sampling effort. The only information about these is by Matishov and Dzhenyuk (2007) (2007) 517/8963 G TA B L E 3 Summary of species data and references to data sources. Matishov, Goryaev, and Ishkulov (2013) and Svetochev and Svetocheva (2008) have only positive abundance information and lack the information of total survey area. Months of seasons are as following: winter (12, 1, 2), spring (3, 4, 5, 6), summer (7,8,9), autumn (10, 11) covariates on both density metrics follow the same functional form.
Also, the differences in expected relative and absolute densities between locations or time points are proportional to each other. The analysis of expected relative densities allows us to assess the effects of environmental covariates on species' densities and the spatiotemporal trends of species' densities, but it does not allow us to assess the expected population sizes.
The spatiotemporally constant parameter for sampling effect ε j adjusts for the variability in the abundance counts originating from the varying sampling methodologies between different surveys (that is, different data sources). The sampling effect was modelled as independently and identically distributed Gaussian random variables, where σ 2 s is the variance that governs the variation in the sampling effects. The effect is included in the distribution models of polar bears and walruses. Seal observations originate from a single source and, hence do not vary depending on the survey. The unstructured random variation in sampling effect in grid cells within a single survey is modelled with the overdispersion of Negative Binomial distribution.
The log relative density was modelled with an additive function where α is a constant intercept for the areal and temporal average, T is an N × 1 vector of coefficients and g(s,t) is a spatiotemporal random effect which captures spatiotemporal variation that cannot be explained by the covariates (Gelfand et al., 2005;Vanhatalo, Hosack, & Sweatman, 2017). We standardized all covariates to have zero mean and standard deviation of one in order to help the assessment of their relative importance for explaining the data. The vector of covariates, x s,t , included all the covariates and their squares so that the responses along covariates were assumed to be quadratic. This is justified as the studied species may have favourable conditions in the middle of the environmental gradients and thus their responses would follow a hump-shaped form (Elith & Leathwick, 2009). A spatiotemporally varying random effect is given a Gaussian Process (GP) prior. GPs are a family of stochastic processes, which define probability distribution over functions. They are a flexible tool for modelling dependency between observations in space, time and covariate space (Golding, Purse, & Warton, 2016;Rasmussen & Williams, 2006;Vanhatalo, Veneranta, & Hudd, 2012;Vanhatalo et al., 2013). A GP is defined by its mean and covariance function. Here we used mean zero and a separable covariance function that is a product of squared exponential spatial and exponential temporal covariance functions where σ 2 ST is the process variance and l i , i = 1, 2 and l s are the lengthscale parameters governing how fast the correlation between g(s,t) and g(s′,t′) decreases (Rasmussen & Williams, 2006).
In addition to abiotic effects, we explained the log relative density of polar bears also with the maximum a posteriori (MAP) estimate of relative density of seals. The relative density of seals was treated in the model as a spatiotemporally varying covariate. To some extent, polar bears follow the distribution of seal lairs (Pilfold, Derocher, Stirling, & Richardson, 2014). Although, polar bears may reduce the seal population (Stirling & Oritsland, 1995), according to previous studies, the spatial correlation between the two species is dictated merely by polar bear presence being dependent on seal presence and not the other way around (Ferguson, Taylor, & Messier, 2000). Hence, we assumed that the species interaction works only in one direction. We modelled the effect of seals' relative density on polar bears' log relative density with a Michaelis-Menten function, f(x s,t ) = ax s,t /(b + x s,t ), which is commonly used in ecology for responses that first increase or decrease and then saturate. It defines an asymptotic response between the log relative density of polar bears and the relative density of seals where a is for saturation level and b for half-saturation point.
The last level of hierarchy is the parameter model which defines the prior distributions for the parameters of the process model (Wikle, 2003). We gave vague priors for the intercept and regression coefficients encoded by mutually independent zero mean Gaussian distributions with large variance; that is, β i ~ N(0,10) for all i and a ~ N(0,10). The variance of study effects and the process variance were given weakly informative half Student-t priors, σ 2 s ,σ 2 ST ∼ Student − t + (0,1). Similarly, the inverse length-scales of the spatiotemporal random effect were given Student-t priors, 1/l i ~ Student − t + (0,0.1) which favours smooth spatiotemporal trends. The overdispersion parameter of the negative Binomial distribution was given a gamma distributed prior with r ~ Gamma(2,.1).
The half-saturation point of the Michaelis-Menten function was given a Gaussian prior b ~ N(0,10).
The models' hyperparameters were estimated with Markov chain Monte Carlo (MCMC) sampling using the GPstuff toolbox (Vanhatalo et al., 2013). The convergence of Markov chains was analysed with the Gelman-Rubin Potential Scale Reduction Factor (PSRF). The models were validated with posterior predictive checks and cross-validation (Gelman et al., 2014). In addition, we compared two polar bear models using leave-one-out cross-validation; the one described above and another where the relative density of seals was removed leaving only environmental covariates.
The models were used to predict the relative density of the species in the Kara Sea in each month in the years 1997 to 2013.
In order to assess the effect of spatiotemporal random effect, we made two separate predictions: one with the full model and another based solely on the covariate effects (for discussion on this kind of separate predictions see e.g., Vanhatalo et al., 2017). If the spatiotemporal random effect has a significant effect there should be difference between these two predictions. We summarized these predictions by calculating the average relative densities in four seasons (December-February, March-June, July-September, October-November) by averaging the expected values of relative densities over the Kara Sea over the months of a specific season. We also made a comparison between average relative densities in spring season between the first (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) and the second half (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) of the study period. This comparison was done using the predictions based solely on the covariate effects in order to estimate the effect of changing environment. The spring season is best represented by observations of all species and moreover, it is a denning and hunting season for seals and polar bears, respectively (Stirling & Derocher, 2012).

| RE SULTS
The posterior predictive checks did not show significant deviations between predicted and observed species' abundancies. All the model parameters identified well with MCMC sampling and the model results were consistent in cross-validation tests where leaving subsets of data out did not alter the results significantly. At least for polar bears and seals, the results are in line with earlier knowledge on the distributions and hence, the data used in this study was adequate for inferring their responses to environmental covariates.
The standard deviation of the spatiotemporal random effects was at the same order of magnitude as the variation of the log relative density along the environmental covariates ( Figure 2; Table 4) which indicates that there were significant deviations in the species' distribution patterns from that predicted only by the environmental covariates.
The spatial length-scale parameters were in the order of tens to hundreds of kilometres indicating smooth spatial random effects across the study region. However, the temporal length-scales were significantly <1 month which indicates that the spatiotemporal variations did not contain temporal trends. The variance of the study effect term in the polar bear and walrus models was of the same order of magnitude as that of the spatiotemporal random effect which indicates significant The posterior of species' responses to environmental covariates with the 95% confidence intervals. The responses are plotted as changes of log relative density over the range of covariate values in the data. The curves are scaled to start from zero. In case of quadratic responses, the location with width zero confidence interval corresponds to the empirical mean of covariate values in the data where the unscaled curve would cross the zero TA B L E 4 The posterior mean and 95% confidence interval (in parenthesis) of hyperparameters of the distribution models. Spatial length scales are in kilometres and temporal length-scales in months. The directions for the length scales are shown in Figure 1 Model variation in detection or reporting probability between studies. Lastly, the overdispersion parameter was small indicating significant overdispersion compared to a Poisson model.
The log relative density of polar bears was explained the most by the relative density of seals. We compared models with and without the seal parameter with a leave-one-out cross-validation using the log predictive density statistics (Vanhatalo et al., 2012) (−0.1053 with the predicted relative density of seals and −0.1137 without one). The cross-validation log predictive density measures how well the model predicts test data and has the greater value the better the prediction is. The response of polar bears to seals was estimated to saturate around the seals' relative density of 4.3 (Table 4). However, this estimate was outside the range of predicted seals' relative densities for which reason the response is almost linear in that range. The log relative density of seals was explained mostly by ice concentration (Figure 2). The response to ice concentration peaked around 70% ice cover with high certainty.
The inference of the log relative density of walruses suffered from the lack of data, and the estimates of responses to covariates came with high uncertainty. Walruses' log relative density was explained mostly by distance to the coast (Figure 2) as their density drops further than 70 km from the coast.
According to the predicted hindcast of each species, the relative densities of seals varied the most between the spring seasons of the first and the second half of the study period (Figure 3). Their relative densities decreased in the Northeastern Kara Sea and increased close to the eastern coast of Novaya Zemlya. The changes in relative densities of seals and polar bears had a similar spatial pattern, but for polar bears, the changes were smaller. The relative densities of walruses decreased slightly across the Kara Sea.

| D ISCUSS I ON
The hierarchical Bayesian model framework provided us with tools to treat the assumed inaccuracies in the heterogeneous data.
F I G U R E 3 The average relative densities of polar bears, walruses and seals averaged over the winter (months 12, 1, 2), spring (3, 4, 5, 6), summer (7,8,9)  Applying a linear model on logit transformed presence and absence observations would have exposed the response estimates on bias originating from spatiotemporal autocorrelation and varying survey effort. Utilizing data from many seasons allowed us to track the seasonally varying species' densities across the Kara Sea, which is of interest for conservation actions.
In previous studies, polar bears' abundance and population trends have mostly been related to ice cover and type (Lunn et al., 2016;Regehr, Lunn, Amstrup, & Stirling, 2007) and attempts to relate polar bears' occurrence to seals have not proved a spatial dependency between them . However, according to the estimated covariate responses and model comparison, the effect of relative density of seals outweighs the effect of ice concentration on the log relative density of polar bears. Even though the effects of the relative density of seals and ice concentration cannot be fully disentangled, as seals are heavily dependent on ice, the results indicate that the relative density of seals has a clear positive effect on the relative density of polar bears.
When excluding the seal covariate, polar bears' response to ice concentration becomes hump-shaped peaking around 70% of ice cover (see Appendix S1), which is similar to the responses found by Durner et al. (2009) and Lone et al. (2018). Thus polar bears and seals follow a similar hump-shaped relationship to ice concentration. When the relative density of seals is included in the model, the response to ice concentration is linearly increasing and assumingly shows the independent effect of ice concentration. This indicates that in a simple trophic system or in case of a highly specialized predator, the occurrence of a prey species is a more informative covariate than an environmental variable for predicting the density of a predator. This has also been recognized in more diverse marine systems (Reisinger et al., 2018). Assessments of polar bears' distribution could be improved by linking areal estimates of seals' density to a RSF (Resource Selection Function) of polar bears (Durner et al., 2009;Lone et al., 2018;Wilson et al., 2014). Reviewers of the manuscript were interested in the reasons to leave bathymetry out from the models of relative densities of polar bears and seals.
We assumed that bathymetry would not have a strong effect on their relative densities, as the effect has been recognized mostly in the shelf breaks and our study area does not cover a shelf break zone. We also carried out a model comparison, which supported choosing the models not having bathymetry as a covariate.
The estimated covariate responses of seals support the earlier hypothesis about seals' habitat characteristics in spring season.
Their utilized ice habitat varies from stable land fast and pack ice to more unstable and productive polynyas and leads depending on their sex and offspring (Krafft, Kovacs, & Lydersen, 2007;Stirling, 1997). This is adequately shown by the positive response to moderate and high values of ice concentration. The estimated effect of distance to the coast speaks for seals inhabiting more pelagic than coastal sites, which is made possible by the wide land fast ice zone in the continental shelf (Pavlov & Pfirman, 1995).
The responses of walruses to covariates come with higher uncertainty than the responses of polar bears or seals. The estimated responses of walruses support the assumption that the Atlantic walruses stay mostly in coastal shelf areas, where they feed on benthic vertebrates (Lydersen, Chernook, Glazov, Trukhanova, & Kovacs, 2012). Hence, ice concentration does not affect much walruses' density pattern in the Kara Sea.
With the covariate responses, we can estimate past changes in the expected relative densities. Arctic wide vulnerability assessments do not consider region specific distributional changes which may actually support species relocation inside the area instead of disappearance (Stirling & Derocher, 2012;Wilson, Regehr et al., 2017). We assume that the shrinking ice cover has caused the decrease of polar bears' relative density in coastal regions (Figure 3).
This highlights the sensitivity of polar bears to changes in ice conditions and supports earlier studies (Durner et al., 2009;Lunn et al., 2016). The slight increase in the relative density of polar bears in the Western and Eastern pelagic Kara Sea may be due to the increased relative density of seals in those regions.

Seals have had opposing trends in the Eastern and Western Kara
Sea due to the lowering ice concentration in both regions. In the Western Kara Sea, the average ice concentration has been lowered close to the optimum of seals' habitat characteristics, whereas in the Eastern Kara Sea ice concentration has dropped below the optimum level. Seals are hypothesized to be less susceptible to suffering from shrinking ice cover as their habitat requirements are more flexible than those of polar bears (Laidre et al., 2008), which is supported by our results. In addition to polar bears, also walruses are hypothesized to be site specific species and thus sensitive to decrease in ice cover. Coastal habitats may maintain small walrus populations, which may be the case in the Kara Sea (Laidre et al., 2008). However, the lowering ice concentration has also decreased the relative density of walruses in the coastal regions. According to our results and the forecasted decline of the average ice concentration (Wang & Overland, 2012), each AMM may have distributional changes ahead as the Southern Kara Sea becomes ice free for a longer season in the future.
The challenge of analysing incomplete and heterogeneous biological data was overcome by thinning point process in relation to unknown survey effort and by explaining relative densities with spatiotemporal random effects. In general, random effects can be used to correct for possible biases in fixed effect estimates in cases where data do not have clearly defined or reported survey effort. Properly defined random effects capture the excess variability in species' relative density that is not explainable by environmental covariates and hence, improve also the estimates for the covariate responses (Ovaskainen et al., 2017). In the Kara Sea, ice type may affect species' densities in such a way that it cannot be explained solely by ice concentration. However, we can expect that such a variable has a spatiotemporally structured effect which can be dealt with by using a spatiotemporal random effect (Ovaskainen, Abrego, Halme, & Dunson, 2015;Vanhatalo et al., 2012). Other possible sources of spatiotemporal variation are prey availability for seals and walruses and species' seasonally varying behaviour Jay, Fischbach, & Kochnev, 2012;Mauritzen et al., 2003). Spatially smooth random effects indicated that there is some environmental variation that has not been included as a covariate in the model. The low temporal variation of the random effect might be due to the strong temporal variation of ice concentration, which keeps the species' densities constantly moving.
The other random effect accounted for the survey bias originating from varying survey protocols (Dorazio, 2014;Fithian et al., 2015). Sampling bias is typically induced by presence-only obser- The procedure of creating species observations from tables and maps created some inaccuracy in data. We estimated the digitizing error by calculating the width of the transect line on the source map. The error is 16 km in Matishov et al. (2014) and Matishov and Dzhenyuk (2007) and 33 km in Glazov et al. (2013), which make a width of three and seven study cells, respectively. The temporal information of cruises was presented as the start and end dates and many cruises covered periods from 1 to 3 months. We consistently chose the central point of the time frame to represent the cruise transect, which may create some temporal error. The spatial and temporal uncertainties related to transects is not expected to produce large errors; the spatial error was in the same order of magnitude as the resolution of the original ice concentration data and covariates did not vary significantly within the spatial error or along survey transects during the survey periods. The studies without survey transects (Table 3) were included in our analysis as the recorded observations contained numbers of observed animals which is informative for the response of the log relative density to environmental covariates even when we do not have information about survey transects (Dorazio, 2014). The benefit of combining systematically and opportunistically surveyed data is that detection probability can be estimated with the former data and hence the latter data can be used along with other data for estimating a model's fixed effects (Dorazio, 2014;Giraud et al., 2016;Pacifici et al., 2017

| CON CLUS ION
We demonstrated how several, heterogeneous, open source data sets can be jointly analysed within the PPP framework to produce new information on AMMs' distributions. Our results suggest that the relative densities of polar bears and walruses have decreased or stayed close to constant in the Kara Sea during the last 20 years and that the distribution of seals has shifted from the Eastern to the Western Kara Sea. The decrease in the average ice concentration across the study region has driven these changes. The spatial dependence of polar bears on seals was significant. This demonstrates that in a simple trophic system, modelling the density of a top predator benefits from taking into account the density of prey species compared to using environmental variables.
Combining open data from different sources created a fairly large but heterogeneous data set for analysing AMMs' distributions. Due to heterogeneity in the data sources and uncertainty concerning sampling techniques and effort, the complex spatiotemporal variation of the data needed to be modelled with care. After accounting for those uncertainties, we were able to produce useful new knowledge on AMMs' distributions during a 17-year-long period. The approach is cost-efficient as it allows the analysis of the vast amounts of existing environmental data. Hence, our approach provides important advances for conservation efforts in these areas by providing a method to build improved information on distributional changes from opportunistic studies.