Reindeer green-wave sur ﬁ ng constrained by predators

. Migratory large herbivores in seasonal environments are known to follow the onset of new growth during spring, so-called green-wave sur ﬁ ng. This ensures prolonged access to forage with an optimal balance between forage quality and quantity. Many studies have focused on herbivores ’ ability to follow the spring ﬂ ush, but without considering potential constraints to sur ﬁ ng the green wave. The presence of predators is likely to be such a limitation, which could force herbivores to deviate from the optimal movement patterns in terms of forage access. We compared how well 319 reindeer ( Rangifer tarandus ) from seven different populations followed the green-up at different population densities of brown bear ( Ursus arctos ). We found that reindeer at higher bear densities selected movement paths with lower access to high-quality forage and deviated more in time from following the peak of the green wave, thus missing out on valuable forage. In addition, reindeer generally moved faster at higher bear densities, but this pattern was more consistent in forest habitats. Our results indicate that reindeer are forced to deviate from following the spring ﬂ ush and alter their movement pattern in areas with high bear densities, which may lead to reduced body condition for reindeer experiencing high predation risk. With the recent recolonization of large carnivores in northern ecosystems, it is critical to understand the direct and indirect effects of predators on large herbivores in order to assess effects on population dynamics and potentially cascading consequences on ecosystem function.


INTRODUCTION
Understanding animal foraging is important as it determines how populations are limited and distributed. The marginal value theorem predicts that the optimal patch residence time depends on the balance of depletion and renewal of resources across habitat patches (Charnov 1976). Similar simple principles are underlying most of habitat selection theory and methods, such as resource selection functions (Manly et al. 2002), and allow ranking of habitats depending on their profitability. However, the habitat profitability may change rapidly when seasonality leads to resource waves of abundant food that changes across spatial scales (Armstrong et al. 2016). For migratory herbivores, forage maturation is known to be an important driver of movement and habitat selection, predicting that they benefit from following the onset of fresh new growth along environmental gradients Sinclair 1988, Hebblewhite et al. 2008). Likewise, animals should respond to spatiotemporal variation in plant phenology between habitat patches within their seasonal ranges (Van Moorter et al. 2013).
The term "green-wave surfing" is now established to emphasize how animals are expected to follow resource waves with an optimal balance of forage quantity and quality (Merkle et al. 2016), and it fits within the broader context of energy landscapes (Wilson et al. 2012). How such waves are utilized by herbivores is less well studied (Armstrong et al. 2016), as it was previously logistically infeasible to measure both the resources and the animals at the scales that large herbivores operate on. The recent emergence of satellite-derived measures of plant green-up (the normalized difference vegetation index, NDVI) and animal tracking devices (global positioning systems, GPSs) has opened a new era, allowing quantification of the temporal and spatial distributions in forage quality across vast scales. Bischof et al. (2012) used NDVI time series to calculate the instantaneous rate of green-up (IRG) relative to speed of migratory deer using GPS, which was a major step forward (Fryxell and Avgar 2012). The development of this methodology has enabled detailed tracking of how well large herbivores follow the green wave during the growing season while moving from their winter to their summer ranges (Bischof et al. 2012, Merkle et al. 2016, Rivrud et al. 2016, Aikens et al. 2017. So far, studies of green-wave surfing have focused on how well migratory large herbivores track the peak of onset of new growth, without considering potential limitations caused by external disturbances. Animal foraging behavior is often the result of trade-offs between different demands, frequently leading to habitat selection processes varying across spatiotemporal scales (Senft et al. 1987). With the recolonization of large carnivores in North America and Europe (Chapron et al. 2014), it is important to understand the herbivore trade-off between safety and access to high-quality forage in northern boreal ecosystems (Atwood et al. 2009). Several herbivore movement studies have investigated predation risk-forage accessibility trade-offs. For example, in Isle Royale, USA, predation risk led female moose (Alces alces) with offspring to postpone the decision to migrate (Edwards 1983), parturient bighorn sheep (Ovis canadensis) ewes in Alberta, Canada, sacrificed forage quality for increased safety from predators when migrating to higher elevations (Festa-Bianchet 1988), and migratory caribou (Rangifer tarandus) in Newfoundland, Canada, selected habitats with abundant forage and reduced black bear (Ursus americanus) predation risk during calving .
The constantly growing literature, evaluating the influence of large predators on ungulate movement behavior and habitat selection, suggests a variety of antipredator behavioral mechanisms, depending on species, environmental conditions, and temporal and spatial scales of investigation. At landscape scale, ungulates may try to reduce predation risk by avoiding areas and habitats associated with higher risk of predation through migratory movements, or movements within their seasonal ranges Merrill 2007, Laundr e et al. 2010). The term landscape of fear is coined for such responses (Laundr e et al. 2014). Furthermore, ungulates may increase movements (Mitchell andLima 2002, Singh andEricsson 2014), move faster through risky habitats (Frair et al. 2005), or decrease movements (Panzacchi et al. 2009) to reduce predation risk. Antipredator behavior may also vary at finer scales, as an immediate response to predator presence (Valeix et al. 2009, Latombe et al. 2014, Basille et al. 2015 regulating vigilance behavior and flight responses (Lingle andPellis 2002, Creel et al. 2014), changing time allocated in different habitat types (Creel et al. 2005, Valeix et al. 2009), or change in movement speed (Fischhoff et al. 2007). In most cases, these behavioral changes lead to less time for feeding (Colman et al. 2003), higher energy expenditure (Parker et al. 1984), and feeding in less profitable habitats (Hern andez and Laundr e 2005, Hebblewhite and Merrill 2009). However, to our knowledge, there has been no evaluation of how ungulate green-wave surfing varies across a gradient of predation risk. During spring migration, risk perception may force ungulates to select less profitable paths. Moreover, migrating animals exposed to predation risk may alter movement and foraging patterns, reducing the ability of optimal foraging along the way, and females may reach the calving grounds too early, leading to a mismatch with vegetation growth (Post andForchhammer 2008, Sawyer et al. 2013).
We take advantage of a unique dataset of 319 GPS-marked reindeer (Rangifer tarandus) from seven herds across a population density gradient of brown bear (Ursus arctos) in Sweden. Recent work suggests that up to 30% of reindeer calves can be predated by brown bears during calving (Sivertsen 2017). The main predation period is concentrated to the first weeks following parturition (Sivertsen 2017), coinciding with the onset of spring green-up. Calving is an energetically demanding period for reindeer females, with high costs of gestation and lactation (McEwan andWhitehead 1972, Crête et al. 1993). Furthermore, reindeer females must utilize the short growing season in spring and summer to restore body reserves after depletion during winter (Parker et al. 2009). It is still unclear how reindeer females respond to brown bear predation risk at different spatial scales and whether behavior to reduce risk comes at a cost of optimal foraging during this critical period.
Our work unifies two largely independent conceptual frameworks: the energy landscapes (NDVI) and the landscape of fear (bear predation risk). It has been argued these frameworks collectively provide an enhanced mechanistic basis for understanding the spatial ecology and decision-making of animals (Gallagher et al. 2017), but empirical cases are few. Here, we hypothesize that different predation regimes, that is, different bear densities, affect reindeer's ability to maximize energy gain during spring and summer. Brown bear kill rates on reindeer calves vary with reproductive status and age, and increase with the time the bear spends within the calving ground (Sivertsen 2017). Thus, brown bear predation rates on reindeer neonates are to a large degree a function of number of bears on the calving ground, and we therefore assume that bear density is a representative index of predation risk. We use two metrics to quantify how well the reindeer are surfing the green wave: the cumulative IRG (CIRG) and the number of days from peak IRG. The IRG curve consists of two parameters: a scale parameter, which measures the number of days it takes to reach peak green-up (speed of spring) and the day green-up peaks (see Appendix S1: Fig. S1 and Aikens et al. 2017). For a given number of days away from the peak IRG, different pixels can have different IRG values depending on the speed of green-up; that is, if green-up is slow, the IRG value will be larger than when green-up is fast. Thus, the CIRG measures the total amount of the green-up that each individual experiences over a given period (Bischof et al. 2012), and is a combined index of access to high-quality forage and animal behavior. The number of days from peak IRG shows how well the individual surfs the green wave, and is thus a measure of behavior. Using these metrics, we predict (P 1 ) that reindeer have lower access to high-quality forage (lower CIRG) when bear density is high, as they must prioritize predator avoidance, and (P 2 ) that higher bear densities may cause reindeer to deviate more from surfing the green wave (increased number of days from peak), due to altered movement paths (moving to higher elevations and areas with snow cover) or due to changes in movement speed or stopover time in foraging habitats, to reduce predation risk. Furthermore, we predict that reindeer move (P 3 ) faster at higher bear densities due to flight responses from bear encounters or as an antipredator strategy to be less predictable in space in response to an actively searching predator (Lima andDill 1990, Fischhoff et al. 2007). Although brown bears are generalist feeders, previous studies have suggested that grizzly bears may actively locate areas above the tree line where caribou aggregate to calve during the short calving period (Reynolds and Garner 1987). The higher calf densities in semi-domesticated herds compared to forest-living wild herds may promote active hunting by brown bears on neonates (Sivertsen et al. 2016) and thus favor speed as a tactic to reduce predictability. Finally, as predation pressure on reindeer calves has been shown to be strongest during three weeks during and after peak calving (Sivertsen 2017), we predicted (P 4 ) more pronounced responses to bear density during the calving season, compared to the post-calving season.

Study area
S ami reindeer husbandry in Sweden is a pastoral system, covering 55% of the Swedish land base (Sandstr€ om 2015). In this study, locational data from reindeer were gathered from four forest herding districts (G€ allivare, Mal a, Udtja, and € Ostra Kikkejaure) and three mountain herding districts (Hand€ olsdalen, Njaarke, and Sirges; Fig. 1; Appendix S1: Table S1). The forest district calving and post-calving ranges are characterized by undulating boreal forests interspersed with mires and lakes. The mountain district calving ranges are all situated in the mountain region and mainly above the tree line. Reindeer density within the summer ranges (all-year land) during the years studied ranged from 1.3 to 3.4 reindeer/km 2 , with the highest densities in the mountain districts (Appendix S1: Table S1). In April, the reindeer herders initiate spring migration to the calving ranges, which takes place by foot or with trucks if migration routes are unavailable. To promote a safe travel for both themselves and the reindeer, herders strive to start migration before snowmelt. Migration by foot may take 1-3 weeks, depending on distance, terrain, and weather. Generally, reindeer moved by trucks will arrive earlier to the mountain area. Irrespective of whether reindeer arrive by foot or truck, they are left in lowland forest areas until the progression of spring has come far enough for the reindeer to move on their own to the calving grounds. Except for the gatherings for calf marking in the summer, the reindeer roam freely within the borders of the calving grounds to autumn ranges until the snow arrives and autumn migration is initiated.

Reindeer data
In total, 555,740 positions were gathered from 319 GPS-collared reindeer females, representing 409 individual reindeer-years over nine years (2003 and 2008-2015; Appendix S1: Table S1). Outliers in the reindeer GPS data were removed both by manual screening of the data and following Bjørneraas et al. (2010). We removed all data when the animals were moved in trucks. In total, 0.03% of the data was removed as erroneous/ truck positions. The fix success rate was 98.5% (range: 94.5-100% in the different herding districts).
Our focus period was from 11 May to 31 August and represents the growing season. Although the mean start day of spring growth (based on NDVI) in the area was 29 April (standard deviation [SD] = 18 d), these dates were chosen based on when we had continuous data for all seven herding districts. Individuals with 95% coverage in the focus period, and with a maximum gap in the logging sequence of 4 d, were retained for further analyses. The focus period was further split into two sub-periods, based on the risk of calf predation by bears. Calving in reindeer is concentrated to a few weeks, with a peak in mid-May (Espmark 1971). The first subperiod (11 May-9 June; hereafter termed calving) represents the reindeer calving season where calves are assumed to be under higher predation risk (Sivertsen 2017), than in the second sub-period (hereafter termed post-calving), representing the post-calving and summer season (10 June-31 August). We calculated the movement speed (in km/h) between two successive reindeer GPS locations as the Euclidean distance divided by the time between the locations. We discarded all locations where the individuals had moved at unlikely distances or speeds (>40 km/h or >10 km between fixes; 0.11% of the data). The mean daily movement speed per individual was calculated from all individual speed entries for each Julian day.

Plant phenology
Plant phenology was quantified using the satellite-derived normalized difference vegetation index (NDVI; Pettorelli et al. 2005). The NDVI has been shown to function well as a proxy for ungulate forage quantity (Hebblewhite et al. 2008, Hamel et al. 2009, Garroutte et al. 2016, although a short mismatch between understory and forest cover green-up can be expected. However, as we work with the rate of green-up and not the NDVI values directly, this is unlikely to bias the results in this study (Appendix S1: Fig. S1). MODIS TERRA MOD13Q1 satellite images covering the study area were downloaded from the NASA Land Processes Distributed Active Archive Center (LP DAAC 2000; https://lpdaac.usgs.gov/data_acce ss/daac2disk). The images were taken every ❖ www.esajournals.org 16 d, with a resolution of 250 m. Details on processing steps including removal of unrealistic values and noise reduction can be found elsewhere (Hird and McDermid 2009, Bischof et al. 2012, Rivrud et al. 2016. A double logistic curve was fitted to each pixel's annual NDVI time series (Beck et al. 2006, Hird and McDermid 2009, Rivrud et al. 2016. From this, the instantaneous rate of green-up (IRG) can be extracted by taking the first derivative of the part of the logistic curve that covers spring (Bischof et al. 2012). In early growth stages, grasses have the highest nutritional value (Van Soest 1994) and are also more easily digestible (Langvatn and Hanley 1993). Thus, the IRG gives an index of when forage is both abundant and of high quality. For each individual, we calculated the mean daily IRG experienced based on the corresponding GPS locations. From this, we calculated the CIRG for each individual in each of the two sub-periods as a measure of the total amount of high-quality forage experienced. Further details on IRG estimation can be found in Bischof et al. (2012). To evaluate how close in time the individuals were to the peak of the green wave, that is, how well they were surfing, we extracted the day of peak IRG for each pixel visited by the reindeer, calculated the difference in days from when they were present to the day of peak IRG for that pixel (Aikens et al. 2017), and averaged this within each subperiod for each individual.
We calculated greenscape metrics following Aikens et al. (2017), to account for underlying green-up patterns that may influence reindeer movement. For each individual, we simplified their movement trajectories within each sub-period as follows: (1) using a correlated random walk (Johnson et al. 2008) to interpolate points between GPS locations, (2) removing tortuosity to avoid oversampling in heavily used areas by applying the Visvalingam line simplification algorithm to the correlated random walks, and (3) sampling locations at 1-km intervals along the simplified movement paths. For each of the sampled 1-km locations, we extracted the day of peak green-up and the spring scale, where the latter is a measure of the time it takes to reach the peak green-up, that is, the speed of spring. We used the reciprocal of the spring scale to avoid issues with skewed distributions. The greenscape metrics calculated for use in the analyses were the rate, order, and duration of greenup. Rate of green-up was calculated as the average spring scale for all 1-km locations per individual each sub-period. Green-up order was calculated as the Spearman's rank correlation between the order in which the individuals visited the 1-km interval locations and the consecutive order in which these locations reached peak green-up. Finally, we extracted the mean Julian day of peak green-up for a buffer zone of 500 m around the first and last location in each simplified movement trajectory, and the duration of green-up was calculated as the difference in days of peak green-up between these two buffer zones. Further details can be found in Aikens et al. (2017).

Home range estimation and landscape features
Home ranges corresponding to the two subperiods (11 May-9 June and 10 June-31 August) were estimated for each individual by calculating 95% adaptive local convex hull (a-LoCoH) polygons using the "adehabitatHR" package in R (Calenge 2006, R Core Team 2016. LoCoHs allow for linear movement patterns, which are common when covering animals in constant motion, such as reindeer. When the a-LoCoHs were calculated, the parameter a was chosen according to Getz et al. (2007), stating that an avalue exceeding the two maximum distances between individual locations should always give the 100% density isopleth. In the case of failure, we increased the a-value using the sum of the three longest distances, up until the sum of the five longest distances to ensure success. The mean number of locations used per individual home range was 400 (range: 28-1439) during calving and 960 (range: 83-2983) during post-calving.
To control for landscape features that are known to influence plant growth (Mysterud et al. 2001 and forage access, we extracted maps of terrain ruggedness, slope, and aspect derived from a digital elevation model (DEM) covering Sweden (www.lantmateriet.se). Terrain ruggedness maps were made in R ("raster" package; Hijmans and van Etten 2015), while maps of slope and aspect were made using ArcMap 10.3.1 (ESRI, Redlands, California, USA). Aspect was converted to northness (cosine transformed) ranging from À1 (south facing) to 1 (north facing). Elevation was extracted directly from the DEM. Maps of minimum distances from power lines, railways, and large (>5 m wide) and small roads (<5 m) were made in ArcMap 10.3.1. All maps were rasterized with a resolution of 100 m. The individual a-LoCoH home range polygons for the two sub-periods were spatially joined with the different maps, and the corresponding means of each landscape feature were extracted using R. For the daily scale analysis of movement speed, we extracted the same covariates for each GPS location and averaged them by day.

Bear density index
The bear density was estimated from the latest scat survey (2006-2010; non-invasive DNA) conducted in each County ( Fig. 1; Appendix S1: Table S1), with data from the Scandinavian database of large carnivore surveys (www.rovbase. no; Bellemain et al. 2005, Kindberg et al. 2011. We used all bear scats where the individual bear had been identified and calculated density of scats with the density tool in ArcGIS Desktop (ESRI) giving a raster map with a 1000 m resolution. We calculated the mean of the bear density index within individual reindeer 95% a-LoCoH home ranges within each sub-period and for each location averaged by day for the movement speed analyses.

Statistical analyses
To investigate whether any spurious correlations related to bear density and movement patterns or green-up patterns were present, we estimated correlations between bear density and spring green-up patterns within herding districts, and between bear density and individual movement trajectories. We also added the movement trajectories of a random sample of individuals to a map showing the bear density gradient, as well as the green-up patterns within herding districts to visualize direction of movement.
We ran three models with different response variables to answer questions about movement in relation to bear density: access to high-quality forage within each sub-period (CIRG), average difference in days (absolute value) from the day of peak IRG within each sub-period (days from peak), and mean daily movement speed. For CIRG and days from peak, we ran linear mixedeffects models using "lme4" fitted with REML (Bates et al. 2015) in R. A generalized additive model (GAM) was run using "mgcv" (Wood 2011) in R when modeling mean daily movement speed.
Candidate predictors for CIRG and days from peak models were bear density index, sub-period (categorical: 1 or 2, representing calving and post-calving, respectively), rate, duration and order of green-up, elevation (m a.s.l.), terrain ruggedness index, slope (°), northness, reindeer herding district habitat type (mountain or forest), and the minimum distances to power lines, railways, and large and small roads (all in m), as well as the interaction between sub-period and bear density index. The same predictors were used as candidates for the model of mean daily movement speed, except sub-period, but with the addition of the non-linear Julian day (smooth term) and the interaction between the non-linear Julian day and bear density index, and between bear density index and habitat. In addition, for the speed model we included a fixed effect for GPS schedule (categorical), as the logging schedule varied between herds. To correct for autocorrelation in the mean daily movement speed model, we added lag variable with lag = 1 of the residuals as a fixed effect.
The residuals of the models were inspected, and to ensure model assumptions were met, the rate of green-up was log-transformed as it caused patterns in the residuals. We checked correlations of all pairs of candidate predictors before running models, and only one of the variables in each pair was retained when Pearson's r > |0.4|. After correlations were determined, we checked the remaining candidate variables for non-linearities using GAM plots, and models were checked for multicollinearity using the variance inflation factor. The need for random factors was determined using the Akaike Information Criterion. Models with year, individual reindeer id, reindeer herding district, or any combinations of these as random intercepts were compared. A total of 815 observations from 319 individuals were used in the analysis of CIRG and days from peak. The same individuals were available for the analyses of mean daily movement speed with N obs = 45,332.

RESULTS
Maps of individual movement trajectories in relation to bear density gradient and green-up gradient within herding districts can be found in the Appendix S1: Figs. S2-S3. The maps showed no indication of spurious correlations between bear density and green-up patterns or reindeer movement. Explicit tests of correlations showed no patterns between green-up parameters and bear density for pixels visited by the reindeer; r = À0.013 for start of spring and bear density, r = 0.089 for end of spring and bear density, r = 0.075 for speed of spring and bear density, and r = 0.036 for peak green-up and bear density. The correlation between the reindeer movement trajectory (Julian day) and bear density was r = À0.017. For all pixels within the herding districts, the correlation between start of spring and bear density was r = À0.043.

Access to high-quality forage
The final model explaining access to high-quality forage (CIRG) included rate (log-transformed) and duration of green-up, northness, bear density index, reindeer herding district habitat, subperiod, and the interaction between bear density and sub-period (Table 1a). Year and herding district were fitted as random intercepts. The proportion of variance explained by the fixed effects (marginal R 2 ) was 0.39, and the total variance explained by both fixed and random effects (conditional R 2 ) was 0.49. As predicted (P 1 ), increased bear density led reindeer to select movement paths with lower access to high-quality forage, but contrary to our predictions (P 4 ), there was no significantly different effect of bear density between sub-periods (Table 1a, Fig. 2a); the predicted CIRG decreased by 7.0 points from the lowest to the highest bear density during calving, and with 5.1 points from lowest to highest bear density during post-calving. Overall, the access to high-quality forage was higher during calving (Fig. 2a), and reindeer had increased access to high-quality forage when the spring green-up rate was slower (Table 1a). The model showed no significant effect of green-up duration, northness, or reindeer herding district habitat on CIRG (Table 1a).

Days from peak green-up
The final model explaining days from peak IRG included rate (log-transformed) and duration of green-up, northness, bear density index, reindeer herding district habitat, sub-period, and the interaction between bear density and subperiod (Table 1b). Year and herding district were fitted as random intercepts. The marginal R 2 was 0.82, and the conditional R 2 was 0.85. As predicted (P 2 ), the reindeer deviated more from following the green wave; that is, number of days from peak increased with increasing bear density (Table 1b; Fig. 2b), but contrary to predicted (P 4 ), there was no difference between sub-periods. Reindeer were predicted on average 10.5 d closer to the peak green-up at the lowest bear density compared to the highest. Reindeer were closer to the day of peak IRG in the calving period than in the post-calving period (Fig 2b), and in mountain herding districts they were closer to the peak green-up compared to forested districts (5.5 d; Table 1b). There were no significant effects of rate and duration of green-up, or northness (Table 1b).

Mean daily movement speed
The final model explaining mean daily reindeer movement speed (log-transformed) included bear density index, reindeer herding district habitat, rate (log-transformed) and duration of green-up, terrain ruggedness, northness and GPS collar schedule, the interaction between bear density and herding district habitat, as well as the non-linear Julian day main effect and in interaction with bear density index ( Table 2). The adjusted R 2 was 0.56.
As predicted (P 3 ), reindeer moved faster at higher bear densities, and the effect depended on the non-linear effect of Julian day; that is, there were periods when reindeer moved faster at higher bear densities, but also other days with no difference in movement speed depending on bear density (Table 2, Fig. 3). However, the effect was habitat-specific and more consistent in forest herding districts than in mountain herding districts (Table 2, Fig. 3). Overall, reindeer were predicted to move 4.6-8.3% faster at the highest bear density compared to the lowest during calving in forested district, while they moved 5.2-15.3% slower in mountain herding districts. During post-calving, they moved 3.9-12.6% faster at the highest bear density during postcalving, while in mountain herding district they could move from 17.1% slower to 7.1% faster at the highest bear density. The greenscape also influenced movement speed, with faster greenup rate leading to increased movement speed, while reindeer moved slower when duration of green-up was longer (Table 2). For landscape features, decreased terrain ruggedness led to faster daily movement speed, but there was no effect of northness (Table 2).

DISCUSSION
Several studies have documented that herbivores follow the green wave during migration, increasing their access to high-quality forage (Hebblewhite et al. 2008, Bischof et al. 2012, Merkle et al. 2016, Aikens et al. 2017). Here, we show how predation risk may limit the herbivores' ability to follow this spring flush of nutritious forage, causing a trade-off between access to forage and avoiding predation. We found that reindeer experienced lower access to high-quality forage when bear density was high (supporting P 1 ), but similarly during calving and post-calving (not supporting P 4 ). More specifically, the reindeer were forced to deviate from the green wave by being further in time from the peak green-up, when bear densities were higher (supporting P 2 ). In addition, reindeer generally moved faster at higher bear densities (supporting P 3 ) throughout the growing season, but most consistently in forest habitat.
Earlier studies have shown that ungulates may shift to less profitable habitats to reduce predation risk (Lima andDill 1990, Godvik et al. 2009).
Predation risk may also trigger migration, that is, selection of calving areas far away from predators (Bergerud et al. 1984, Bergerud and Page 1987, Hebblewhite and Merrill 2009, or postpone migration to assumed better quality foraging areas due to higher predation risk (Edwards 1983). Likewise, our results suggest that reindeer females shifted to less valuable habitats for foraging, deviating from following the green wave, when predation risk from brown bears increased. Reindeer may utilize higher elevations and areas with snow cover, where spring phenology is delayed, to reduce predation pressure, as documented in systems with long-term coexistence of Rangifer and large carnivores (Bergerud and Page 1987, Barten et al. 2001, Gustine et al. 2006, Pinard et al. 2012. Animal response to predation risk and human disturbance is comparable (Frid and Dill 2002). Recent studies suggest that in addition to alteration of movement routes, disturbance can cause migratory ungulates to increase movement speed and reduce stopover time (Lendrum et al. 2013, Sawyer et al. 2013, Skarin et al. 2015. For example, mule deer spend 95% of the time during migration in stopover areas taking advantage of elevational gradients Notes: Reference for reindeer herding district habitat is forest. Period represents the two sub-periods within the growing season: calving (11 May-9 June) and post-calving (10 June-31 August). The calving period is reference. SE, standard error. Random-effects standard deviations: year = 2.09 (a) and 3.71 (b); herding district = 0.95 (a) and 1.35 (b). N obs = 815. in forage quality (Sawyer et al. 2013). Thus, reduction in stopover time can potentially lead to less optimal foraging along the migration path, whereas increased movement rates can lead to phenological mismatch both along the route and when arriving at the seasonal range. The increased movement speed of reindeer with higher bear density in our study populations indicates that increased speed may have contributed to the greater phenological mismatch experienced by reindeer at higher bear densities. The question remains whether reindeer also reduce stopover time as a response to predation risk. Reindeer are considered to be capital breeders (Taillon et al. 2013), implying a strategy where the females rely on body reserves for gestation and early lactation (Stephens et al. 2009, Albon et al. 2017. It is hence possible that they are adapted to cope with low forage quality during the fairly short calving period, and are therefore less prone to adjust to the green wave during this period. We predicted a stronger effect of bear density on access to high-quality forage and movement speed during the calving season, as the actual predation risk peaks during this period (Sivertsen 2017). During post-calving, the calf is usually large enough to escape predation from brown bears (Adams et al. 1995), and brown bears likely concentrate more on other food items. We therefore expected reindeer to change foraging strategy in relation to predation risk among seasons, to reduce costs of predator avoidance (Pierce et al. 2004), as shown for several deer species under temporal variation in predation risk (Lone et al. 2017). However, the relationship between bear density and total access to the green wave over the growing season, and between bear density and days from peak IRG, did not differ between periods. Several factors may explain why the negative effect of brown bears on forage access and green-wave surfing, and the higher movement speed, persisted in the post-calving period in our study. It could indicate an innate and strong general risk avoidance behavior, as suggested for other ungulates (Byers 1997). Even in populations devoid of predators, females show characteristic risk-reducing behaviors around calving and high fidelity to calving sites (Tremblay et al. 2007). Nonetheless, Latombe et al. (2014) documented adjustments in habitat selection by caribou in response to temporal variations in predator presence, and Barten et al. (2001) reported that caribou females moved down from higher elevations to areas with higher predator densities, when the calf had developed sufficient locomotive skills (around 10 weeks old). An alternative explanation may be that bear density has been confounded with other environmental factors during post-calving. Alpine habitats may experience less insect harassment and have more areas for relief than forest habitats (Helle and Aspi 1984). However, alpine habitats with low insect harassments are Fig. 2. Predicted access to (a) high-quality forage (measured as cumulative instantaneous rate of greenup; CIRG) and (b) predicted average number of days (absolute value) from the day of peak forage production in spring (day of peak instantaneous rate of greenup) in relation to bear density and sub-period. The two sub-periods calving (11 May-9 June) and post-calving (10 June-31 August) have presumed different calf predation risk, with higher risk in the first period. Points are raw residuals. N obs = 815. Notes: Reference for reindeer herding district habitat is forest and for GPS schedule reference is 12 h. Period represents the two sub-periods within the growing season: calving (11 May-9 June) and post-calving (10 June-31 August). The calving period is reference. SE, standard error; df, degrees of freedom. N obs = 45,332. Fig. 3. Predicted mean daily movement speed in relation to Julian day and bear density index, based on a generalized additive model. Predictions are made for the mean bear density within home ranges of all individuals within each reindeer herding district. The reindeer herding district habitat is shown with solid (forested) and dashed (mountainous) lines. The vertical dashed line shows the two sub-periods calving (11 May-9 June) and post-calving (10 June-31 August), with presumed different calf predation risk (higher in the first period). N obs = 45,332. often of low forage quality (i.e., snow patches and glaciers), which might cause or enhance a mismatch with green-up (M€ orschel and Klein 1997, Hagemoen and Reimers 2002, Skarin et al. 2010. Bergerud and Luttich (2003;and earlier works) argue that predation risk is the most important driver of female behavior during calving and that insect avoidance is more important during summer. Thus, different exposure to insect harassment may have contributed to the observed trends of access to forage and movement speed in the post-calving period.
Habitat structure and forest cover are important for understanding prey responses to predation risk (Murray et al. 1995, Mysterud andØstbye 1999). Our setting is unique in spanning over both mountain and forest reindeer populations, although with mountain populations typically wintering in forested areas. In our system, brown bears preferably reside in the forest (May et al. 2008. Reindeer with access to mountain habitats may advance ahead of green-up into alpine calving ranges. Calving higher up might be a tactic to reduce predation by brown bears during calving. In mountain regions, according to reindeer herders' observations, brown bears may overcome their usual avoidance of open areas and take occasional visits from the forest and up to the open calving grounds to catch reindeer calves (Sikku and Torp 2008). In contrast, the forest reindeer are restricted to calving grounds completely overlapping with the brown bear home ranges (Sivertsen et al. 2016, Sivertsen 2017. Reindeer in the forest generally moved faster compared to those in the mountains, which may reflect higher brown bear predation pressure in the forest compared to mountain calving habitats. With only small pockets of safer habitats to maneuver between (Sivertsen 2017), forest reindeer could be forced to move more between forage and cover habitats to hide from predators (Mysterud and Østbye 1999), and flee when encountering a predator. There is also more human infrastructure intersecting the calving and summer ranges of the forest districts, which may lead to increased disturbance and higher movement rates among the reindeer (Skarin et al. 2015). In addition, insect harassment may explain the higher movement rates in the post-calving period (Hagemoen andReimers 2002, Skarin et al. 2010). Insect harassment is also a possible explanation of the gradual increase in daily movement rate from calving until early August, interrupted by annual calf marking in July.
The bear density index used in this study is based on spatial contrasts, but lacks temporal resolution. Scat surveys are done during late summer/fall typically every fifth year, and the reindeer data used in this study are collected from a diverse set of projects already terminated. An optimal measure of bear density would be annual maps or even indexes following the temporal scale of reindeer GPS locations; however, this was not available. Nevertheless, recent analyses suggest that the presence of individual brown bears in reindeer calving ground areas remains relatively stable throughout the growing season . However, since the bear density index interpolates density based on scat surveys without integrating bear habitat selection patterns, the index may overestimate density in habitats less preferred by brown bears, but situated close to areas with high bear density. The calving ground in Njaarke is located in the mountain, thus in a habitat presumably less susceptible to bear predation, compared to the forested calving grounds. This is supported by statements from the reindeer herders in Njaarke, reporting few problems with brown bear predation. However, Njaarke calving ground is located close to areas with high brown bear densities further north, and accordingly, the risk of brown bear predation may be overestimated from the brown bear density index. Also, lakes and rivers are intersecting the area between the calving ground in Njaarke and the core area with high bear density further north, probably causing a barrier for brown bear movements into the calving ground. Such overestimation of bear predation risk from density maps is most likely not a problem in the other herding districts, located in areas with a more uniform distribution of habitats, and less local spatial variation in bear predation risk. We thus consider the bear density index used in this study sufficient to infer on limitations on reindeer green wave tracking due to potential brown bear presence.
Lower access to greenness combined with higher and more variable speed affects the energy budget and should lead to reduced body conditions, as shown for red deer (Bischof et al. 2012) and caribou (Couturier et al. 2009). White (1983) showed that just a small change in forage digestibility (14%) more than doubled the predicted body weight gain in reindeer. In addition, increased energy expenditure is only part of the effect of increased activity. More important is rather the loss of foraging time, which reduces energy intake (Colman et al. 2003). The recolonization of large carnivores in northern ecosystems recently (Chapron et al. 2014) requires a thorough understanding of population dynamics and potentially cascading consequences on ecosystem function. To be able to evaluate the effect of predators on prey population dynamics, it is important to understand both direct and indirect costs of predator presence (Lima 1998, Lima andBednekoff 1999). Our study indicates that increased predation risk is a factor limiting herbivores' ability to surf the green wave, which may result in negative consequences for the herbivores, such as decreased body condition. This adds to a growing literature of indirect costs of predation (Schmitz 1998, Creel andChristianson 2008), with potential repercussion for ecosystem function (Ripple and Beschta 2004, Creel and Christianson 2009, Ripple et al. 2014, and provides a case for the utility of a unification of the frameworks of energy landscapes and landscape of fear (Gallagher et al. 2017).

ACKNOWLEDGMENTS
Inger Maren Rivrud and Therese Ramberg Sivertsen contributed equally to this work. We greatly acknowledge the support of Centre for Advanced Study in Oslo, Norway, that funded and hosted our research project "Climate effects on harvested large mammal populations" during the academic year of 2015/16. Our study was made possible thanks to the sincere cooperation of the Sami reindeer herders in Sirges, Njaarke, Hand€ olsdalen, G€ allivare, Udtja, € Ostra Kikkejaure, and Mal a reindeer herding districts. Data on reindeer numbers for the study years were kindly provided by the Sami Parliament of Sweden.