We re-evaluated the relationship between territorial Svalbard rock ptarmigan male presence and ecological relevant variables related to vegetation, terrain and snowmelt, building on the lessons learned from a previous regional habitat multiscale model, to predict breeding habitat suitability of this high-arctic, endemic ptarmigan on a large spatial scale. We used 11 years (2000–2010) of presence/absence data of territorial males, a multi-scale generalized linear modelling framework (glms) and recent advances in digital satellite based vegetation mapping. The final habitat model contained four significant predictors related to vegetation, terrain (elevation and slope) and a heat load index, as a proxy for snowmelt. Increasing amount of one particular habitat type, ‘established dense Dryas heath’ influenced habitat suitability positively at a small scale, while gentle sloping landscapes of intermediate steepness (10–25°) and elevation in the upper southwest facing part of the mountain slopes characterized occurrence at the landscape scale. The best model attained a relatively high explanatory power with a good ability to discriminate correctly between used and available sites. We extrapolated the habitat model to all parts of Spitsbergen with a current growing season sufficiently long for ptarmigan to complete a breeding cycle. The model indicates that only a small proportion of the vegetation covered land area (∼3.9%) is highly suitable for ptarmigan. In Svalbard the ptarmigan is an attractive game species, appears in low densities with low annual survival, has restricted availability of breeding habitat and is likely vulnerable to climate change. In such contexts, we suggest to use our predictive habitat model in harvest management of the species so that hunting quotas and efforts may be adjusted to habitat suitability, indicating the reproductive potential, in the hunting areas.
Understanding species spatial distribution across large geographic scales is important to conserve, manage and monitor species effectively. Extrapolating species distributions beyond areas of known presence, using predictive species distribution models, is a tool to gain information about continuous habitat suitability. This is particularly important in areas where surveys may be difficult owing to low abundance or rareness of the species or to remoteness of the study locations making surveys logistically challenging and costly (Jensen et al. 2008, Sharma et al. 2009, Speed et al. 2009, Revermann et al. 2012, Zohmann et al. 2013, Hacohen-Domene et al. 2015). Digital spatial data commonly covers large land areas that makes it possible to develop maps of habitat suitability and ecological forecasts (Lawler et al. 2011). Recently studies have pointed out the importance of multi-scale approaches to understand species distributions at various ecological spatial and temporal scales, especially when predicting habitat use in a changing climate (Storch 2002, Graf et al. 2005, Mayor et al. 2009, Revermann et al. 2012, Bean et al. 2014). Geographic information systems are readily available to a wide community of managers and conservationists. Multi-scale statistical modelling approaches enable researchers to identify environmental predictors and their spatial scales relevant to habitat use of the study species (Mayor et al. 2009, Bean et al. 2014). Thus, by combining knowledge on species habitat use and environmental variables relevant to spatial ecology, managers have the tools for making land use decisions and scientists more complex spatial models to explain ecological processes to predict future conditions or distributions and patterns in locations that have not yet been sampled (Drew et al. 2011).
The high-arctic archipelago, Svalbard, Norway, houses an endemic sub-species of the rock ptarmigan Lagopus muta, the Svalbard rock ptarmigan L. m. hyperborea. This non-cyclic ptarmigan species is among the northernmost ptarmigan in the world, vulnerable to climate change (Henden et al. 2017; see also Pernollet et al. 2015) and the most important small game species within the archipelago (Løvenskiold 1964, Pedersen et al. 2012, Soininen et al. 2016). The Svalbard rock ptarmigan persists at low densities (range 1.2 to 4.7 males km-2 in spring) with limited inter-annual population size variability (Pedersen et al. 2012, Soininen et al. 2016). Lack of breeding habitat is likely a limiting factor for the population size (Unander and Steen 1985). A recent study found surplus birds in spring (Pedersen et al. 2014b), but the size and possible trends of this population component is not known (Soininen et al. 2016). A better knowledge of the quality and distribution of ptarmigan breeding habitat would hence be very valuable both for managing the population in relation to hunting and for assessing likely effects of a changing climate on the population.
Habitat suitability models for territorial Svalbard rock ptarmigan males have previously been developed (Pedersen et al. 2007), but model predictions were limited to a section of central Spitsbergen (Nordenskiöld Land) due to lack of archipelago-wide digital vegetation data at that time. Pedersen et al. (2007) found increasing amount of vegetated xeric sloping areas at local scale (foraging; spatial resolution 200 × 200 m) and terrain variables (elevation at site and slope-aspect ruggedness index at spatial resolution 1000 × 1000 m) at a somewhat larger scale (> territory) to positively influence male presence in spring. Model predictions were made using both a detailed (only regional availability) and coarse DEM, but models performed equally well. Their study was, however, not able to address fine-scale variations in terrain exposure and sun radiation, both known to be highly influential for snow accumulation, distribution and melting patterns (DeBeer and Pomeroy 2009). Such characteristics are regarded important determinants of territory quality because early snowmelt allows for earlier onset of the breeding season (Unander and Steen 1985).
In this paper, we re-evaluated the relationship between Svalbard rock ptarmigan presence and environmental variables, building on the lessons learned in Pedersen et al. (2007), and expanded the habitat model predictions to cover the majority of the Svalbard archipelago. We used a longer time-series (11 years, 2000–2010) of population abundance monitoring data of territorial Svalbard rock ptarmigan males, and recent improved advances in digital satellite based vegetation mapping covering the entire archipelago (Johansen et al. 2009, 2012), as well as an updated digital elevation model. As Pedersen et al. (2007), we explored variables related to vegetation and terrain characteristics, relevant to ptarmigan breeding biology and fine scale habitat selection, using a multi-scale generalized linear modelling approach. Additionally, we also evaluated a heat load index as a proxy for snowmelt and spring onset within this framework. We discuss our results with focus on how this new and more comprehensive habitat model can guide future management of this low-density, endemic ptarmigan species.
Methods
Study area
Svalbard is a high-arctic Norwegian archipelago (74–81°N, 10–35°E; 62 700 km2), encompassing more than 500 islands, with the largest being Spitsbergen, Edgeøya, Barentsøya and Nordaustlandet. Deeply indented fjords characterize the areas along the west and north coast. Mountains, with summits reaching up to around 1700 m a.s.l., cover much of Spitsbergen. Svalbard spans three bioclimatic zones. The middle Arctic tundra zone covers central Spitsbergen, the northern Arctic tundra zone covers large areas along the west coast of Spitsbergen and the polar desert zone covers areas in the northern and eastern parts of the archipelago (Elvebakk 1999, 2005). Svalbard comprise 85% glaciers, barren and sparsely vegetated areas and only 15% of the land areas is vegetated. Of this about 22% is classified as rich vegetation (i.e. moss tundra, mires, fens, swamps and grassy heaths) and 78% as heaths and polar desert (Johansen et al. 2009). Svalbard has a relatively mild climate, shifting from typically humid, oceanic in the west to colder and drier climate conditions in the northern and eastern parts (Førland et al. 2012).
The study area for the annual population abundance monitoring of territorial Svalbard rock ptarmigan males was located on the largest island, Spitsbergen, in the northeastern part of Nordenskiöld Land. It encompassed approximately 1200 km2, and two large valleys, Adventdalen and Sassendalen (78°15′N, 17°20′E), surrounded by peaks reaching 1200 m a.s.l. (Fig. 1). Glacial valleys with wetland, ridge, and heath vegetation rarely above 5–20 cm in height dominates the study area (Elvebakk 1999).
Ptarmigan data
Ptarmigan monitoring data were available from annual point count surveys (2000–2010) of calling territorial Svalbard rock ptarmigan males. Pedersen et al. (2007, 2012) documented the monitoring design and survey protocol followed, and we only briefly summarize them here. Male ptarmigan establish territories in April, at which time continuous daylight prevails. The distinct territorial behaviour of ptarmigan males (body displays and territorial calls) can therefore be observed throughout the 24 h day (Unander and Steen 1985b). Point counts were conducted only under optimal weather conditions (calm, clear days) by trained observers for 15 min on each site 2–3 times per season mainly in the month of April (from 2 to 30 April and in two years up to 6 May). During each visit, presence or absence of male ptarmigan was recorded to determine the size of the pre-breeding population. During 11 years (2000–2010), 208 sites had been visited (range 1–29 visits; median = 22). For the present analysis, we consider only sites visited ≥10 times, resulting in a total sample size of 148 sites. Ptarmigans were detected at least once on 141 of these 148 sites, i.e. absence was only observed at seven sites (Fig. 1). To represent the landscape potentially available for ptarmigan males, we used a sample of 368 random sites (pseudo-absence sites), also used in Pedersen et al. (2007), and added the seven monitoring sites with absence of ptarmigan males. We selected the pseudo-absence sites from the same landscape (i.e. study area) as the presences count sites, and located them randomly (see e.g. Wisz and Guisan 2009 for recommendations) on ice free areas at least 500 m from both ptarmigan count sites and other random sites (Fig. 1). Thus, the pseudo-absences sites most likely represent similar ranges of values in the explanatory variables as the count sites. Due to the landscape characteristics of the valleys within the study area (e.g. narrow, steep river valleys), we were not able to randomly allocate such a large number of pseudo-absences sites as recommended by Barbet-Massin et al. (2012).
Figure 1.
Map of the study area for the annual population abundance monitoring of Svalbard rock ptarmigan (Lagopus muta hyperborea) males (smaller quadrat on Svalbard map). The black circular dots indicate point count monitoring sites (n = 148; 141 presence and seven absence sites; 2000–2010) and the smaller grey circular dots the pseudo-absence sites generated by Pedersen et al. (2007) (n = 368). The inset show the Svalbard archipelago with the survey region (black quadrat) and the extrapolation region (black line) for predictive habitat modelling.
Digital spatial data layers
We calculated and extracted all environmental variables using ArcMap with the Spatial Analyst extension (ver. 9.3; Environmental System Research Inst.). Pedersen et al. (2007) identified two significant scales for ptarmigan breeding habitat selection, a local scale around 200 m for vegetation variables, and a larger scale of 1000 m for terrain variables. Unless otherwise indicated we hence chose to evaluate vegetation variables at scales ranging from 100–500 m, and terrain variables ranging from 100–1200 m. Variables were calculated using neighbourhood estimates in a quadratic moving window with side length equal to the scale (Table 1, Supplementary material Appendix Table A2).
Vegetation
We used a digital vegetation map of the Svalbard archipelago that was based on unsupervised classification of 11 Landsat Thematic mapper scenes (TM/ETM+), from July and August during the time-period from 1987 to 2000 (spatial resolution of 30 × 30 m; Johansen et al. 2009, 2012 [see their Table 1 for details on the scenes]). This map was not available when Pedersen et al. (2007) developed the first regional habitat model. The digital vegetation map consisted of 37 vegetation types which were grouped into 19 habitat types (Johansen et al. 2009; Supplementary material Appendix Table A1) and later re-classified into 18 habitat types (Johansen et al. 2012). This re-classification did not affect the habitat types relevant to ptarmigan spatial distribution (Supplementary material Appendix Table A1). We selected nine habitat types, which were vegetated or contained important landscape elements such as boulder fields and rocky areas, known to be important for ptarmigan breeding (Unander and Steen 1985, Unander et al. 1985, Pedersen et al. 2007, 2014a, Wilson and Martin 2008, Zohmann and Wöss 2008, Revermann et al. 2012, Nelli et al. 2013). The proportion of habitat types and proportion of vegetated area was calculated for each of five spatial scales and values were extracted for both presence and pseudo-absence sites (Table 1, Supplementary material Appendix Table A2).
Terrain
We obtained terrain variables (altitude, aspect, slope and ‘Vector Ruggedness Measure’ [VRM]; Sappington et al. 2007) for both presence and pseudo-absence sites from the DEM of Svalbard which is based on aerial photographs (1:50 000) from 1990 and onwards. The spatial resolution of the DEM is 20 × 20 m and estimated uncertainty is 5–10 m (Norwegian Polar Institute 2010). Elevation was extracted as a local variable at each site. Slope and VRM were calculated as mean values for each of six spatial scales (Table 1, Supplementary material Appendix Table A2). Aspect was categorized in eight directional groups (north = 337.5–22.5; north east = 22.5–67.5 etc.; Table 1) and calculated at increasing spatial grid cell size as above.
Table 1.
Observed ranges of the environmental variables assessed for probability of presence of territorial Svalbard rock ptarmigan males In spring. Observed ranges are from presence and pseudo-absence sites at all spatial study scales. See Supplementary material Appendix 1 Table A2 for mean ± SD of presences and pseudo-absences sites for variables Included In the top ranked habitat model.
To describe sections within the valleys with high incoming radiation and hence likely early snowmelt, we calculated a ‘heat load index“ (termed HLI; Parker 1988). This index accounts for that the southwest side of a slope in northern land areas receives the most heat load during 24 h. HLI allowed for a distinction between the different parts of the characteristic valleys in the study area that contain landscape features, relevant to ptarmigan ecology, which is not accounted for by slope and/or aspect alone. HLI was based on the extracted slope and aspect values from the DEM by calculating cosine (aspect-225) × tan(slope) was aspect was expressed as degrees azimuth and slope angle is expressed in degrees (Parker 1988). We calculated HLI as mean values at increasing spatial grid cell size with the survey site in the middle of a quadrat (four spatial scales with sides from 200 to 800 m; Table 1, Supplementary material Appendix Table A2). We therefore categorised HLI into five classes capturing the upper north east hill slope, lower north east hill slope, valley bottom, lower south west hill slope, and upper south west hill slope.
Data analyses
Multi-scale habitat modelling approach
We followed the protocols of Johnson et al. (2006) and developed resource selection functions based on a useavailability design as described in Pedersen et al. (2007). We developed logistic generalized linear regression candidate model sets (glm's in R < www.r-project.org >) considering the probability of used sites (n = 141 with presence of males) versus available sites (n = 375; 368 with pseudoabsence and seven with absence based on point count data) as a response variable. We included environmental variables related to vegetation (proportion of habitat types and proportion of vegetated area at increasing spatial grid cell size), terrain variables (altitude, aspect, slope and VRM) and the ‘heath load index’ as predictors of male ptarmigan presence in a multi-spatial scale modelling framework to identify the most relevant ecological scale for ptarmigan breeding habitat. See previous section for details on predictor variables and spatial scales. To arrive at the initial list of candidate models, we combined vegetation variables at each of five given scales with terrain variables at each of seven different scales yielding a total of 35 candidate models. Each of these were again combined with HLI at each of four different scales, yielding a total of 140 candidate models. Based on this, we first explored importance of the habitat type predictors (n = 9 ecological relevant types; Supplementary material Appendix Table A1) at the smallest spatial scale (i.e. 100 × 100 m spatial scale; vegetation characteristics are predicted to be most important at ptarmigan fine scale habitat selection and foraging scales) in the candidate model sets by assessing ‘variable importance’ (VI). This statistic was calculated by summarizing the AICc weight values for the 40 best models (Anderson 2001). VI ranges from 0 to 1 with the value of 1 characterizing the most important variable. The habitat type ‘established dense Dryas heaths’ (Supplementary material Appendix Table A1) was the overall most important habitat type (VI = 1.00). All other habitat types had VI < 0.12, thus we decided only to include this habitat type in the subsequent development of candidate model sets. Second, we checked for co-linearity between the continuous predictor variables, as listed in Table 1, at each spatial scale using Pearson’s correlation coefficient. Mean slope and mean VRM were highly correlated at all spatial scales (maximum r was mostly above 0.70; see Dormann et al. 2013) and the variable ‘slope’ was retained for further development of candidate model sets for predictive habitat modelling. Additionally, we plotted the categorical variables, aspect and HLI, for visual inspection of correlative patterns. These variables appeared to slightly correlate, but the variable ‘aspect’ was not retained in any of the final model sets evaluated for predictive purposes and was subsequently ignored. Third, we tested all possible combinations of the retained vegetation and terrain predictors, relevant to ptarmigan breeding ecology, by using the dredge function in the R package MuMin (Barton 2009). We described non-linear relationships for the retained predictor variables, elevation and slope, using second order polynomials. Interaction terms (i.e. multiplicative terms for continuous variables, centred to mean 0) altitude × slope was only tested in the final candidate model set, however, there was no evidence for this interaction and we subsequently ignored it.
We selected the final habitat model using AICc and differences in delta AICc (Burnham and Anderson 2004). When difference in delta AICc was less than two units, we retained the simplest model. We also report AICc weights and residual deviance for the ranked models, and model coefficients for the selected habitat model. We checked the final model for spatial autocorrelation following Pedersen et al. (2007) using the package ncf and the function ‘correlog’ for R (Bjørnstad 2009), and calculated the degree of spatial autocorrelation in the model residuals at increasing spatial distance from zero to 30 000 m (increment 2000 m). The spatial autocorrelation coefficient decreased with increasing spatial scale and was always below 0.12 (range 0–1), thus we subsequently ignored spatial autocorrelation in the final predictive habitat model. We showed the effect of each predictor variable in the final model on the predicted probability of territorial male presence by marginal model plots (95% confidence intervals), based on predicted values and marginal residuals, for each level of the HLI. To do so we let the predictor variable take a set of values from the data while we held the other variables constant at an average values.
Habitat model evaluation
To assess discriminatory ability of the final habitat model, we calculated the threshold-independent ‘receiver operator characteristic’ (ROC) curve and the ‘associated area under the curve’ (AUC) (Pearce and Perrier 2000) using the presence—absence library for R (Freeman and Moisen 2008) and the ‘true skill statistic’ (TSS) (Allouche et al. 2006) using the R BIOMOD package. The AUC describes the discrimination capacity ranging from 0.5 (models with no discrimination ability) to 1.0 (models with perfect discrimination) (Pearce and Ferrier 2000). The TSS range from —1 to + 1 where the latter indicates perfect agreement and values of zero and less indicate a performance no better than random. Since we lacked an independent data set for ptarmigan male presence, we followed the protocols developed by Pedersen et al. (2007, 2014a) that internally cross-validated predictive accuracy (i.e. proportion of observations correctly classified in a random sample of data) using a 25 fold internal cross-validation implemented in the library DAAG for R (Maindonald and Braun 2012). Cross-validated estimates are presented as mean ± standard deviation of the iterative runs.
Goodness of fit of the final model was assessed by calculating Nagelkerke's R2 which quantifies the power of explanation of the model by comparing the null model to the fitted model (i.e. the smaller this ratio, the greater the improvement and the higher the R-squared) (Nagelkerke 1991). We also checked for over-dispersion in the final model (i.e. presence of greater dispersion in a data set than would be expected according to the statistical model in use) by calculating phi from the R library MASS (Venables and Ripley 2002). A phi-value (sum residuals [model, type = “pearson”]ˆ2)/ df.residual (model)) close to one indicates no over-dispersion. We found no lack of fit in our final model (phi = 0.96).
Predictive habitat suitability map
The Svalbard archipelago covers 7 degrees latitude and 20 degrees longitude, and consequently large gradients in growing season characteristics (e.g. onset of spring, length of the growing season; Karlsen et al. 2014). The possibility for ptarmigans to complete a breeding season is hence highly variable across the archipelago. Thus, models of breeding habitat suitability must be extrapolated with care into bioclimatic regions not represented in the survey data (Braunisch and Suchant 2010, Bain et al. 2015). We therefore chose to limit the extrapolation of the habitat model to areas where; 1) the growth season starts earlier than 7 July, and 2) the growth season is longer than 35 days (Karlsen et al. 2014, Karlsen et al. unpubl.). We based our choice on described mean dates for onset of incubation (range 15 June to 4 July; Steen and Unander 1985) and the fact that chicks forage only on bulbils of Bistorta vivipara that needs to be present at the time of hatching and during the first few weeks of rearing (Unander et al. 1985).
We calculated the habitat suitability raster map for territorial males (spatial resolution 100 × 100 m) using the model-building tool in ArcMap (ver. 9.2). Before the calculation we developed a land filter (total 8148 km2), based on Pedersen et al. (2007), consisting of land areas below 600 m (anticipated to be potential ptarmigan habitat) without glaciers and moraines (a buffer of 1 km was set around moraine areas). The map describes habitat suitability on a scale from 0 (not suitable) to 1 (highly suitable) and for visualization we categorized habitat suitability into four classes (marginal 0–0.099, low 0.10–0.399, medium 0.40–0.699, high 0.70–1.00). The flat valley bottoms and mountain plateaus (< 600 m) were natural borders between the lowest habitat suitability classes (classes marginal and low), and the remaining values from 0.1–1.0 were divided into three habitat suitability classes.
Results
The selected habitat model for probability of territorial male presence contained four significant environmental predictors variables related to vegetation at the smallest spatial scale (foraging) and terrain characteristics at landscape scales (> territory; see Unander and Steen (1985) for territory sizes). All other competing models had a ΔAICc larger than 2 (Table 2). The best model contained the habitat type ‘established dense Dryas heaths’ at the smallest spatial scale (100 × 100 m), elevation at ptarmigan survey site (m a.s.l.), slope (1000 × 1000 m) and the ‘heath load index’ (600 × 600 m) at landscape scales (Table 2, Supplementary material Appendix Table A2). Increasing amount of the habitat type ‘established dense Dryas heaths’ influenced habitat suitability positively (Fig. 2, Supplementary material Appendix Table A2). Male presences decreased with elevation and the highest occurrence was in gentle sloping landscapes of intermediate steepness (10–25°) in the upper southwest facing part of the mountain slopes that receive the most heat load (Fig. 2, Supplementary material Appendix Table A2). The probability of male presence was approximately 3.5 times higher in upper southwest facing slopes compared to the lowland valley bottom (odds ratio south-west upper hill slope / valley = 3.49, 95 CI [1.56, 7.80]; other odds ratios are not reported since they were not significant).
Table 2.
List of the five best candidate habitat models for probability of presence of territorial Svalbard rock ptarmigan males in spring according to differences in Akaike's information criterion corrected for small sample size (ΔAICc) and AICweights. The spatial scales of the variable in meter are in parentheses. Np = number of parameters estimated. The model in bold was selected for predictive habitat modelling.
The selected habitat model attained a Nagelkerke's R2 of 0.41, which indicates that the model has relatively high power of explanation. The AUC score of the model (0.84) indicates that the model has a good ability to discriminate correctly between the presence (used sites) and pseudo-absence (available sites) (i.e. 84 % of the sites were classified correctly) (Fig. 3). The TSS (0.52) of the best model also supported this result. The proportion of correctly classified observation (internal cross-validation for predictive accuracy) ranged from 0.7968 to 0.8183 (mean [± standard deviation] for 25 iterative runs = 0.8070 [± 0.0053]).
Figure 2.
Marginal plots for the fitted logistic regression model for probability of Svalbard rock ptarmigan male presence. The effect of each predictor variable on the predicted probability of male presence is shown by letting the predictor take a set of values from the data while the other variables are held constant on average values. Effect sizes are shown by shaded 95% confidence intervals for each level of the categorical ‘heat load index’ (HLI) variable. Valley = valley bottom; NE_up = upper northeast hill slope; NE_low = lower northeast hill slope bordering the valley bottom; SW_up = southwest upper hill slope; SW_low = lower northeast hill slope bordering the valley bottom.
Figure 3.
Summary statistics for model evaluation of the selected habitat model for predicting probability of presence of territorial Svalbard rock ptarmigan males. Left panel: observed proportion of sites with presences plotted against predicted proportion of sites with presence. Right panel: the ‘associated area under the curve’ (AUC) between 0.8 and 0.9 indicates good discrimination capacity (Pearce and Ferrier 2000). The absence and pseudo-absence sites were distributed as follows: marginal 0–0.099, presence (P) = 6% (note that this might happened by chance be due to our way of selecting habitat suitability classes) and absence (A) = 94%; low suitability 0.10–0.399, P = 24%, A = 76%; medium suitability 0.40–0.699, P = 51%, A = 49%; high suitability 0.70–1.00, P = 78%, A = 22%.
Extrapolating the best habitat model to most of Spitsbergen (Fig. 4) yielded only a small proportion of the vegetation covered land area (˜3.9%) to be highly suitable for ptarmigan breeding territories (Table 4).
Discussion
The predictive habitat model for probability of territorial Svalbard rock ptarmigan male presence at the scale of Spitsbergen had a high predictive ability (AUC = 0.84), and patterns of habitat use were consistent with the previous regional habitat model of Pedersen et al. (2007). Our refined habitat model was, however, able to better link presence of males to a particular habitat type, ‘established dense Dryas heaths’ at a smaller spatial scale related to forage plants, and to discriminate between terrain exposure which influence snow accumulation and snowmelt patterns important to territory establishment and breeding on the landscape scale (Fig. 2, Table 2). The best habitats were confined to narrow vegetated areas, in southwest facing slopes of intermediate elevation and steepness. Only a small proportion of the vegetation covered land area (˜3.9%) was of high-quality breeding habitat (Fig. 4, Table 3).
We were able to identify one particular habitat type at a smaller spatial scale (100 × 100 m) comparable to the previous regional habitat model by Pedersen et al. (2007) that identified increasing amount of vegetated xeric slopes at a larger scale (200 × 200 m) to be an important predictor of ptarmigan breeding habitat. Vegetation cover in form of ‘established dense Dryas heat’ (Fig. 2) increased the habitat suitability relative to the overall available habitat types across all sections of the landscape. While poor in terms of plant species diversity (Johansen et al. 2009, 2012), it contains high quality forage for the ptarmigan in the winter months (Unander and Steen 1985, Unander et al. 1985). When males return to the territory in late March and early April they ingest the small greenish leaves of Dryas octopetala along with plant parts of the wide spread Salix polaris and Saxifraga oppositifolia (Unander et al. 1985). These plants are accessible along the hill slopes and up towards the exposed ridge habitat with less snow accumulation (Elvebakk 1994). Our multi-scale modelling approach identified that the terrain attributes, elevation at count site and slope (1000 × 1000 m) which represented a landscape level spatial scale (> territory scale), were better predictors of ptarmigan presence than more local scale estimates (Fig. 2). This corresponds to findings in Pedersen et al (2007) and is reflected in habitat studies from elsewhere in the distribution range that demonstrate the importance of rocky, sloping areas in medium to high altitudes for ptarmigan spatial distribution (Zohmann and Wöss 2008, Schweiger et al. 2012, Pedersen et al. 2014a). In contrast to Pedersen et al. (2007), we selected a predictor based on slope rather than terrain ruggedness. Both variables were highly correlated. Terrain ruggedness is likely closer linked to terrain heterogeneity and probably associated vegetation diversity (Kudo 1991), which is important in rock ptarmigan habitat selection (Revermann et al. 2012, Schweiger et al. 2012), slope might better capture the topographic terrain characteristics predominant in our study area, and indeed in most of the larger glacial valleys on Spitsbergen. Here the landscape is dominated by e.g. large valley systems, with an elevational habitat gradient spanning from tundra mires in the valley bottom to exposed ridges in the upper part of the hill slope, interspersed by glaciers (Elvebakk 1994). The significant second-order polynomial term suggested a peak in habitat suitability at lower elevations in gentle sloping landscapes of intermediate steepness (10–25°; not present at flat homogenous areas nor in the steepest parts of the mountain slopes; Fig. 2). This complies with the studies of Unander and Steen (1985) and with rock ptarmigan habitat requirements in other distant parts (e.g. the Alps and northern Sweden), although with contrasting environmental conditions of the rock ptarmigan distribution ranges (Revermann et al. 2012, Schweiger et al. 2012, Nelli et al. 2013, Pedersen et al. 2014a). In Svalbard, Unander and Steen (1985) described the breeding territories to be on slopes, including ridges and one or more small gulleys, with at least one steep rocky area. These terrain characteristics combined with a good view of the surroundings, along with shelter, cover and good drainage makes the nesting grounds.
Figure 4.
Predicted habitat suitability (i.e. presence of territorial males in spring) of Svalbard rock ptarmigan on Spitsbergen, Svalbard, Norway. The black line show the delineation of the extrapolation region.
Table 3.
Parameter estimates and standard error (SE) for the predictor variables in the selected habitat model for probability of presence of territorial Svalbard rock ptarmigan males in spring . Reference level was set to the categorical class ‘Valley bottom’ of the heat load index. The estimates are differences (contrasts) between the intercept and the estimated effect. Estimates (SE) in bold are statistical significant (*p < 0.01; **p < 0.001; ***p = 0). See Fig. 2 for marginal plots of the fitted logistic regression model for probability of Svalbard rock ptarmigan male presence.
The inclusion of the ‘heat load index’ by Parker (1988), as a proxy of snow-melting patterns, enabled us to relate and refine habitat suitability to sections of the landscape with different levels of sun exposure (Fig. 2). The categorical HLI variable at this scale (600 × 600 m) likely captured the topographic characteristics of the wide glacial valleys (e.g. sloping mountain sides, ridges and v-shaped side river valleys and scars), and allows a relevant ecological distinction between the lower and upper parts of the mountain slopes and the valley bottom with different exposure. This is an improvement relative to the previous regional habitat model (Pedersen et al. 2007) where there was no difference in habitat suitability in landscapes with different terrain exposure. The increasing habitat suitability in upper southwest facing slopes is likely related to the importance of snow free patches and early thaw for territory establishment, as qualitatively described by Unander and Steen (1985) to favor territory establishment. They found a higher proportion of territories occupied by pairs and higher reproductive output in areas with early thaw as compared to areas with late thaw, and suggested that the quality of the male territory likely is the most important selective criterion for the hen and thus a good signal of the reproductive fitness of the male (Unander and Steen 1985). Our habitat model indeed supports their results and underlines the importance of snow cover (i.e. snow free patches) and spring onset. Such factors are major structuring determinants of habitat and forage availability and reproduction of many arctic species (Meltofte et al. 2008, van der Wal and Hessen 2009, Callaghan et al. 2011, Jensen et al. 2014). Ultimately, they may determine the spatial and temporal distribution of the Svalbard rock ptarmigan.
The selected habitat model had a high power of explanation of male rock ptarmigan presence in spring (Nagelkerke's R2 = 0.41), and the predictive ability of the model was good (˜8.5 out of 10 presence or absence sites were predicted correctly; Pearce and Ferrier 2000). The current model attained better predictive performance (AUC = 0.84, Fig. 3; TSS = 0.52) and predictive accuracy (mean [± standard deviation] for 25 iterative runs = 0.8070 [± 0.0053]) than the previous regional habitat model by Pedersen et al. (2007; see their Table 4), although methods for evaluation statistics differed. Despite the high predictive power, we limited model extrapolation to Spitsbergen due to the lack of independent validation data from other regions. Certain areas outside Spitsbergen have comparable climatic and environmental conditions to the survey area, in particular parts of Edgeøya, likely contain ptarmigan habitats. Indeed, considering the predicted changes in the growth season (Xu et al. 2013) and the close relations between plant biomass and summer temperature (van der Wal and Stien 2014), a larger proportion of lowland areas also on the eastern islands on Svalbard will experience similar climatic conditions to the ones found in our western study area today. Future collection of field data from regions outside the extrapolation region is required to validate the model and to predict across larger areas.
Table 4.
Summary of land cover statistics for probability of territorial Svalbard rock ptarmigan male presence In the extrapolation region on Spitsbergen, Svalbard, Norway. Habitat suitability classes are calculated from the predictive habitat model (Fig. 4).
Extrapolating the habitat model to Spitsbergen identified only a small proportion (˜3.9%) of the land area highly suitability as breeding habitat (Table 4), which is consistent with the previous regional habitat model (Pedersen et al. 2007). This implies that breeding habitat availability is limited also on a very large scale and may restrict the spatial distribution of this high-arctic ptarmigan sub-species. Several lines of evidence point to the Svalbard rock ptarmigan is limited by shortage of territories sufficiently attractive for breeding (Unander and Steen 1985, Pedersen et al. 2012, 2014b). This implies that habitat limitation, besides predation (Steen and Unander 1985, Prestrud 1992) and climatic variability (Steen and Unander 1985, Hansen et al. 2013), is an important factor in determining the size of the breeding population.
Our multi-scale habitat modelling approach identified relevant predictors (Fig. 2, Table 2) and their important ecological scales that ultimately determine the spatial distribution of suitable breeding habitats across a large high-arctic island (Fig. 4). Recent studies have demonstrated the use of predictive habitat suitability models for rock ptarmigan as tools to guide and inform science-based planning, management and monitoring (Booms et al. 2011, Pedersen et al. 2012, Revermann et al. 2012, Zohmann et al. 2013). For instance, Revermann et al. (2012) predicted future changes in habitat availability and consequences for rock ptarmigan distributions and demonstrated the importance of adapting conservation strategies to a dynamic environment. Zohmann et al. (2013) suggested combining habitat suitability with abundance data of rock ptarmigan from different points in time to evaluate changes in habitat suitability and species distribution over time. Booms et al. (2011) assessed changes in rock ptarmigan and gyrfalcon habitat suitability (fundamental niches) over a 200-year time-period and detected substantial environmental changes influencing long-term population viability and predator—prey interactions. Finally, Pedersen et al. (2012) stratified survey sites for long-term population abundance monitoring according to habitat suitability. These examples from different rock ptarmigan populations demonstrate the applicability and importance of having a basic understanding of spatial distribution at several spatial and temporal ecological scales to address conservation, management and monitoring questions.
The Svalbard rock ptarmigan is a highly attractive endemic game species where most of the offtake comes from the local breeding population (Soininen et al. 2016, Fuglei et al. 2017), and the current size of the surplus fraction of the population is unknown (Pedersen et al. 2014b, Soininen et al. 2016). It also appears in low densities (Pedersen et al. 2014b, Soininen et al. 2016) with low annual survival (Unander et al. 2015), has restricted availability of breeding habitat (Pedersen et al. 2007 and this article) and is likely vulnerable to climate change (Henden et al. 2017). In such contexts, we suggest our predictive habitat model, combined with hunting statistics (Soininen et al. 2016), to be used in harvest management of the species so that hunting quotas and efforts are adjusted to habitat suitability (i.e. indicative of reproductive potential) at the scale of a hunting area. A spatial harvesting approach, with protected core areas, could avoid potential accelerating of over harvest, which is commonly characteristic of quota-based harvest systems (McCullough 1996, Jonzen et al. 2001).
Acknowledgements
We thank especially Stein-Rune Karlsen for allowing access to partly unpublished digital data on the start and end of the growing season and scientific advice, Nigel G. Yoccoz for valuable statistical assistance, Anders Skoglund for GIS support, Virve Ravolainen for detailed information on vegetation types and Leif Binar Støvern for use of a photo.
Funding — The Svalbard Environmental Protection Fund, the Norwegian Polar Inst, and the Norwegian Institute for Nature Research provided funding for this study.
References
Notes
[1] Supplementary material (available online as Appendix wlb00239 at < www.wildlifebiology.org/appendix/wlb-00239>). Appendix 1.
