Data collection and study areas

Line-transect surveys were conducted in August from 1996 to 2011 in up to 60 survey areas across south-central and eastern Norway (Fig. 1). Four areas were surveyed from 1996, and new areas were subsequently added to the study design throughout the period. For practical reasons (weather, illness, shortage of voluntary field workers etc.), surveys were not conducted in all survey areas in all years and not all transects were sampled in a survey area every year. Because of the sub-alpine distribution of Willow Ptarmigan, survey areas were geographically clustered within five mountain regions (Fig. 1). Volunteer dog handlers with pointing dogs walked along predetermined transect lines and the free-running dogs searched the area on both sides of the line following the procedure of distance sampling (Buckland et al. 2001; Pedersen et al. 1999, 2004; Warren and Baines 2011). At each encounter, the number of birds (chicks, adult males, adult females and birds of unknown age/sex) and perpendicular distance from the transect line to the observed birds (m) were recorded. Pedersen et al. (2004) provide a detailed description of the sampling protocol. The number of years with data in each survey area varied between three and 15 (median = 7); the number of transects per survey area varied between two and 39 (median = 11), giving total transect lengths varying between 7.6 and 107 km (median = 33 km); and the number of encounters per year per survey area varied between four and 179 (median = 30). The lowest total transect length, number of transects and number of encounters were independent of each other (not from the same area and year). Data from different transects were pooled per site and year.

Fig. 1 Study areas (filled polygons) within mountain regions (open circles) in south-central Norway. RS Rondane, DF Dovre and Folldal, FH Forollhogna, GNE Glomma northeast, GSE Glomma southeast). Filled stars and filled triangles are the positions of meteorological stations and rodent trap sites, respectively Full size image

We defined recruitment as the number of juveniles per pair in a given survey area a given year, based on the encounters described above. To obtain estimates of recruitment, we estimated the proportion of juveniles (PJ) from the raw data. We used data from all survey areas and years, but included only transect lines with recorded encounters, and only encounters where the sex and age class (i.e. no observations with unknown sex or age, c.f. above) were noted, including pairs without broods. The total number of observations was 16,468 (per mountain region: DF = 2,321, FH = 6,672, GNE = 1,735, GSE = 2,192 and RS = 3,548, c.f. Fig. 1). To estimate the proportion of juveniles in each survey area each year, we used generalized mixed effect models with a logit link function for each mountain region separately (Crawley 2007), with number of juveniles/adult in each encounter as the dependent variable and a variable linking survey areas to year (called survey area-year) fitted as a random intercept. Then, we estimated the proportion of juveniles for each mountain region in each year by fitting a random intercept linking mountain region to year (called mountain region-year). This allowed us to estimate the proportion of juveniles from each encounter for each year in all survey areas and mountain regions separately. Large clusters are easier to detect than small clusters, and dogs spend more time searching close to the transect line than farther away (Pedersen et al. 2004). This might result in a size bias where average cluster size becomes larger at long distances compared to distances close to the transect line, and consequently, estimates of cluster size might be overestimated. As there is a positive correlation between cluster size and recruitment, we included distance from the transect line to the observation as a covariate in the models (Buckland et al. 2001). Consequently, we assumed that the effect of detection distance on cluster size was linear on the logit scale. While other relationships are also possible, low sample sizes in some areas/years would preclude more complex modelling of the relationship.

To estimate the recruitment of juveniles (number of juveniles/pair) we first estimated PJ, the proportion of juveniles in the sample estimated at the intercept (i.e. the back-transformed logit-value at the intercept). This corresponded to the proportion of juveniles at zero distance from the transect line, where detection probability is assumed to be one (Buckland et al. 2001). The number of juveniles/pair was then estimated as: \({\text{PJ}}/\left[ {\frac{{1 - {\text{PJ}}}}{2}} \right]\). The total number of estimates was 464 and 60 for the survey area and mountain region scales, respectively.

Weather and rodent data

The breeding season was divided into three time periods: Pre-incubation (PRE-INC), Incubation (INC) and Brood (BROOD). Based on an average hatch date of 24 June (Erikstad et al. 1985), we backdated 21 days of incubation (Westerskov 1956) and defined this period (3–24 June) as the incubation period (INC). The period prior to incubation was defined as the pre-incubation period (PRE-INC), and included laying and pre-laying days (1 May–2 June), and we defined the period after hatching (25 June–15 July) as the brood rearing period (BROOD). At the end of this period, most chicks are fledged and mortality is reduced compared to the preceding periods (Erikstad 1985a).

Local weather data

We obtained data on mean daily temperature (°C) and daily precipitation (mm) from local meteorological stations located >600 m above sea level. Not all stations recorded both temperature and precipitation, and many stations were opened or closed during our study period. Thus, we selected the five stations close to our survey areas with the most complete time series that included both temperature and precipitation data (Fig. 1). We measured distance between meteorological stations and the centre points of the survey areas and mountain regions. Ptarmigan data were then linked to data from the nearest meteorological station at both spatial scales. As a measure of temperature, we estimated the mean of all daily mean-temperatures (T) in all periods (T PRE-INC , T INC , and T BROOD ). Further we summed all daily precipitation in millimetres (RR) to obtain a measure of total precipitation in each period (RR PRE-INC , RR INC , and RR BROOD ). All local meteorological data were obtained from the open access database of the Norwegian Meteorological Institute at: http://www.eklima.met.no/.

Onset of plant growth

The OPG in spring is related to weather conditions such as snow-cover and temperature (Wielgolaski et al. 2011; Odland 2011). Variation in the timing of plant growth can possibly affect recruitment of juveniles through its effect on maternal nutrition and prey availability (Steen et al. 1988a; Moss and Watson 1984; Erikstad and Spidso 1982). To obtain estimates of OPG, we first used Geospatial Modelling Environment (Beyer 2012) to create minimum convex polygons (MCPs) for each mountain region, based on the centre points of survey areas within each region. Then, we extracted OPG from MODIS satellite data from 2000 to 2011 separately for each mountain region. The time-series of MODIS data have been atmospheric corrected and the measurement of OPG is well correlated with field observations of the onset of leafing (Karlsen et al. 2009, 2012). Because of data deficiencies, the OPG estimates for Rondane and Glomma southeast were only based on parts of the mountain regions. Nonetheless, we believe the data were adequate since the general year to year variation was present, and the focus of this study is the temporal, rather than spatial variability in driving factors. The mean start of the growing season across all years and regions was 3 June, while the earliest mean start of the growing season was 28 May in 2011, and the latest mean start was 10 June in 2005.

Large-scale climate variation

Large-scale climatic variability, such as the North Atlantic oscillation (NAO), is known to impact on population dynamics and ecological processes in birds (Forchhammer and Post 2000; Stenseth et al. 2002; Barnagaud et al. 2011). The NAO gives an index of the difference in atmospheric pressure over the North Atlantic and, during winter, it strongly influences temperature and precipitation in Northern Europe (Hurrell 1995). The focus in this paper is climatic variability during the breeding season (cf. 1 May–15 July, above), thus we choose to use a seasonal station-based NAO-index for the period May, June and July (NAO MJJ ) (Hurrell 2013) obtained from an open-access database at: https://climatedataguide.ucar.edu/guidance/ hurrell-north-atlantic-oscillation-nao-index-station-based.

Rodent abundance data

Steen et al. (1988b) demonstrated that recruitment of juveniles was strongly related to variation in rodent abundance. Abundance of rodents can function as an index of predation rates if the alternative prey hypothesis (Kjellander and Nordstrom 2003; Hagen 1952) is valid. We obtained long term rodent trap data from two sites in our study area; Åmotsdalen from 1991 to 2011 (Framstad 2012; Selas et al. 2011) and Fuggdalen from 1974 to 2009 (Selas et al. 2011) (see Fig. 1). Rodents were caught in snap-traps in September and abundances were indexed as number of rodents caught per 100 trap nights. The dynamics of rodent populations is complex, but one important determinant is the winter climate (Cornulier et al. 2013; Ims et al. 2008; Kausrud et al. 2008), where favourable conditions during winter can result in high densities in early spring and vice versa. There is often a close relationship between spring and autumn densities of rodents (Kausrud et al. 2008); hence, data collected in September are likely to provide a good index of rodent abundance throughout the Willow Ptarmigan breeding season. We linked Willow Ptarmigan data from survey areas and mountain regions to the nearest rodent trapping site.

Since the OPG data were restricted to the period 2000–2011, we used this period as the time-frame for further analyses. Then we omitted 36 estimates from the survey areas including two survey areas that were lacking data after 2000. For the mountain region scale we omitted six estimates. Hence, when assessing spatial synchrony in the period 2000–2011, the data consisted of 428 (57 survey areas) and 54 (five mountain regions) estimates of juveniles/pair at the survey area and mountain region scale, respectively. Further, as there were missing records in the meteorological and rodent data series as well (c.f. above), the dataset used for investigating climatic and predation effects was additionally reduced to 330 (57 survey areas) and 40 (five mountain regions) estimates of juveniles/pair with corresponding predictor variables at the survey area and mountain region scale, respectively.

Statistical analysis

We assessed spatial synchrony in recruitment rates by constructing matrices of pair-wise Pearson cross-correlations, both between survey areas and between mountain regions. Because of a lack of statistical independence of pair-wise cross-correlations, we calculated mean cross-correlation coefficients and confidence limits with a bootstrap procedure (Kvasnes et al. 2010). Pair-wise cross-correlation coefficients were then sampled with replacement to generate 100,000 matrices of randomly drawn correlation coefficients (Crawley 2007). This distribution was then used to estimate the mean, together with 2.5 and 97.5 % percentiles from the original matrix of pair-wise cross-correlation coefficients. We estimated bootstrapped means and percentiles across all survey areas, across survey areas within mountain regions and across mountain regions. We also assessed the level of synchrony in the rodent trap data by calculating a Pearson cross-correlation between the two trap sites.

The effect of climatic conditions and predation (indexed by rodent abundance) on recruitment of juveniles was modelled with linear mixed effect models at the survey area and mountain region scale (c.f. Fig. 1). We only considered additive effects and did not combine confounded variables. The local and regional climatic variables were modelled separately. At the survey area scale we included survey area, mountain region and year nested within mountain region as random effects, and at the mountain region scale we included mountain region and year as random effects. From the set of candidate models we used an information theoretic approach (Burnham and Anderson 2002) to select the most parsimonious model explaining the variation in recruitment of juveniles at survey area and mountain region scales, respectively. Because of the low sample size (~12 years), we used AICc as the selection criteria. ∆AICc values of <2 suggest that the models are equally parsimonious, but in such cases we selected the simplest model. As the amount of variance explained (R 2) by the explanatory variables can be of biological interest (Nakagawa and Schielzeth 2013), we estimated R 2 of the fixed effects from the most parsimonious models following the guides in Nakagawa and Schielzeth (2013).

To investigate how local conditions (weather variables and OPG) were related to the NAO index, we fitted linear mixed effects models with local variables as dependent variables and NAO MJJ as fixed effect. Since local weather variables and OPG data were derived from different locations (five meteorological stations and five mountain regions, respectively), we considered two models for each variable: one with additive effects of location and NAO MJJ , and one with the interaction between the two terms. All models were fitted with year as random effect. We used an information theoretic approach (as described above) to select the most parsimonious model (Burnham and Anderson 2002), and we calculated bootstrapped confidence intervals and used these to evaluate if the slopes relating NAO MJJ to the climate variable of interest from the selected models were different from zero.

All statistical analyses were carried out using R (R-Core-Team 2012). For the mixed effect models we used the lmer function in the lme4 package (Bates et al. 2011), and for the model selection procedure we used the MuMIn package (Barton 2013).