Breeding den selection by Arctic foxes (Vulpes lagopus) in southern Yamal Peninsula, Russia

Selecting the right location for a den during the breeding season is a type of habitat selection in the Arctic fox (Vulpes lagopus) that is likely to affect its reproductive success. A den’s suitability likely depends on its ability to provide shelter, as well as its proximity to prey resources. Depending on the different relative risks that Arctic foxes may face across their broad circumpolar range, Arctic foxes may place different emphases on selection for shelter and prey resources in different ecosystems. Understanding the different requirements for reproduction under different ecological conditions is highly relevant to conservation efforts in areas where Arctic foxes are threatened by rapid environmental changes. Here, we investigated the relative selection for shelter and prey resources in southern Yamal Peninsula (Russia) using data from 45 dens collected over a 13-year period. Arctic foxes preferred to breed in dens with more den entrances; an indicator of shelter quality. Arctic foxes also preferred dens surrounded by more prey resources (quantified by the amount of river valley habitat), but this result was less conclusive. These results complement the findings reported from other study areas, illustrating that Arctic foxes in ecosystems with diverse predator communities may put emphasis on selection for shelter quality. In less productive ecosystems, Arctic foxes may rather put emphasis on selection for prey resources. As tundra ecosystems become more productive and generalist predators move north, the reproductive requirements and habitat selection of Arctic foxes may change accordingly, depending on the species’ ability to adapt.


Introduction
Habitat selection is the disproportionate use of available habitat types to positively influence the survival and reproduction of an individual (Block and Brennan 1993). Selecting the right location for a nest or den during the breeding season is one type of habitat selection that is likely to affect reproductive success (Doligez and Boulinier 2008). The breeding period is a critical time for many animals, because of the need to shelter their vulnerable and often immobile offspring and provide access to resources for the growing young. Understanding how an animal selects its breeding den can therefore provide insight into a species' requirements for reproduction and the relative importance of different limiting factors in a specific ecological context.
The Arctic fox (Vulpes lagopus) is a small predator found over a wide circumpolar range in Arctic and alpine tundra habitats (Audet et al. 2002). Throughout its range its dietary niche varies from populations being specialised on lemmings (e.g. Scandinavia) to being generalist and relying on marine resources (e.g. Svalbard) (Elmhagen et al. 2000;Eide et al. 2012). In some areas, the species is threatened due to historic hunting pressure (Lönnberg 1927), a northward expansion of the red fox (Vulpes vulpes) and changes in small rodent dynamics (Hersteinsson and MacDonald 1992;Ims and Fuglei 2005;Killengreen et al. 2007;Angerbjörn et al. 2013;Elmhagen et al. 2017). As an endemic species to the Arctic tundra, it has been selected as one of ten "flagship species" expected to suffer most from future climate change (IUCN 2009).
Arctic foxes depend on dens for shelter during the breeding season (Eberhardt et al. 1983;Reichman and Smith 1990). Due to difficulties of digging in frozen Arctic soil, these dens are often located on southerly exposed slopes and mounds with little snow accumulation (Smits et al. 1988;Prestrud 1992a;Nielsen et al. 1994;Dalerum et al. 2002;Szor et al. 2008). Dens can be up to several hundred years old (Macpherson 1969). Interestingly, only 30 to 50% of the available dens is actually used by breeding Arctic foxes (Anthony 1996;Dalerum et al. 2002;Szor et al. 2008;Gallant et al. 2012;Ehrich et al. 2017), suggesting that not all dens are equally preferred or suitable for breeding (Tannerfeldt et al. 2003). Suitable breeding dens may thus be a limiting resource in some populations (Smits and Slough 1993;Eide et al. 2012).
In addition to providing shelter, the location of a den within the surrounding habitat affects distances to prey resources and feeding efficiency for central place foragers (Orians and Pearson 1979). The suitability of a den for breeding is therefore expected to depend on its ability to provide shelter against predators and environmental conditions, as well as on its proximity to surrounding prey resources, which may fluctuate from year to year (Tannerfeldt et al. 2003). Foxes, in their selection of breeding dens, thus have to make choices in the relative emphasis they place on minimising predation risk, minimising energy expenditure in unfavourable environmental conditions, and minimising starvation risk. The relative emphasis placed on these factors may change with yearly fluctuating abundances of prey and predators, affecting the relative risks of predation and starvation.
On Bylot Island in the high Arctic, proximity to prey resources was shown to determine breeding den selection, in that Arctic foxes either denned close to a goose colony or in better lemming habitat (Szor et al. 2008). During years with high lemming abundance, proximity to the goose colony was less important to denning Arctic foxes (Giroux et al. 2012). In low Arctic and Scandinavian study areas with red foxes and other predators like wolverines (Gulo gulo), Arctic foxes preferred dens with more entrances which are likely to provide better shelter (Dalerum et al. 2002;Gallant et al. 2014). This illustrates the different emphases that Arctic foxes may place on minimising different relative risks in different environments.
Better understanding of the different relative risks that breeding Arctic foxes may face in different environments may benefit management and conservation efforts in areas where Arctic foxes are threatened. Because Arctic foxes are found over a wide circumpolar range, studies conducted in differing ecosystems should contribute to this aim. As proposed by Gallant et al. (2014), studies conducted in ecosystems with differing risks from large predators, differing environmental conditions, and differing abundances of prey resources, could help explore the relative importance of a den's ability to provide shelter versus the proximity to preyrich habitats.
In this study, we investigate breeding den selection by Arctic foxes in southern Yamal peninsula, Russia. Our study area is interesting because it represents an intermediate situation compared to previous studies. Namely, in terms of starvation risk, it is located in the low Arctic, so it is more productive in terms of prey than high Arctic study areas. In terms of predation risk, however, red foxes are relatively scarce compared to other low Arctic study areas (Sokolov et al. 2016). Still, the number of other potential predators like ravens (Corvus corax), birds of prey, wolverines, dogs (Canis lupus familiaris), and humans is greater than in high Arctic study areas. Using a 13-year dataset on the use of 45 dens, we evaluate the three competing hypotheses that Arctic foxes place greater emphasis on selecting for (1) shelter quality, (2) proximity to prey resources, and (3) place emphasis on selection for both shelter and prey.
Because our study area is likely to represent an intermediate situation compared to previous studies in terms of predation and starvation risk, we predict that den selection in our study area is best explained by taking into account both the quality of the shelter and the proximity to prey resources (hypothesis 3). Yearly fluctuating abundances of small rodents (voles and lemmings), a key prey resource in our study area, will likely affect the relative risks of starvation and predation on a yearly basis. Accordingly, during years with a high abundance of small rodents (and a lower starvation risk) we predict a relatively greater emphasis on selection for shelter quality. During years with a low abundance of small rodents (and a higher starvation risk) we predict a relatively greater emphasis on selection for prey resources.

Study area
The study area is located at the Erkuta Tundra Monitoring Site in the southern part of Yamal Peninsula, Russia (68.2° N, 69.2° E) (hereafter Erkuta). Lying on the border between erect dwarf shrub tundra and low shrub tundra (Walker et al. 2005), Erkuta's topography is mostly flat, with hills up to 40 m high. Dense thickets of willow (Salix spp.) and sometimes alder (Alnus fruticosa) grow along numerous rivers and lakes. In low-lying areas, spring flooding can be widespread. Occasional steep sandy cliffs and slopes are found along riverbanks, lakes and hillsides (Sokolov et al. 2012;Ehrich et al. 2012). The mean temperature is − 24.9 °C in January and 14.5 °C in July (World Meteorological Organization 2019), and a stable snow cover is usually present between early October and early June (Sokolov et al. 2016).

3
The soil is predominantly sand and clay and provides good conditions for digging dens; permafrost is continuous (Pavlov and Moskalenko 2002).
Before the turn of the century, lemmings (Lemmus sibiricus and Dicrostonyx torquatus) were common in the study area, and their population size was characterised by high amplitude fluctuations every 3-5 years (Shtro 2009). The last lemming peaks were recorded in 1999 and 2004 (Fufachev et al. 2019). Currently, the small rodent community is dominated by voles (mainly Microtus gregalis and M. middendorffii) showing low amplitude population fluctuations (Sokolova et al. 2014). Medium sized herbivores include willow ptarmigan (Lagopus lagopus), mountain hare (Lepus timidus), and muskrat (Ondatra zibethica) ). In addition, 40 species of passerines, waders, and waterfowl breed in the study area (Sokolov et al. 2012). Nomadic Nenets reindeer-herders use the study area as pastures for their reindeer (Rangifer tarandus) herds in all seasons (Sokolov et al. 2016).

Habitat map
A habitat map of the study area was made by doing a supervised classification of a Landsat-8 satellite image, courtesy of the U.S. Geological Survey, using the QGIS SCP plugin (QGIS Development Team 2016). The image has a pixel size of 30 by 30 m and was downloaded from earthexplorer. usgs.gov (scene identifier; LC81650122016227LGN01, acquisition date; 2016/08/14). The habitat composition of the study area was categorised into eight different habitat types ( Table 1). The habitat types were defined as clearly recognisable landscape elements which were considered likely to represent different habitat qualities for herbivores. More detailed descriptions and photographs of the habitat types is given in Online Resource 1. A detailed description of the habitat map classification and accuracy assessment is given in Online Resource 2.
In addition to the aforementioned eight habitat types, the study area was divided into river valley (land below 6 m altitude) and upland (land above 6 m altitude) using the ArcticDEM dataset in ArcGIS (Porter et al. 2018). This approximately corresponds to upland and lowland mapped in the study area by Sokolov et al. (2012). All willow thicket pixels adjacent to other habitat types except water were classified as "willow thicket edge" habitat. Observations showed that these linear edge habitats were used by hunting Arctic foxes, in correspondence with observations of linear habitat use by other canids (Dickie et al. 2017;Bischof et al. 2019).

Prey species abundance indices per habitat type
Arctic fox diet in Erkuta consists mostly of small rodents and birds, particularly waterfowl but also waders and passerines ). There is a variety of other potential prey species that are common in the study area including hare and ptarmigan, which are also regularly used by the foxes. Abundance indices of small rodents per habitat type were created using trapping data and abundance indices of hare and ptarmigan per habitat type were obtained from faeces counts, as described below. Waterfowl and other birds were assumed to be generally more abundant in the river valley (Sokolov et al. 2012).

Small rodents
Small rodents were trapped from 2007 to 2019 following the small quadrat method of Myllymäki et al. (1971) as part of the long-term monitoring programme in Erkuta (see Ehrich et al. 2017 for more details). For the prey species abundance indices, we used the years 2012-2019, during which trapping was conducted on 54 small quadrats in three parts (spatial units) of the study area. The trapping quadrats were originally placed in three different landscape elements: "wet tundra", "dry tundra", and "willow thicket edge" (Sokolova et al. 2014). For use in this study, these quadrats were revisited and categorised into the habitat types used in the habitat map (Table 1).
No trapping quadrats were located in "birch shrub". Here, small rodent abundances were estimated from trapping with a trapline method from 2010 to 2015 in the centre of the study area (Fufachev et al. 2019). To account for differences with the quadrat method, the number of rodents estimated with the trapline method was multiplied by 0.67 following Fufachev et al. (2019) and Fauteux et al. (2018). Online Resource 3 shows an overview of the trapping effort used in this study and its distribution over the habitat types. Abundance indices of small rodents were calculated by taking the mean number of trapped rodents per 100 trap nights per habitat type.

Hare and ptarmigan
Hare and ptarmigan faeces have previously been used as indices of abundance and habitat use (e.g. Ims et al. 2007;Henden et al. 2011a). Faeces were counted within 50 × 50 cm squares grouped into transects of 20, in late July 2017. In each transect, squares were placed two metres apart in two parallel rows of 10, which were placed 10 m apart. Faeces counts were conducted on three transects per habitat type in three separate parts of the study area (spatial units) located between two and five kilometres apart. This resulted in a total of nine transects per habitat type "wet lowland", "mesic tundra", "dry tundra", "willow thicket edge", and "birch shrub".
The locations of transects were chosen by randomly selecting 3 out of 10 possible locations per unit per habitat type using the habitat map. Possible locations were selected based on a relatively large and homogeneous cover of the habitat type in question. When the location proved unsuitable upon arrival, transects were moved to the nearest suitable site instead. All sampled transects of the same habitat type were located 350 m or more apart.
Only intact faeces that were not in advanced stages of degradation were counted. Because faeces were only sampled once, the counts represent the accumulated relative habitat use over several seasons.

Arctic fox dens
As far as possible, all known dens were visited yearly between late June and late July. Upon finding evidence of reproduction, an automatic camera or visual observations were used to ascertain the presence of pups. Moving of pups between dens has been observed three times. Because these movements can be caused by disturbance, the first den in which pups had been observed was taken to be the selected breeding den.
Since 2007 the study area has been gradually enlarged from 130 to 230 km 2 , increasing the number of known dens from 12 (2007) to 60 (2019). Since 2007, reproduction has been observed in 26 different dens, the smallest of which had four entrances. Dens with more than 3.5 entrances on average (measured over several years) were thus considered to be large enough for reproduction (45 dens in total). Additional shelter quality characteristics were measured on 38 dens large enough for reproduction in 2017.

Shelter quality
Shelter quality may be described by several variables. A greater number of den entrances may enable more escape routes from predators (Dalerum et al. 2002;Gallant et al. 2014). The number of open den entrances was measured during most yearly visits and an average was taken (variable: entrances). A den's topographical position in the open tundra landscape may also affect the quality of shelter. A southerly slope exposure may provide a favourable microclimate (Tannerfeldt et al. 2003). On the other hand, dens situated below steep slopes can be expected to accumulate more snow, hindering access to the den at the start of the breeding season. Southerly slope exposure may also affect snowmelt. Slope exposure was measured using a compass with dens on flat hilltops counted as exposed to all directions (variable: south exposure). The slope angle above the den was measured using a Nikon Forestry PRO Laser Hypsometer (variable: upper slope).

Prey indices
A prey species index per den was calculated by multiplying the cover of surrounding habitat types by their respective prey species abundance indices. The resulting value was then scaled to a mean of 0 and a standard deviation of 1 over the different dens.
Arctic fox summer and spring home range sizes described in literature range from 4 to 60 km 2 Prestrud 1992b;Angerbjörn et al. 1997;Anthony 1997;Landa et al. 1998;Eide et al. 2004;Poulin et al. 2021). This corresponds to a circular radius of 1.1 to 4.4 km from the den. In our study area, home ranges per den overlap greatly when considering a radius of more than 1.5 km. Furthermore, the core area of the home range is arguably used more intensively than the margins. Following Gallant et al. (2014), using a 1.5-km home range radius was thus considered most appropriate to keep habitat measures per den are as independent as possible whilst maintaining a biologically plausible size. The cover of surrounding habitat types within this 1.5-km radius from each den was calculated using the habitat map and the ArcGIS Buffer and ZonalHistogram tools (ESRI 2017).

Yearly fluctuating prey resources
The yearly abundance of small rodents and mass reindeer mortality have previously been shown to determine the breeding activity of Arctic foxes in our study area ). These yearly conditions were thus taken into account in the analysis. The yearly abundance of small rodents in spring (June) was estimated by snap trapping as part of the long-term monitoring programme in Erkuta from 2007 to 2019. Only data from plots that were surveyed during the whole period were used   (Fig. 1). Mass reindeer mortality due to rain-onsnow events occurred twice during the study period, in the winters of 2006and 2013-2014(Bartsch et al. 2010Forbes et al. 2016). The resulting unusually high abundance of reindeer carcasses during the spring of 2007 and 2014 was quantified as a binary variable, indicating whether mass mortality had occurred or not (Fig. 1). Because feeding of pups with reindeer carrion from the late winter has not been observed in our study area, only small rodent abundance was hypothesised to affect den selection. A summary of all variables considered in the analysis and their predicted relationships to den selection are shown in Table 2.  Table 2 Overview of variables considered in the analysis of breeding den selection by Arctic foxes (V. lagopus) in southern Yamal Peninsula (Russia) a Indices were scaled to a mean of 0 and a standard deviation of 1 over the different dens b Hare habitat quality index was correlated to the lemming habitat quality index (Pearson's r = 0.72) and was therefore not tested separately c The small rodent abundance in spring was scaled to a mean of 0 and a standard deviation of 1 over the different years

Statistical analysis
Breeding den selection was analysed by constructing models that tested three competing hypotheses describing a selection for (1) shelter quality, (2) proximity to prey resources, and (3) a selection for both shelter and prey resources. This analysis was conducted in two steps.
First we selected the variables that best described the preferred shelter, and the variables that best described the preferred prey resources. These variables were selected by comparing several candidate models with variables considered to reflect shelter quality and prey resources (Tables 2,  3). We used generalised linear models with a quasibinomial error distribution describing the number of successes (times a den was occupied) and failures (years visited minus times occupied). Unless indicated otherwise, all statistical analyses were carried out in R statistical software (R Core Team 2019). A quasibinomial distribution was used to account for overdispersion. The most suitable models were selected based on ΔQ-AICc values computed using the overdispersion parameter value (c-hat) of the most overdispersed model with R package "AICcmodavg" (Mazerolle 2019). Q-AICc values are the quasi-likelihood counterparts of AICc and are used instead of AIC values to correct for small sample sizes and a quasibinomial distribution. Models with ΔQ-AICc < 2 were considered equally supported by the data. However, as the goal of this step was to select a reduced number of variables for shelter and prey to be used in the subsequent analysis, we selected the model with the lowest ΔQ-AICc, unless a simpler model was within a ΔQ-AICc score of 2, in Table 3 Overview of the tested models describing breeding den selection of Arctic foxes (V. lagopus) in southern Yamal Peninsula, Russia Models are sorted by (Q-)AICc values and were constructed to test competing hypotheses describing a selection for shelter quality, prey resources, or both shelter and prey. Interaction terms were used to test whether the relative emphasis on selection for shelter quality or prey resources varied in years with different small rodent abundances a See Table 2 for variable descriptions and hypothesised relations to den selection b Using a c-hat value of 2.1 c Using a c-hat value of 1.9 d Den ID was used as a random factor in the analysis e Arctic fox breeding activity in the study area is shown to be related to yearly small rodent abundance and reindeer mortality in late winter ). These variables are therefore taken into account in all yearly models which case the simpler model was selected (Burnham and Anderson 2002).
In the second step, the selected variables were used to construct models describing a selection for shelter quality, prey resources, or both. Mass reindeer mortality and the yearly abundance of small rodents were included as variables in all models to take into account their effect on Arctic fox breeding activity in the study area . Additional models included interaction terms to test whether the relative emphasis on selection for shelter quality or prey resources varied in years with different small rodent abundances (Table 3). These generalised linear mixed effect models (GLMMs) constructed in the R package lme4 (Bates et al. 2015), used a binomial error distribution and had each yearly visit to a potentially occupied den as one data point. Den ID was used as a random factor in the GLMMs. These final models describing the competing hypotheses were compared using AICc values, where models within ΔAICc = 2 were considered equally supported by the data (Burnham and Anderson 2002).
Because several models received equal support from the data, model averaged effect sizes of explanatory variables in models within Δ AICc = 2 are reported as odds-ratios with 95% confidence intervals. Correlated variables (Pearson's r ≥ 0.7) were not used in the same models because of collinearity. Neither were correlated variables tested separately in similar models to reduce the number of tested models and the chances of finding spurious result (Table 2) (Burnham and Anderson 2002). Model fit was assessed graphically by plotting the residuals and assessing the importance of outliers. The residuals of all models were tested for spatial autocorrelation with the Moran's I autocorrelation coefficient in R package "ape" (Paradis et al. 2019). For the GLMMs, Moran's I was calculated for the estimates of the random intercepts. Temporal autocorrelation was assessed graphically with the Auto-and Cross-covariance and -correlation function estimation in R. No significant spatial autocorrelation (p > 0.49 for all models) or temporal autocorrelation was observed.

Habitat map
The classified habitat map showed that the most common habitat types are wet lowland, water bodies, and dry tundra, representing, respectively, 44%, 18%, and 15% of the study area (Fig. 2). The overall kappa index of agreement (Rossiter 2014) assessed by random points was 0.74. This has previously been considered to indicate "fair to good agreement beyond chance" (Fleiss et al. 2003). The lowest accuracy Fig. 2 Habitat map of the study area including den locations in southern Yamal Peninsula, Russia. The map was classified from a Landsat-8 image courtesy of the U.S. Geological Survey and the ArcticDEM dataset (Porter et al. 2018) was obtained for willow thickets (Table 4). Especially small willow thickets located on hillsides-which are known to be rich in prey species (Sokolova et al. 2014)-were not always recognised on the habitat map. The overall tau index of agreement (Rossiter 2014) assessed using set-aside training area pixels was 0.96. This indicates that the map more accurately recognises large areas with homogeneous cover. The confusion matrices of both methods are given in Online Resource 2.

Prey species per habitat type
Voles were most abundantly trapped in willow thicket edge (mainly M. gregalis) and wet lowland (mainly M. middendorffii). Lemmings (mainly D. torquatus) were most abundantly trapped in birch shrub and dry tundra, but they were considerably less abundant than voles (Fig. 3). Hare faeces were most abundantly counted in willow thicket edge and dry tundra, and ptarmigan faeces in willow thicket edge and birch shrub (Fig. 4).

Breeding den selection
The first step in the model selection for shelter quality resulted in two models within ΔQ-AICc = 2. Both included the number of entrances as explanatory variable, with the best model suggesting an additional positive effect of southern exposure (ΔQ-AICc = 0.97; Table 3). The model selection for prey resources also resulted in two models within ΔQ-AICc = 2. These models used the amount of river valley habitat and the lemming habitat quality index (ΔQ-AICc = 1.54; Table 3) as explanatory variables. Based on these results, we selected number of entrances and the amount of river valley habitat as explanatory variables to be used the final step of the analysis.
Using these variables to test our competing hypotheses resulted in three models that were all within ΔAICc = 2 (Table 3) and could thus be considered equally supported by the data. These models included the number of den Table 4 Users and producers accuracy of the habitat map of the study area, per habitat type. Accuracy was assessed using 97 random points 1 Fraction of pixels of a habitat type on the classified image that accurately represent this habitat type on the ground 2 Fraction of pixels of a habitat type on the ground, that appear as the correct habitat type on the classified image  . 3 The abundance indices of small rodents (voles and lemmings) trapped in five different habitat types in the study area in southern Yamal Peninsula, Russia. Error bars indicate 95% confidence intervals generated by a generalised linear mixed model with a negative binomial distribution and log-link. The models were constructed in SPSS (IBM 2013), and spatial units were taken into account as random factor. Small rodents in birch shrub were trapped using a different method and were therefore reduced by a factor 0.67 to enable comparison Fig. 4 The abundance of hare (L. timidus) and ptarmigan (L. lagopus) faeces sampled in five different habitat types in the study area in southern Yamal Peninsula, Russia. Error bars indicate 95% confidence intervals generated by a generalised linear mixed model with a negative binomial distribution and log-link. The models were constructed in SPSS (IBM 2013), and spatial units were taken into account as random factor entrances (shelter quality), the amount of river valley habitat (prey resources), and an interaction between the amount of river valley and the abundance of small rodents in spring as explanatory variables (Table 3). Number of entrances was included in all three best models, the amount of river valley was included in two of them, and the interaction term in one of them. The model averaged effect sizes of the variables included in the final models within ΔAICc = 2 indicates that the odds of den use were 1.07 times greater for each added entrance (CI 1.02-1.13, Table 5). Furthermore, the odds of den use were 1.23 times greater for a 1-km 2 increase in the amount of surrounding river valley habitat within 1.5 km of the den, yet the confidence interval of this effect size included the null value of 1 (CI 0.98-1.55, Table 5). The odds of den use increased during years with mass reindeer mortality, and during years with a higher abundance of small rodents, although the confidence interval for the effect of small rodent abundance included the null value of 1 (Table 5). The effect sizes also indicated a positive interaction between the yearly abundance of small rodents and a preference for dens surrounded by more river valley habitat, yet the confidence interval for this positive interaction also included the null value of 1 (Table 5).

Discussion
In line with our predictions, Arctic foxes in the low Arctic tundra of Erkuta preferred to breed in dens with more den entrances, our variable used as an indicator of shelter quality. The data also supported our prediction that Arctic foxes preferred dens surrounded by more prey resources, quantified by the amount of river valley habitat. The results therefore suggest that Arctic foxes in our study area place emphasis on both shelter and prey in den selection (hypothesis 3). That said, prey was only included in the top models when taking into account shelter quality. Furthermore, the confidence interval of the effect of prey on den selection included the null value of 1. The strength of evidence for the importance of prey resources in den selection in our study areas is therefore limited. The data most strongly supported the importance of shelter quality because it was included in all of the top models and the confidence interval of the effect size did not include the null value of 1.
The data did not support our prediction that Arctic foxes in our study area put relatively greater emphasis on selection for prey resources during years with a lower abundance of prey species. On the contrary, one of the top models included a positive interaction between the yearly abundance of small rodents and a preference for dens surrounded by more prey resources. Because the confidence interval of this interaction included the null value of 1, its strength of evidence is limited. Below we discuss each aspect of den selection separately.

Selection for shelter quality
Several structural den characteristics were considered to positively affect the quality of shelter. The number of den entrances best described the preferred shelter, and was positively correlated to a den's use for reproduction (Tables 3, 5, Fig. 5a). This is in accordance with the findings of previous studies, the usual interpretation being that larger dens with more entrances have more possible escape options, and therefore provide greater shelter against predators (Anthony 1996;Dalerum et al. 2002;Gallant et al. 2014). Possible predators of Arctic fox pups in Erkuta are rough-legged buzzard, white-tailed eagle, skuas (Sterconarius spp.), raven, and gyrfalcon (Falco rusticolus) Booms and Fuller 2003;Chevallier et al. 2016). Wolverine, red fox, domestic dogs, and people also pose a threat to adults (Frafjord et al. 1989;Tannerfeldt et al. 2002). One could alternatively argue that more frequently used dens experience more digging activity, and therefore acquire more entrances (Dalerum et al. 2002). Such a reverse cause-effect relationship between den use and the number of entrances is possible, especially because the sandy soil throughout our study area provides relatively easy digging conditions. Still, the energetic and time investments of digging extra den entrances is expected to yield fitness benefits, such as increased shelter against predators.
Previous studies have found that dens are often constructed on sites with less snow cover and better drainage (Smits et al. 1988;Prestrud 1992a;Tannerfeldt et al. 2003;Szor et al. 2008), and often face towards the south (Smits et al. 1988;Prestrud 1992a;Nielsen et al. 1994). These favourable microclimatic conditions may improve living conditions whilst the den is occupied (Gallant et al. 2014), but also improve digging conditions whilst the den is under construction. The number of den entrances may therefore

Selection for prey resources
Several species and habitat types were considered to be important prey resources to denning Arctic foxes in Erkuta. The amount of river valley habitat best described the preferred prey resource surrounding dens, and was included as explanatory variable for den selection in two of the top models within 2 ΔAICc (Table 3). Nevertheless, confidence interval of the effect size of river valley habitat on den selection included the null value of 1 (Table 5, Fig. 5b). Furthermore, river valley was only included in the top models within 2 ΔAICc when taking into account the number of den entrances. This suggests that prey resources are also of importance to denning Arctic foxes in Erkuta, but less conclusively than shelter quality. The diet of Arctic foxes in the study area consists mostly of small rodents and waterfowl, but there is a large variety of other potential prey species that are consumed; including hare and ptarmigan but also muskrat, waders and passerines, reindeer carcasses, and even fish . A possible interpretation of the results is therefore that Arctic foxes sustain themselves on a variety of prey resources and do not select for dens near any particular species' of prey. Nevertheless, Arctic foxes may prefer breeding in dens that are surrounded by comparatively productive habitat found in the river valley. The river valley in the study area is generally more productive than the upland. It is composed of numerous rivers, lakes, and shrubs of higher stature compared to the more exposed upland, and it hosts a higher density of breeding birds (Sokolov et al. 2012). Shrubs, especially willow thickets, provide important ecological functions in open landscapes such as shelter, food, and breeding sites for numerous insects, birds, and mammals, including voles (den Herder et al. 2004;Ripple and Beschta 2005;Henden et al. 2011b;Ehrich et al. 2012). The river valley can therefore be expected to host a richer and more diverse assembly of prey species to Arctic foxes compared to upland areas.
Previous studies that showed an effect of prey distribution on Arctic fox breeding den selection have mostly been conducted in the high Arctic. In Svalbard, dens near bird cliffs (a rich and stable source of prey) were more likely to be used for reproduction (Eide et al. 2012). On Bylot Island (Canada), Arctic fox dens in better lemming habitat and dens closer to a goose colony were more likely to be used (Szor et al. 2008). The bird cliffs and the goose colony described in these studies represent highly localised prey resources. In such high Arctic, less productive ecosystems there may be a greater dependence on a few prey-rich habitats, perhaps resulting in clearer selection pressures for specific prey species than in more productive, low Arctic areas.
Besides being generally more productive, the habitat structure of our study area also represents a patchy mosaic (Fig. 2) instead of being highly localised. Each den is therefore surrounded by more or less the same habitat types, differing only in their relative proportions. There may therefore be less variation in habitat quality between the dens in this study compared to previous studies, resulting in weaker selection pressures between dens. Furthermore, faeces counts failed to take into account seasonal differences in habitat use. For example, hares use the upland tundra in winter, when this more exposed habitat type is wind-blown and vegetation is within relatively easy reach, but much less during the Arctic fox breeding season. Decomposition rates of faeces in the drier uplands are also likely to be lower and detectability higher. These factors may have led The relation between the number of den entrances (averaged over the years) and the probability that the den was used for reproduction. b The relation between the amount of river valley habitat within 1.5 km of the den and the probability that the den was used for reproduction. The solid lines show the model aver-aged estimate for the generalised linear mixed effect models that best described den selection (within ΔAICc = 2). The dotted lines show the upper and lower 95% confidence intervals of the model averaged estimate. The jittered circles show single data points (den visits) that indicate whether a den was in use or not to inaccuracies in estimates of prey resources around dens during the breeding season.

Effect of yearly fluctuating prey on den selection
The yearly fluctuating abundance of small rodents was hypothesised to affect the relative risks of starvation and predation, and thus predicted to affect the relative emphasis placed on shelter quality or proximity to prey resources in den selection. Our results show no statistical support for this prediction. On the contrary, the results suggest an opposite trend, where Arctic foxes may prefer to den in the productive river valley during years with higher small rodent abundance (Table 5). A possible explanation is that the yearly small rodent abundance in the study area (Fig. 1) varied too little to affect the relative emphasis on selection for shelter or prey resources to a significant degree. Furthermore, Arctic foxes in our study area may not be dependent on small rodents to the same extent as other mainland populations and use several alternative prey resources . A high yearly variation in lemming abundance did affect den selection by Arctic foxes in the simpler ecosystem on Bylot Island (Giroux et al. 2012). By counting reindeer carcasses, Eide et al. (2012) were able to show that inland dens in Svalbard were used more frequently in years with more reindeer carcasses. These dens were located further away from the rich marine food-web. In our study, selection pressures may comparatively have been too weak to see yearly differences.
Although we found no significant temporal autocorrelation in den use between consecutive years, a certain degree of breeding site fidelity may still reduce the yearly flexibility Arctic foxes have to select a den. Marking of individual foxes and telemetry studies would provide a greater insight into the degree of breeding site fidelity, territories, and habitat use of Arctic foxes.

Future implications
The results of this study suggest that Arctic foxes on the low Arctic tundra of southern Yamal Peninsula select for dens that provide better shelter quality, but also take into account the den's proximity to prey resources. This finding is a valuable contribution to what we know from other low Arctic and Scandinavian study areas where large dens likely provide shelter against a diverse community of potential predators (Dalerum et al. 2002;Gallant et al. 2014). In the high Arctic, with little predators and fewer prey species, resource availability appears as the most important factor (Szor et al. 2008;Eide et al. 2012).
As tundra ecosystems become more productive due to climate change (Myers-Smith et al. 2011), and generalist predators find their way north (Sokolov et al. 2016), the relative risks of predation and starvation are expected to shift. This may spur further changes in the reproductive requirements and habitat selection of Arctic foxes. Indeed, Gallant et al. (2014) already showed that the presence of red foxes can affect den selection by Arctic foxes. This is echoed in a comparison between our findings from Erkuta, and observations from Varanger Peninsula in Norway where red foxes are more abundant; whilst Arctic foxes in Erkuta may prefer to den in the productive river valley, in Varanger it has been suggested that a large population of red foxes excluded Arctic foxes from the more productive landscapes at lower altitude (Killengreen et al. 2007;Ims et al. 2017). Thus, by comparing results from study areas that differ in environmental conditions, predators, and prey resources we gain a deeper insight into the adaptive abilities and fate of circumpolar Arctic fox populations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.