Is pre-breeding prospecting behaviour affected by snow cover in the irruptive snowy owl? A test using state-space modelling and environmental data annotated via Movebank

Tracking individual animals using satellite telemetry has improved our understanding of animal movements considerably. Nonetheless, thorough statistical treatment of Argos datasets is often jeopardized by their coarse temporal resolution. State-space modelling can circumvent some of the inherent limitations of Argos datasets, such as the limited temporal resolution of locations and the lack of information pertaining to the behavioural state of the tracked individuals at each location. We coupled state-space modelling with environmental characterisation of modelled locations on a 3-year Argos dataset of 9 breeding snowy owls to assess whether searching behaviour for breeding sites was affected by snow cover and depth in an arctic predator that shows a lack of breeding site fidelity. The state-space modelling approach allowed the discrimination of two behavioural states (searching and moving) during pre-breeding movements. Tracked snowy owls constantly switched from moving to searching behaviour during pre-breeding movements from mid-March to early June. Searching events were more likely where snow cover and depth was low. This suggests that snowy owls adapt their searching effort to environmental conditions encountered along their path. This modelling technique increases our understanding of movement ecology and behavioural decisions of individual animals both locally and globally according to environmental variables.


Background
Tracking individual animals using satellite telemetry has improved our understanding of their movement ecology considerably by providing databases that otherwise would be impossible to obtain. Indeed, the last several decades have nurtured a rise in the number of studies using the Argos system to analyse the movements of individual animals through time [1]. However, even with those ever-growing datasets of animal locations, we are still limited in our ability to address some broad questions pertaining to movement ecology of organisms. An important limiting factor for understanding the functional hierarchy of decisions underlying movements of individuals is our narrow set of analytical tools that can couple raw locations of individuals with any of the four basic mechanistic components of movement (i.e. motion and navigation capabilities as well as internal and external states), as described by Nathan et al. [2].
Indeed, thorough statistical treatment of Argos datasets and interpretation of specific animal behaviour during movement are hampered by three major limitations. First, those datasets have a coarse temporal resolution and often have irregular location estimates because most transmitters use programmed duty cycles and are dependent on the communication strength with moving satellites. This temporal irregularity increases the complexity of the data and prevents the use of common statistical analyses. Second, the locations themselves only come with a quality indicator, derived from the number of signals received, but do not provide clues as to what behaviour the individual animal was performing other than where it was. Third, the difficulty in linking any set of environmental conditions that are heterogeneous in source and format, and often obtained at different spatiotemporal scales than movement data, with any set of locations of an individual, often prevent users of answering rather basic questions as to what is affecting movements of the tracked organisms (but see [3,4] for examples).
State-space hierarchical switching modelling (hSSSM; [5]) can circumvent some of the inherent limitations of a typical Argos dataset, such as the restricted temporal resolution of locations and the lack of information pertaining to the behavioural state of the tracked individuals at each location. Indeed, this Bayesian approach allows the estimation of the most probable locations at fixed time steps based on the previous and forthcoming locations while taking into account the accuracy of each location provided by CLS Argos (see [5,6] for further details; CLS stands for Collecte Localisation Satellite). One can thus estimate locations at regular time steps and use them to perform subsequent movement analyses. Moreover, for each estimated location, the model can assess the probability of an individual being in one of two (or more, see [5]) predefined contrasting behavioural states according to proxies such as travel speed and turn angles, providing useful behavioural classifications (see [7] for a review). Such behavioural assessment can improve our understanding of underlying processes occurring at the individual level and, ultimately, of movement ecology. In addition, the recent development of the Environmental-Data Automated Track Annotation System of Movebank (Env-DATA; [8]) now allows automatic annotation of movement trajectories with instantaneous environmental conditions using large volumes of environmental data. This publicly available system eliminates the technical difficulties related to the annotation process such as data acquisition and interpolation, as described by Dodge et al. [8]. By coupling those two promising techniques (SSM and annotation of locations), we here show how movements of animals tracked with the Argos system can be influenced by environmental conditions encountered.
As one of the main predators of the tundra [9,10], the snowy owl (Bubo scandiacus) shows some of the most spectacular irruptive movements in boreal areas [11,12], concentrating in high density at some sites in given years and deserting them completely in subsequent years. Recent tracking studies were able to document breeding dispersal at the individual level and reported an almost complete lack of fidelity to breeding sites [11,12]. Those movements have been correlated with the abundance of lemmings (Lemmus and Dicrostonyx sp.; [13,14]) since snowy owls feed almost exclusively (>95%) on those small mammals during summer and rely on them for reproduction [10,12,15]. Those small rodents show highamplitude variations in abundance across years at a given site [16,17]. Such fluctuations can be spatially unpredictable [18,19], triggering long-distance movement of owls between breeding attempts to find sites harbouring high lemming densities [12]. While prospecting during the pre-breeding period (March -May), owls have been reported to reduce their speed and increase their turn angles in certain areas with more directional movements in between these areas [12]. The zigzagging movements shown by owls in some areas during prospecting are likely used to assess nesting conditions, and especially whether lemming density is high enough to entice owls to settle in that area for breeding. This prospecting occurs when the ground is still almost 100% snow-covered and the precise mechanisms used by owls to assess lemming densities remain unknown but likely include auditory and/ or visual cues [14].
We hypothesized that pre-breeding prospecting movements of snowy owls are affected by environmental conditions encountered during this period. For instance, a thick snow cover may prevent owls from hearing or seeing lemmings or their signs (e.g. tracks) or denuded sites at the margin of snow banks can facilitate lemming catching by owls. In addition, snowy owls are facing a time constraint during the pre-breeding season and individuals should benefit from finding a suitable breeding site as early as possible. Indeed, competition for territories is likely high as nesting sites with good food resources are probably limited. Moreover, an earlier onset of breeding is likely to increase reproductive success in this harsh environment characterized by a short summer, as observed in many other Arctic species [20,21]. Therefore, we predicted that snowy owls should spend less time searching in areas with a thick or extended snow cover, and we expected owls to settle where snow cover and depth are lower than average. Here we show how the combination of SSM with the Env-DATA System can help addressing this hypothesis.

Results
We tracked 9 of the 12 marked female snowy owls during the first year and 7 of them for an additional 2 years. Over the years, tracked snowy owls began their prebreeding prospecting movements on (± SD) 5 April ± 17 days and continued to do so for 36 ± 25 days. During that time, birds covered 1251 ± 1175 km. The birds often alternated between the moving and searching states, patrolling several distinct searching areas before eventually settling for the summer. Owls actually settled for summer on 12 May ± 21 days.
We used a total of 181 moving and 480 searching locations on land in the analyses ( Figure 1). Average snow cover encountered by the tracked birds varied significantly over time throughout the pre-breeding prospecting period in only one year (2009: F 1,261 = 25.7, P < 0.01) out of three (2008: F 1,233 = 0.02, P = 0.9; 2010: F 1,167 = 2.2, P = 0.14). The decrease in snow cover in that year (2009) was negligible over the prospecting period (−0.33% per day or 12% over the whole period). The overall average snow cover encountered by owls during all years was 78 ± 36%. Average snow depth however showed light but significant downward trends over time in all three years of tracking (2008: The probability of being in a searching state was not constant, as we found a significant negative relationship between searching probability and both snow cover (slope ± SE = −1.43 ± 0.37, variance explained = 20.5%) and depth (slope ± SE = −2.27 ± 0.60, variance explained = 20.5%; Table 1, Figure 2). Those covariates were however highly correlated (R = 0.65, P < 0.01, df = 659) and we therefore tested their effect on searching probability separately in 2 different models. Across the range of snow cover encountered by owls during their prospecting movements, the probability of being in searching state increased from 0.60 to 0.87 when birds passed from an area with 100% of snow cover to 0%. Similarly, the probability of being in searching state decreased from 0.86 to 0.41 when birds passed from an area with no snow to 0.80 m of snow. The average snow cover and depth (± SD) when owls settled at their breeding site was 69 ± 42% and 0.23 ± 0.18 m, respectively ( Figure 3).

Discussion
Our study showed that in snowy owls, there is a strong relation between prospecting behaviour (switching from moving to searching state) during pre-breeding movements and environmental conditions encountered and more specifically snow cover and depth. Indeed, birds tended to pass quickly and follow a direct path over areas with more complete or thicker snow cover and  In all models, individual's ID and year (nested in animal ID) were included as random factors (n = 661 locations).
were more prone to enter into searching behavioural state when encountering a thinner or reduced snow cover. Although the owls' diet during the pre-breeding period is not known precisely, lemmings are likely to be the dominant prey during that time, as it is the case during the breeding period, and snow cover can be a refuge for small mammals against avian predators. Indeed, hunting success was reduced by snow depth in wintering snowy owls feeding on small mammals [22]. A thick snow cover is therefore likely to reduce snowy owls' ability to detect signs of lemming presence as well as their capacity to actually catch them. In order to maximise resource acquisition early during the breeding season, owls may thus increase their searching behaviour in areas where their capacity to find and catch prey is high. Reproduction in snowy owls is tightly linked with lemmings [10] and those cyclic populations vary tremendously in abundance from one year to another at a given location. Their reproductive success is thus driven by their capacity to cope with the spatial unpredictability of food sources when finding a suitable nesting site. Indeed, most sites within the breeding range of snowy owls are actually vacated when lemming numbers are low (e.g. [9,12,14]) and only transient or non-breeding birds are  occasionally observed in those years ( [23], N. Lecomte unpubl. data; F. Doyle pers. comm.). These highly mobile predators are thus covering hundreds of kilometres during the short time window of pre-breeding movements to find a site harbouring a high density of lemmings [12]. Finding a good nesting territory as early as possible likely increases their reproductive success by ensuring access to a high abundance of food to raise their chicks and allowing young to reach the critical period of independence not too late during the summer, before the onset of the harsh winter conditions. Patrolling over a site with a thick snow cover in spring should hamper their ability or increase the time required to accurately assess lemming abundance. Therefore, it makes sense for owls to concentrate their searching behaviour in areas where this assessment is easier and quicker, such as areas with little snow, considering the potentially high cost of making a poor or a late decision in their settlement.
It could be argued that in snowy owls, the male is the sex selecting a nesting site and therefore should display such searching behaviour [14]. In our study, tracked females nonetheless exhibited extended prospecting behaviour every year, as described elsewhere [11,12]. This offers a new perspective on how pre-breeding strategies are expressed between sexes in this species. Females are likely building up their body reserves as well as looking for an optimal site for egg-laying, which can be more easily achieved where snow cover is thin or absent. If males adopt a similar strategy during the pre-breeding period and also concentrate their activity in those areas, this could also increase the chances of having a successful rendezvous when looking for a mate. Moreover, the possibility that the two sexes are travelling together during this period remains to be studied. Whether searching behaviour of males and females differ in terms of timing and intensity remains an open question and requires males equipped with transmitters during the pre-breeding season.
Even though owls concentrated their searching behaviour in areas with a relatively low and thin snow cover, they eventually settled in areas with moderate snow cover and depth (around 70% of cover and 0.25 m of depth) and not in areas with little snow. In winter, lemmings concentrate in areas of extensive and especially deep snow cover [24,25], a good winter habitat for them where their survival and perhaps even reproduction is high [26]. Therefore, owls may be facing a trade-off when they settle as areas with more snow have the potential to harbour a higher density of lemmings but the quality of these sites may be more difficult to assess and prey may be more difficult to access for a longer period of time, i.e. until snow-melt.
Snow cover and depth will be affected by the expected environmental changes occurring presently in the Arctic and over the next decades. Recent climate models predict reduction in snow cover and increased snow thinning in most arctic regions in spring, and such a reduction has already been documented in northern Canada [27]. A thinner and sparser snow cover may facilitate lemming detection by owls, but it might also reduce winter habitat quality for lemmings, thus reducing their overall availability. Other issues, such as increasing rain-on-snow events creating crusty snow that prevent herbivores from reaching their food under the snow [28,29] may also negatively affect accessibility to prey for predators like owls. Such events are likely to affect the prospecting behaviour of owls and they may need to increase the spatial scale of their searching movements.

Conclusions
The use of state-space hierarchical switching modelling as described by Jonsen et al. [5] allowed us to perform basic movement analyses on locations regularly spaced in time even though the raw Argos data were not. Moreover, this approach allowed us to incorporate estimates of uncertainty around raw Argos locations, while simultaneously providing us with parameter-based estimates of behaviour. With this approach, we were able to extract detailed information pertaining to the behavioural state of the tracked individuals. By doing so, we categorized critical behaviour (searching and moving) that related to key life cycle events (pre-breeding prospecting and summer settlement) of snowy owls. Behavioural discrimination using state-space hierarchical switching models has been used in the past in totally different species and environments (see [5,30,31]) for examples with seals, turtles and thrushes respectively). This versatile approach is however likely to be limited in the number of different behavioural states that it can successfully discriminate by the accuracy and temporal resolution of any dataset [3,30]. Finally, by combining this method with the annotation of environmental variables using the Env-DATA System of Movebank [8], we merged two modelling approaches to generate both detailed maps and relevant biological relationships. We here provided detailed evidence of how an avian predator of the tundra adapts its behaviour according to the conditions encountered during its pre-breeding prospecting movements. This coupling of techniques opens the possibility of using predicted variations in snow fall due to climate change to generate predictive maps of future sensitivity to snow attributes for a key arctic predator, early in the season. This approach thus increases our understanding of movement ecology and behavioural decisions of individual animals both locally and globally according to environmental variables.

Methods
Field activities took place on the southern portion of Bylot Island, Nunavut, Canada (73°N, 80°W). From 27 June to 11 July 2007, we captured 12 adult female snowy owls on their nests using a bow-net trap. We equipped owls with 30-g satellite transmitters (Microwave Telemetry, MD, USA; PTT-100) attached as a back-pack with a harness made of Teflon ribbon (see [32] for details). Tracking these birds for up to 3 years suggested no impact of the transmitter on their subsequent survival or reproduction [32]. During the first spring and summer (March to July 2008), transmitters were programmed to transmit for 5 h and then turn off for 49 h. Cycles of 4 h of transmission and 142 h off were programmed for the remaining battery life of the transmitters, including subsequent spring migration. Since we were interested in movement patterns on a broad temporal scale, those variable settings were programmed to increase transmitter life expectancy. We received locations of marked owls via the Argos system [6]. Each location was assigned by CLS Argos to a class (0, 1, 2, 3, A, B, or Z) according to its estimated accuracy, using the least squares filtering process. The accuracy of location classes 0, 1, 2, and 3 are > 1000 m, < 1000 m, < 350 m, and < 150 m, respectively [6]. There is no accuracy estimate associated with the remaining classes (A, B, and Z). We applied a state-space hierarchical switching model (hSSSM, [5]) to estimate daily locations of individual birds during pre-breeding prospecting movements for each year of monitoring (n = 3), using almost all locations (see below) provided by CLS Argos (see [12] for details). For each location, hSSSM assigns the probability that a bird was in a moving (probability close to 0) or searching (probability close to 1) behavioural state according to its speed and turning angle. hSSSM estimations were made using the R package bsam [33] under the R 2.15.2 environment (R Core Team) with the Markov Chain Monte Carlo (MCMC) sampler of JAGs [34]. Basic movement analyses were done using the trip package [35]. To ensure accurate and faster estimates with the hSSSM, we applied a speed filter (500 km/h, [36]) before running estimations to remove extreme error locations. The state-space hierarchical switching models were fitted to the dataset using 2 chains of 250 000 MCMC samples; the first 100 000 samples were discarded as a burn-in and the remaining 150 000 were thinned out to 3000 samples by retaining only every 50th sample to reduce autocorrelation. Estimations of parameters were based on these final 3000 samples. For each estimated parameter, convergence and absence of autocorrelation were checked, and we also applied the convergence diagnostic of Gelman and Rubin [37].
We defined pre-breeding prospecting movements as those occurring between the first time a bird switched to a searching state over potential breeding habitat after departing from their wintering grounds and the time of settlement for the summer. We defined individual departure date from wintering sites when behaviour state switched for the first time from searching to moving after 1 March (see [12] for details). A bird was deemed to have settled on a potential breeding site when it entered into continuous searching state and stayed in a 5km radius restricted area (i.e. home range) throughout July ( [38]; see also [12,33] for details on how summer settlement and home range were calculated).
Locations estimated with hSSSM were annotated with snow cover (%) and depth (m) at surface for the date of each location from the North American Regional Reanalysis of the National Center for Environmental Prediction (NCEP NARR). Original data were provided by the Physical Science Division of the US National Oceanic and Atmospheric Administration Earth System Research Laboratory. The datasets have a spatial and temporal granularity of 0.3 degrees and 20 minutes, respectively. We used the bilinear interpolation method through the Env-DATA System of Movebank [8].
We compared the snow cover and depth encountered during moving and searching behaviour states to examine behavioural responses to variations in environmental conditions experienced along the track [39]. As multiple dependencies were present in the dataset (measure of searching probabilities at different locations by the same individuals through several years), a general linear mixed models approach was used. For each model, residuals were inspected for linearity and homoscedasticity assumptions [40], and a transformation of the explanatory variable was applied when relevant. Models were fitted using R 3.1.1 (R core team 2014) and general linear mixed models (lme4 R package, [41]), with individual's ID and year (nested in animal ID) as random factors and the response variable (searching probability) was logittransformed. As nest sites and main prey (lemmings) occur exclusively on land, we limited our analyses to locations over land even though owls can travel over the sea during their annual movements [42]. We also evaluated snow cover and depth at the date that owls actually settled at a site to breed. The significance of models with either of the two variables of interest (snow cover and depth, this latter variable was transformed by squaring to reach model assumptions) was assessed using the information theoretic approach based on AIC [43] and comparison to a null model [40].