Simulating changes in polar bear subpopulation growth rate due to legacy persistent organic pollutants – Temporal and spatial trends

a b s t r a c t a r t i c l e i n f o


Introduction
Arctic ecosystems are subjected to many threats induced by human activity. Especially polar bears (Ursus maritimus) have received much attention, as these species are suspected to be significantly impacted by climate change, with sea ice decline hindering these species in finding prey (i.e. seals) on the ice (Jenssen et al., 2015). Additionally, hunting and exposure to toxic pollutants (e.g. persistent organic pollutants; POPs and mercury (Hg)) are also considered to be harmful to polar bear populations Jenssen et al., 2015;Letcher et al., 2010;Nuijten et al., 2016). Polar bears depend on a lipid-rich diet, mainly consisting of ringed seal (Phoca hispida) and bearded seal (Erignathus barbatus), sometimes also including other prey, such as hooded seals (Crystophora cristata), beluga whales (Delphinapterus leucas), narwhals (Monodon monoceros) and walrus (Odobenus rosmarus) Stirling and Archibald, 1977). Due to their diet and the bioaccumulative nature of POPs and mercury, high levels of these compounds have been detected in polar bears in several subpopulations in the Arctic (Letcher et al., 2010). Although concentrations of many conventional POPs have decreased in the Arctic over the past few decades due to restrictions on their production and use, levels of most POPs in biota and environmental compartment remain high since the 1990s or start to increase after years of decrease (Rigét et al., 2019;Wang et al., 2020). This is likely due to an unfavorable combination of chemical properties of POPs (e.g. high thermal stability, high volatility and slow degradation), marine currents and air flows transporting POPs northwards and changing marine food webs due to climate change (AMAP, 2004;Cabrerizo et al., 2018;Laender et al., 2011;Lie et al., 2003;Macdonald et al., 2005;US EPA, 1975;Wania, 2003). Additionally, global warming may cause POPs and mercury deposited in sea water and ice to revolatilize into the atmosphere (AMAP, 2011;Ma et al., 2011). Overall, these processes resulted in POP levels being relatively high in Arctic regions, with concentrations in mammals often exceeding threshold levels for physiological and toxicological effects (Letcher et al., 2010;Sonne, 2010). Previous research has focused on especially Arctic marine mammals (including polar bears, beluga whales and ringed seals), because POP levels are elevated in top predators; the fact that the marine environment accounts for a large percentage of area of the Arctic; and, the fact that many marine mammals are important components of the human diet (De Wit et al., 2004).
POPs and their metabolites, as well as mercury (Hg) have shown to cause adverse health effects in mammals, including disruption of reproductive, thyroid and stress hormone systems, decreasing bone density and immunological effects (Colborn, 2004;Horri et al., 2018;Lair et al., 2016;Nuijten et al., 2016;Rattner, 2009). Although polar bears are able to metabolize several POPs, the metabolites may pose even more severe negative health effects than the parent compounds (Andersen and Aars, 2016;Gutleb et al., 2010;Letcher et al., 2010). While many studies have focused on acute effects of POPs and mercury on individual health effects in polar bears (Oskam et al., 2004;Sonne, 2010), there is a lack of studies focusing on the direct link between POP concentrations and population vital rates. Bechshoft et al., 2018 found that 98% of the papers included in their review dealt with individual health effects of contaminant exposure only, disregarding population level effects (Bechshoft et al., 2018). Especially chronic exposure may give rise to subtle individual effects affecting population dynamics, but toxicological data on chronic exposure are often lacking Forbes et al., 2016;Kannan et al., 2000;Nuijten et al., 2016). To correctly assess adverse effects of POPs on polar bear populations, endpoints included should be relevant to populations, such as survival, reproductive success and population density (Nuijten et al., 2016).
Ecotoxicity endpoints for polar bears are currently grossly lacking in literature. Previously, several methods have been developed to estimate these endpoints for untested species based on known ecotoxicity data (Forbes et al., 2016;Pavlova et al., 2016;Raimondo et al., 2007). In one common method, ecotoxicity endpoints (e.g. EC 50 and LC 50 ) for untested wildlife species are extrapolated from ecotoxicity data from tested species, based on body size (allometric) scaling using chemical-specific scaling factors based on empirical research. This method in general assumes larger species to be less sensitive to persistent pollutants than smaller organisms (Sample and Arenal, 1999). Body size scaling factors, however, vary considerably among chemicals and for many chemicals these factors are lacking, in which case typically an average is applied to extrapolate toxicity among species (Mineau et al., 1996;Raimondo et al., 2007). Additionally, these scaling factors are typically developed for acute toxicity only and their applicability for estimating chronic toxicity data (and thus population effects) is unknown. As for many compounds, toxicological modes of action differ between acute and chronic effects, scaling factors may also considerably differ (Sample and Arenal, 1999). Allard et al. (2010) therefore recommend to not use allometric scaling when extrapolating endpoints between species or when extrapolating acute to chronic endpoints (Allard et al., 2010). Alternatively, species sensitivity distributions (SSDs) based on chronic ecotoxicity endpoints related to survival and reproduction may be used to estimate ecotoxicity of an untested species relative to ecotoxicity data of similar species (Raimondo et al., 2007). A lot of discrepancy exists in literature regarding the relative sensitivity of warmblooded predator species towards long-term chemical exposure, compared to other species (Shore and Rattner, 2001). Therefore, in the present study, we simulated changes in population growth rate by assuming polar bears to be moderately sensitive to POPs and mercury compared to other endothermic species, by taking the HC 50 (the hazardous dietary concentration at which 50% of the species in the SSD is affected). Additionally, changes in population growth rates were modelled in a high-risk scenario and a low-risk scenario, by taking the HC 5 and HC 95 of the SSDs.
In the present study, we investigated the potential impact of POPs on polar bear populations by using a modelling approach as developed by Hendriks and Enserink (1996). To this end, we collected POP dietary ecotoxicity data for endothermic species from the EPA's ECOTOX database (US EPA, 2019). These data were used to construct SSDs based on both EC 50 (for reproductive effects) and LC 50 data, from which the HC 50 s, HC 5 s and HC 95 s act as direct input in a model to estimate the change in intrinsic population growth rate (Hendriks and Enserink, 1996). The model was then applied to nineteen Arctic polar bear subpopulations, assuming polar bears to prey upon seal species in the same area.

Model application
High levels of POPs and mercury have been associated with multiple negative effects in marine mammals. Concentrations in Arctic marine mammals often exceed threshold levels for physiological and toxicological effects based on bioassays on rats (Letcher et al., 2010;Sonne, 2010). Strong relationships between PCBs and cortisol and sexual thyroid hormones have been observed in polar bears in the Barents Sea (Braathen et al., 2004). Although no POP toxicity data related to reproductive success of polar bears exist, data relating POP exposure to population-level effects are available for other marine mammal species (i.e. seals) in Baltic areas (Sonne et al., 2020). Therefore, we assumed POP exposure to also affect reproduction rate and survival rate, and thus population growth rate, in polar bear subpopulations. We applied our model to estimate changes in intrinsic growth rate (r(C)/r(0)) of nineteen recognized polar bear populations in different Arctic areas between 1970 and 2020 (Fig. 1). In the present study, the change in population growth rate (defined as r(C)/r(0)) is calculated according to Hendriks and Enserink (1996) and Hendriks et al. (2005) in which r(C)/r(0) is the change in intrinsic population growth rate, expressed as the ratio between the growth rate in presence of individual POPs or mercury, and the growth rate in absence of these compounds. C is the POP concentration to which the polar bear is exposed through its diet (in mg/kg w.w. marine mammal blubber), LC 50 is the half-maximal lethal concentration (in mg/kg w.w. diet), EC 50 is the half-maximal chronic effective concentration (in mg/kg w.w. diet), pertaining to reproduction, β is the slope assumed to be similar in both exposureresponse curves, ranging from −0.4 to −0.27 (median: −0.33) (Hendriks et al., 2005), and R 0 is the lifetime fecundity, defined as the total number of lifetime offspring, excluding cub survival. Lifetime fecundities in this study were calculated based on generation lengths for polar bear subpopulations separately and maximal intrinsic growth rate (r max ; Table 1), as reported by Regehr et al. (2016) and Regehr et al. (2017), through log(R 0 ) = r max *T g (Steiner et al., 2014). A r(C)/r(0) of 1 implies that there is no impact of the compound on population growth rate (i.e. the population is growing at maximum rate (100% of the maximal intrinsic growth rate)). A r(C)/r(0) of 0 implies a stable population (growth rate = 0.00), and a negative r(C)/r(0) implies population decline.
The combined toxic pressures [r(C)/r(0)] mix , defined as the change in intrinsic population growth rate of polar bear populations induced by multiple compounds per subpopulation and year are expressed as fractions of a maximum possible effect (0% ≤ E ≤ 100%). Since [r(C)/r (0)] can in theory be -∞, depending on the observed concentration, modelled intrinsic growth rate values were normalized between 1 and 0, taking the minimal modelled r(c)/r(0) across all subpopulations and years as a minimum value. [r(C)/r(0)] c values were calculated based on the response addition principle and can be calculated through (Aldenberg et al., 2002;Gregorio et al., 2013) where N i is the number of substances included in calculation of the change in population growth rate. Finally, standardized values were converted to original values. All analyses and simulations were performed in R statistics v3.5.1. (See the Supporting Information for the full R script).

Parametrization of the model
EC 50 s and LC 50 s for polar bears were estimated based on derived species sensitivity distributions for endothermic species. Dietary ecotoxicity data (reproduction-related EC 50 s and LC 50 s) for these species (including rodents, mink and chicken) were obtained from the US EPAs ECOTOX database (US EPA, 2019) (See the full dataset (.xls) in the Supporting Information). Species sensitivity distributions (SSDs) were derived based on species-aggregated medians for EC 50 and LC 50 , assuming a log-normal spread in species sensitivity (Aldenberg and Rorije, 2013;Posthuma et al., 2002). A bootstrapping method was used to determine the HC 50 , HC 5 and HC 95 of these SSDs (the chemical concentration at which 50%, 5% and 95% of the species is affected or the inflection point of the log-normal SSD). Ecotoxicity endpoints for polar bears for individual chemicals were then estimated in three scenarios, based on the HC 50 , HC 5 and HC 95 of these SSDs. In a first scenario ecotoxicity endpoints were based on the median HC 50 of the SSD (Intermediate-risk scenario), assuming polar bear sensitivity to be based on the median sensitivity to POP exposure compared to other endothermic species. In a second, high-risk scenario, estimated EC 50 s and LC 50 s for polar bears were based on the 5th percentile of the SSD (HC 5 ), assuming polar bears to be relatively sensitive to POPs. In a final (low-risk) scenario, endpoints were based on the 95th percentile of the SSD (HC 95 ), assuming polar bears to be relatively insensitive to POPs (Fig. 2). All parameters values included in the present are shown in Table 1.

Data collection
2.3.1. Exposure data Data on POP and mercury residues in marine mammal species (mainly ringed seal (Phoca hispida), spotted seal (Phoca largha), harp seal (Phoca groenlandica), ribbon seal (Crystophora cristata), bearded seal (Erignathus barbatus), walrus (Odobenus rosmarus) and narwhal (Monodon monoceros)), assumed to be the main prey of polar bears in the Arctic (Sonne, 2010), were compiled to calculate potential changes in intrinsic growth rates of polar bear populations. POP concentrations (transformed to mg/kg wet weight (w.w.)) in prey for each subpopulation area were obtained from a literature search using the Web of Knowledge and Google Scholar. As healthy adult polar bears are known to mainly eat blubber of ringed seals (subadults are known to scavenge the kills of others and will feed on muscle tissue), only residues in blubber were included in the present study (Cherry et al., 2011;Dyck and Kebreab, 2009;Stirling and McEwan, 1975). Concentrations on lipid basis were converted to wet weight (w.w.) basis, based on the reported lipid content. If no lipid concentration was reported, a lipid content of 85% was assumed for marine mammal blubber samples. To calculate the toxic equivalency of PCBs, we assumed that the planar and coplanar PCB composition in marine mammal blubber was similar across all sampled individuals. In the present study, the planar PCB composition in blubber was taken from Savinov et al., 2011(Savinov et al., 2011. Although PCB-77, −126 and −169 show to contribute substantially to total TEQs, due to their high TEF values (Table 2) (Letcher et al., 1996), monitoring data on levels of these congeners in seal blubber are lacking in literature or are measured to be below the detection limit. As the way in which we calculated TEQ is highly dependent on data availability for all congeners, toxic equivalencies (TEQs) for each blubber sample were calculated based on the ratio of most dominant dioxin-like PCBs to ∑PCBs (PCB-118/∑PCBs (59.5%) and PCB-105/ ∑PCBs (24.1%)).
All concentration data used in our simulations were collected between 1972 and 2018. In our simulation, we assumed that a polar bear consumers blubber from one seal every six days of a year and each serving consists of 25 kg of blubber (Dyck and Kebreab, 2009). r (C)/r(0)s were calculated for each data record. Subsequently, the medians of these values were taken for subsets of data based on chemical, geographical location (subpopulation) and year.

Ecotoxicity data
In the present study, we modelled changes in polar bear intrinsic population growth rate based on concentrations in fat tissue of their main prey. Effect concentrations for polar bears were based on effect concentrations for other endothermic species. Dietary ecotoxicity data (standardized to mg/kg in diet) for endothermic species (including rodents, mink and chicken) were taken from the ECOTOX database (US EPA, 2019). Only toxicity endpoints pertaining to reproduction and mortality and administered through diet (food or capsule) were taken into account. NOECs (No-Observed-Effect-Concentrations) or NOELs (No-Observed-Effect-Levels) were transformed to EC 50 s and LC 50 s according to Aurisano et al., 2019(Aurisano et al., 2019. If no ecotoxicity data on at least three endothermic species were available, these data were extrapolated based on data from the other endpoint (EC 50 , when no sufficient LC 50 data were available and LC 50 , when no sufficient EC 50 data were available), assuming EC 50 /LC 50 = 0.1 (OECD, 1990). These ecotoxicity data were used in the derivation of Species Sensitivity Distributions (SSDs) for multiple POPs in which effect concentrations are plotted against the potentially affected fraction of species and fitted by a lognormal curve, yielding a cumulative distribution (Posthuma et al., 2002). Although animals at the top of the food chain are endangered due to increased POP exposure (De Wit et al., 2004;Kallenborn, 2006), there is a lack of evidence that sensitivity increases (e.g. the EC 50 or LC 50 decreases) along trophic levels in the Arctic marine food chain, as reported mammalian ecotoxicity data in literature varies considerably between and within species . However, for aquatic ecosystems, Baird and Van den Brink (2007) suggest that, for most chemicals, LC 50 might be higher for predators, albeit reported pvalues being above 0.05 (Baird and Van den Brink, 2007). Jeram et al., Table 1 Parameters used in the equations with typical or default values.

Model performance
Although intrinsic population growth rates are available for some polar bear subpopulations and years (Hunter et al., 2010;Lunn et al., 2016), these data are too limited to establish reliable population trends. Therefore, modelled changes in growth rates were compared to current population sizes. These population sizes were based on data on the most   (Durner et al., 2018) and additional population size data based on mark-recapture analysis from literature (Aars et al., 2009;Amstrup, 1995;Lunn et al., 2016;Regehr et al., 2018;Stirling et al., 2011), and were supplemented with simulated population trends by York et al. (2016)  , based on mark-recapture and aerial survey data for polar bears collected by multiple authors (Obbard et al., 2007;Obbard et al., 2013;Peacock et al., 2012;Peacock et al., 2013;Polar Bear Technical Committee, 2007;Regehr et al., 2007a;Regehr et al., 2007b;Sodhi and Ehrlich, 2010;Stapleton et al., 2016;Taylor et al., 2002;Taylor et al., 2005;Taylor et al., 2008a;Taylor et al., 2006a;Taylor et al., 2008b;Taylor et al., 2009;Taylor et al., 2006b), manually digitized using DigitizeIt (See the Supporting Information for the complete dataset). Temporal trends in model outcomes were compared to trend data for polar bear subpopulations as reported by the Polar Bear Specialist Group (IUCN/SSC Polar Bear Specialist Group, 2019). Additionally, temporal trends in modelled changes in population growth rates for individual polar bear subpopulations were quantitatively compared to relative population sizes (scaled between 0 and 1) for each subpopulation individually, as well as to median relative population size across all subpopulations, due to a severe lack of population trend data. Model performance was evaluated by means of R-squared in R statistics v3.5.1.

Species sensitivity distributions
Dietary species sensitivity distributions (SSDs) were constructed for endothermic species for four compounds and toxic equivalencies based on Aroclor mixtures (Fig. S2). Polar bear effect concentrations (EC 50 and LC 50 s) were then based on the HC 50 (the median or infection point of the SSD), HC 5 and HC 95 (low-and high-risk scenarios, See Table 3). The lowest EC 50 and LC 50 , and thus highest toxicity, was calculated for PCBs (expressed as TEQ dioxin), followed by Mercury (Hg) , o,p'-DDT and p,p'-DDT.

Spatial effects: combined effects
Considerable variation in combined r(c)/r(0) could be observed among the three scenarios (Table 4) (Fig. 3). Highest combined r(c)/r (0) values, and thus smallest effects of POPs, were calculated for Gulf of Boothia (GB, 0.92) and Foxe Basin (FB, 0.85). In our high-risk scenario, all of the 15 subpopulations yielded a negative combined r(c)/r(0), implying that for these subpopulations POP levels in blubber of ringed seal severely affected polar bear populations to such extent that a decrease in population size could be expected. Finally, in our low-risk scenario,  assuming polar bears to be relatively insensitive to POP exposure, for ten subpopulations negative average combined r(c)/r (0)

Spatial effects: individual compounds
When looking at individual chemicals, especially PCBs (expressed in toxic equivalency (TEQ)) showed to contribute to high combined effects, with individual r(c)/r(0) values yielding values below zero for four out of 19 polar bear subpopulations. Next to PCBs, levels of the DDT metabolite p,p'-DDT in marine mammal blubber also showed to be potentially hazardous to five out of 19 polar bear subpopulations, with r(c)/r(0) values around or below 0, indicating stagnation of the population size (Fig. 4).

Temporal effects
No significant temporal trend was observed in both the modelled changes in intrinsic growth rate and relative monitored polar bear population size (subpopulation mean and pooled standard deviations), with R 2 s of 0.17 and 0.04, respectively (Fig. 5). Additionally, no statistically significant correlation between median modelled changes in intrinsic growth rates and median relative monitored subpopulation size was observed (R 2 < 0.1). A statistically significant temporal trend could only be determined for r(c)/r(0) values calculated for p,p'-DDT (R 2 = 0.6, p < 0.01, Fig. 6). Yet, modelled changes in intrinsic population growth rate and monitored relative population sizes for all individual chemicals were not correlated (0.00024 > R 2 < 0.16, p > 0.1, Fig. 6).
In the present study, we investigated the possibility of simulating changes in polar bear population growth rates due to pollution with persistent organic pollutants, as levels of these compounds in Arctic biota remain high, and are negatively correlated to polar bear population density (Nuijten et al., 2016). We therefore expected a decrease in potential population growth rate for the four compounds included in this study (Hg, o,p'-DDT, p,p'-DDT and p,p'-DDD) and toxic equivalencies based on Aroclor mixtures for the nineteen recognized polar bear subpopulations. As chronic ecotoxicity data for polar bears were grossly lacking, several assumptions were made, as commonly applied in chemical risk assessment. (Table S1 in the SI). First, in evaluating spatial differences in risks imposed by POP exposure, we assumed POP levels to remain constant over time . Secondly, toxicity levels, i.e. reproduction EC 50 s and survival LC 50 s for polar bears, were based on toxicity data for endothermic species (including mink, rodents and bird species) (US EPA, 2019). EC 50 and LC 50 values for polar bears were set at the HC 50 , HC 5 and HC 95 of the SSDs. Sensitivity of marine mammals to POPs has been shown to differ across species, locations and chemicals, with PCBs showing to pose the greatest risk for polar bear populations Nuijten et al., 2016). Multiple allometric relationships have been developed for acute ecotoxicity, typically following the principle of larger organisms yielding higher endpoint values (and thus lower sensitivity), as in larger organisms reaching an equilibrium concentration typically takes more time (Hendriks, 1995). However, allometric relationships for chronic ecotoxicity are lacking (Sample and Arenal, 1999). So, as an alternative we used three separate scenarios in the present study. HC 50 s (assumed to be the EC 50 s and LC 50 s for polar bears) for dietary SSDs for compounds derived in the present study were typically 2 to 10 times higher than HC 50 s (or μs), based on dietary No-Observed-Effect-Concentrations (NOECs), calculated for the same compounds (or TEQ in lipid weight) as reported by Schipper et al. (2010) and Korsman et al. (2016) (See Table 4) (Korsman et al., 2016;Schipper et al., 2010). Additionally, the calculated HC 50 for DDT based on dietary LD 50 data for mammals (rats and mice) and birds (chicken, mallard and wild birds) from multiple studies (ATSDR, 2002;Gaines and Medicine, 1960;Gaines and pharmacology, 1969;Hudson et al., 1979;Luttik et al., 1997;Mineau et al., 2001;Schafer et al., 1983;Schafer and Bowles, 1985), calculated by Golsteijn et al., 2012 showed to be 1.5 times higher than the HC 50 calculated for p,p'-DDT in the present study (Golsteijn et al., 2012). However, again, these HC 50 s were based on acute ecotoxicity data, rather than chronic ecotoxicity data. SSD derivation is typically based on limited ecotoxicity data for only a small amount of test species, leading to increasing uncertainty in HC 50 values (Raimondo et al., 2007). However, despite the differences in test species and effects, all HC 50 values reported in literature were within the 95% confidence intervals of the HC 50 s from SSDs reported in the present study.
Thirdly, another assumption, related to the ecotoxicity of POPs in polar bears, is that the range of the slopes (β) of the exposureresponse curves for polar bears in the Monte Carlo iterations was similar for both reproduction (EC 50 ) and mortality (LC 50 ). The steepness of the slope of the exposure-response curve is of high importance in chemical risk assessment, as a steeper slope indicates higher percentages of individuals within a population to be affected by the chemical (Penningroth, 2016). Mortality and reproduction rates used in the present study are assumed to be equally important for changing population growth rates (Eq. (1)). Therefore, an increased steepness in slope in exposure Fig. 4. The ratio of exposed and control population intrinsic rate of increase r(C)/r(0) at several Arctic locations, based on the intermediate scenario in which we assumed polar bears to be among 50% most sensitive endothermic species. Risks were based on monitored POP concentrations in seal species (ringed seal (Phoca hispida), spotted seal (Phoca largha), harp seal (Phoca groenlandica), ribbon seal (Crystophora cristata), bearded seal (Erignathus barbatus), and walrus (Odobenus rosmarus).  response curves for either mortality or reproduction may lead to an equally increasing negative effect on polar bear subpopulations. Although thorough research on differences in slopes between ecotoxicity endpoints is lacking, multiple studies suggest that the heterogeneity in slopes of exposure-response curves for the same compound may be induced by differences in experimental design, such as temperature or exposure duration (Lenters et al., 2011;Samoli et al., 2005). All assumptions in the present study had to be made due to a lack of ecotoxicity and monitoring (POP residues in blubber of prey species) data. This approach would greatly improve by including more data on POP ecotoxicity and exposure, or by incorporating new techniques for estimating chronic ecotoxicity by using e.g. interspecies correlations.

Spatial trends
Lowest combined r(c)/r(0) values, and thus largest effects of POPs, were calculated for Baffin Bay and Kane Basin in the Canadian Arctic, and Eastern Greenland, Barents Sea, and Kara Sea in the Eurasian Arctic. Although r(c)/r(0) values for Baffin Bay, the Barents Sea and Kane Basin were based on a considerable amount of chemical residue data in prey species (Table 4; >100 data records), residue data used in modelling changes in population growth rates for the Kara Sea and Eastern Greenland were relatively scarce (<50 data records). Considerable variation in combined effects were observed among the three scenarios, but the order of vulnerability of the subpopulations remained the same across the three scenarios. As r(c)/r(0) is directly related to the LC 50 s and EC 50 s of the individual chemicals, this may imply that the composition of PCBs, DDT metabolites and Hg in blubber and liver of prey species is similar across all subpopulations. In this case bias may be introduced as we assumed both PCB congener composition in blubber of prey species, and PCB composition in Aroclor mixtures to be similar across all subpopulations and bioassays in calculating toxic equivalency. However, as we use concentrations of dioxin-like PCBs to calculate TEQ values, we expect variation in TEQ, and thus variation in modelled r (c)/r(0) to be relatively small. Overall, considerable differences in modelled combined r(c)/r(0) values were observed between subpopulations (Fig. 3). This may be due to differences in trophic interactions among regions or to different polar bear feeding strategies (Kleivane et al., 2000;McKinney et al., 2009). Although Svalbard, Eastern Greenland and Hudson Bay are known as the hotspots in terms of POP levels in polar bear fat tissue Letcher et al., 2010), Hamilton and Derocher (2019) identify the two Beaufort Sea polar bear subpopulations and the Arctic Basin subpopulation as the most vulnerable, due to the small shelf area, low prey diversity and the consequences of climate change (Hamilton and Derocher, 2019).
Polar bears occur in relatively discrete subpopulations (SWG, 2016). When comparing modelled risks imposed by chemical compounds between subpopulations, we again made a couple assumptions regarding the ecology of polar bears and prey species. First, we assumed that polar bears hunt within their assigned subpopulation area. Secondly, we assumed prey species to reside within these areas, while earlier work shows that, although polar bears tend to stay in discrete subpopulations, movement between these populations is not uncommon (IUCN/SSC Polar Bear Specialist Group, 2019). Given the large geographical home range sizes for polar bears (Auger-Méthé et al., 2016), this may especially be the case for smaller subpopulation areas. Moreover, home ranges of ringed seals, the polar bears most abundant prey species, have shown to be notoriously seasonal (Kelly et al., 2010). A third assumption that may cause bias in model outcomes is the fact we assume that the diet of polar bears is similar across all subpopulations, consisting of mainly one species group, namely marine mammals. Polar bear diet is known to vary considerable spatially and throughout seasons due to sea-ice decline in summer (Gormezano and Rockwell, 2013), which may have major implications for POP accumulation and exposure (De Laender et al., 2009).

Temporal trends
We calculated a decreasing trend in potential effects of POPs on polar bear population growth rate, based on POP and mercury levels in the blubber of prey species, albeit it not being statistically significant. This insignificance is likely due to a combination of insufficient POP and mercury residue data in biota, and insufficient data on monitored population sizes. No statistically strong correlation could be observed between modelled changes in population growth rates and median monitored population sizes (Figs. 5 and 6). The chemicals included in the present study, however, cover only a small fraction of contaminants potentially affecting polar bear population growth in the Arctic. Therefore, full validation of this model shows to be challenging. Concentrations of many POPs have decreased in the Arctic over the past few decades, due to restrictions on their production and use. Subsequently, the composition of the mixture of POPs to which Arctic biota is being exposed has substantially changed since the 1980s, resulting in decreased contribution of DDTs and its metabolites to the total risk imposed by this POP mixture (Villa et al., 2017). This decrease is likely less prominent for PCBs due to difficulties in controlling PCB emissions (Villa et al., 2017). This is concerning, as PCB exposure still shows to negatively affect polar bears (Sonne, 2010). Earlier work by Rigét et al. (2019) shows that an increasing long-term monitoring trend of certain legacy POPs in Arctic biota may be observed in a fraction of locations and species, including β-HCH and PBDE-47 (Rigét et al., 2019). Additionally, emerging persistent pollutants for which ecotoxicity data are lacking, such as perfluoroalkyl substances, are recognized as potential threats for the future (Strobel et al., 2018;Tartu et al., 2018;Villa et al., 2017).
Other anthropogenic stressors may negatively affect polar bear subpopulations in the future. Multiple studies have focused on the effects of global warming-induced sea ice decline on polar bear populations, as this potentially restricts polar bears from hunting on the ice (e.g. (Amstrup et al., 2008;Bromaghin et al., 2015;Hamilton and Derocher, 2019;Jenssen et al., 2015)). Sea ice decline leads to increasing periods of fasting, with polar bears arriving on land earlier each year, likely resulting in decreasing reproductive success and adult survival (Laidre et al., 2020;Molnár et al., 2020). However, earlier work by Nuijten and all showed that PCB and DDT concentrations explained more of the variance in polar bear density than ice cover and hunting (Nuijten et al., 2016). Next to population decline due to starvation, this may also affect POP levels in polar bear adipose tissue. Climate-change induced sea-ice decline may induce changes in the polar bears diet, and therefore may result in reduced fat tissue stores, and increasing levels of lipophilic compounds, such as POPs (Tartu et al., 2017). Additionally, this may affect the degree of bioaccumulation of POPs in Arctic biota (Dietz et al., 2004;Hoondert et al., 2020). Global warming may also cause legacy and emerging POPs deposited in sea water and sea ice to revolatilize into the atmosphere . Revolatilization and subsequent bioavailability of POPs and mercury may also be enhanced by an increasing frequency of global warming-induced extreme events such as melting ice, storms, floods, and forest fires (Teng et al., 2012). Earlier work by Hung et al. (2010) already showed an increasing trend in atmospheric PCB concentrations in the early-to-mid 2000s, potentially due to revolatilization from the ocean due to global warming, at two air monitoring stations; Zeppelin on Svalbard, and Alert in the Canadian Arctic (Hung et al., 2010). Although Rigét et al. (2010) observed an increasing temporal trend in both air temperature and Hg concentrations in Arctic char in a land-locked Greenland lake, concentrations of PCBs and DDTs showed to decrease over time. However, top predators in Arctic food webs (i.e. polar bears) are significantly impacted by climate change, with global-warming-induced sea ice decline hindering these species in finding prey (i.e. seals) on the ice (Jenssen et al., 2015). Increasing longer periods of fasting may result in increasing POP levels, which may lead to an increased mortality and decreasing reproductive success (Jenssen et al., 2015). A final anthropogenic stressor that has impacted polar bear subpopulations in the past is subsistence hunting by indigenous people. This showed to be especially harmful in the East Greenland polar bear subpopulation. However, due to the introduction of hunting quotas and a decrease of hunting because of increasing climate-change induced unpredictable weather events, hunting activities have decreased significantly in this area, allowing the population to increase . Due to these quotas and extreme weather events, hunting will likely play a less significant role in polar bear population decline in other areas in the future (Born et al., 2011).

Spatio-temporal trends
In our diet-based modelling, PCBs (expressed in toxic equivalency (TEQ)) and p,p'-DDT contributed most to the total effect on intrinsic population growth rate. Based on observed polar bear adipose tissues, Nuijten et al. (2016) and Sonne et al. (2009) concluded that concentrations of PCBs and DDTs in polar bear adipose tissue exceed critical levels for liver lesions and testes size in rat (Nuijten et al., 2016;Sonne et al., 2009). Work by Sonne et al. (2009) concluded that polar bears in the Barents Sea are most likely affected by PCBs (Sonne et al., 2009). When looking at spatio-temporal trends, we did not find a strong correlation between our modelled r(c)/r(0) values and observed monitored population trends for the majority of chemical-subpopulation combinations. This is likely due to an insufficient amount of data on capture-recapture data for polar bear subpopulations, in combination with a lack of corresponding data related to POP concentrations in seal blubber. Additionally, in the present study we compared change in population growth rate with monitored relative population sizes. As it takes some time for a population size to reflect a change in population growth rate, it is not unrealistic to expect a delay in r(c)/r(0) compared to the monitored relative population sizes. Trends in modelled population growth rates reported in the present study were inconsistent with monitored population level trends for polar bear subpopulations as reported by the Polar Bear Specialist Group for the Northern and Southern Beaufort Sea and the Southern and Western Hudson Bay (IUCN/SSC Polar Bear Specialist Group, 2019). While modelled population growth rates in this study tend to increase over time, monitored population sizes have decreased for these polar bear subpopulations., which may suggest incorrect estimation of polar bear sensitivity. Another explanation may be that the decreasing risk of POPs in the Arctic may be counterbalanced by increasing importance of other emerging POPs (e.g. organophosphate pesticides or perfluoroalkyl compounds) (Routti et al., 2019;Sühring et al., 2016) or risks imposed by other anthropogenic factors, such as global warming or hunting. A recently published paper by Molnár et al. (2020) modelled decline in cub recruitment and adult survival in subpopulations subjected to prolonging periods of fasting due to sea-ice decline. In this study, Canadian subpopulations at lower latitudes (Western and Southern Hudson Bay, Davis Strait, Foxe Basin and Baffin Bay) were identified as those most sensitive to climate change induced sea ice decline, while not all of these subpopulations showed to be severely affected by POP exposure. When combining the results of Molnár et al. (2020) with our results, we conclude that especially polar bears in Baffin Bay as well as in the Southern Hudson Bay will be severely affected by both POP exposure and sea ice decline. In contrast, subpopulations residing at higher latitudes (e.g. Norwegian Bay) will likely be less affected by both climate change and POP exposure.

Conclusion
To our knowledge, we were the first to investigate the potential of modelling changes in polar bear population growth rate due to POP exposure. We found considerable variation in effects across polar bear subpopulations, with blubber of prey species in multiple subpopulations yielding POP levels high enough to result in polar bear population decline. However, a decreasing temporal trend in POP levels in seal blubber, and thus an increase in potential polar bear population growth rate was observed for TEQ (PCBs) and DDT metabolites, implying that other anthropogenic factors (especially sea ice decline due to global warming) plays an important role in decreasing polar bear populations. In the analysis several assumptions, common in risk assessment, were made such as chronic ecotoxicity data for warm-blooded predators, monitoring data of POPs in prey species, data related to polar bear ecology and especially population size data to validate model outcomes were all grossly lacking. Further research on these subjects is required to further improve models used in estimating population levels effects imposed by POP exposure.

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