Environmental factors in migratory route decisions: a case study on Greenlandic Arctic Terns (Sterna paradisaea)

Abstract Identification and characterization of seasonal migration routes and stopover sites has been recognized as important to the conservation of migratory species. This project utilizes multiple regression models including circular-linear regression to identify associations between route choice, travel speed, and environmental preferences using trajectory data of migratory Arctic Terns (Sterna paradisaea) and environmental data obtained through remote-sensing techniques. Results of this study suggest that route choice on the southward post-breeding migration route may be more dependent on underlying environment than the northward postwintering migration route. In contrast, travel speed was variably associated with underlying environment between southward and northward migrations, including several differences regarding the impact of interactions between environmental variables. These results reveal the importance of using multiple metrics in the estimation of spatial resistance and highlight conflicts between the theoretical resistance framework of GIS and movement analysis methods.

which may or may not be environmentally similar [2]. Characterizing the complex distributions of migratory species requires understanding of not only the breeding and wintering ranges, but also migratory pathways and stopover sites [3]. Animals navigating heterogeneous landscapes face difficult route choice decisions where a linear path is not always optimal; environmental differences in space have been shown to affect movement patterns in many species including Ovenbirds (Seiurus aurocapilla) [4], Hedgehogs (Erinaceus europaeus) [5], Elk (Cervus elaphus) [6], Turkey Vultures (Cathartes aura) [7], and Caribou (Rangifer tarandus caribou) [8] to name a few.
For the past decade, cost-distance approaches have been popular to incorporate the impact of environment on optimal movement pathways [5,7,[9][10][11][12][13][14]. Cost-distances, calculated using geographic information systems (GIS) software, use estimates of travel difficulty across landscapes (resistance surfaces) to calculate the amount of effort required to take a given path across geography; a map of optimal travel corridors can then be generated using diffusion models [10,11]. Resistance surfaces are critical parameters in cost-distance calculation; despite this, there is no standard methodology to generate resistance estimates. Of 24 different cost-distance analyses surveyed by Beier et al. [15] only 9 used empirical data to estimate resistance. The remaining 15 studies relied upon combinations of expert opinion and literature reviews; of the 9 empirical studies, 5 derived resistance estimates from ecological niche modeling (ENM), but none utilized actual movement track data obtained from tracking devices.
The development of modern tracking devices and growing availability of environmental data present an opportunity: by analyzing the movement of individual animals in relation to underlying environments, inferences can be made about the movement preferences of entire populations and even species [1,9]. While there has been a recent growing usage of track data to study animal movement in the literature, how to appropriately analyze tracks remains contentious. While friction models

Introduction
Animal movements occur under diverse circumstances and across spatial scales [1], the most spectacular being seasonal migration. Migratory animals may traverse great distances between breeding and wintering ranges, may be appropriate for environmental resistances such as wind, water currents, or terrain ruggedness, they can misinterpret slow movements associated with foraging or stopover sites [7,[16][17][18][19]. Movement analysis approaches, an alternate methodology, provide a more ecological interpretation of track data where increased residence time indicates increased environmental suitability rather than resistance [18,19]. Speed alone does not determine residence time; a tortuous movement path, turning back on itself frequently, also increases residence time [18,19. The contrasting interpretations of identical movement behavior from resistance-conductance GIS models and movement analysis approaches reveal the need to incorporate multiple metrics when analyzing movement data. In this project, using environmental layers from remote sensing, I analyze trajectory data with two approaches, a novel circular-regression method and a linear regression of travel speed. These are applied to investigate movement decisions and environmental preferences of Arctic Terns (Sterna paradisaea).
Arctic Terns are ground-nesting seabirds which exhibit a broad, circum-polar breeding distribution in high northern latitudes and a wintering distribution in high southern latitudes [20]. Previous researchers have studied Arctic Tern migration [7,21,22], but only recently have complete tracking data of the circum-global migration from Greenland, Iceland, the Netherlands, and Alaska to Antarctic waters been collected [20,[23][24][25]. The massive scale of Arctic Tern migration makes it ideal for studying impacts of large-scale processes (e.g. global climate change) upon migration navigation, resistance factors, and timing choices. Egevang et al. [20] hypothesized that regions of high ocean productivity and prevailing favorable wind currents influence Arctic Tern migration. The importance of food resources and upwelling areas has been well documented for Arctic Terns during migration [24], and for other seabirds [26][27][28][29][30]; the relationship between sea surface winds and route choice has also been documented in other seabirds using predefined resistance models [16,17,27]. While both surface winds and food availability are assumed to be important, the degree to which they influence route choice and travel velocity in Arctic Terns has not been fully explored. The goals of this project are therefore to (1) identify associations between underlying environment and route choice of Arctic Terns; (2) investigate the results of models of route choice and travel speed; and (3) develop methods that utilize track data directly without the need for predefined resistance models.

Tracking Data
Arctic Tern migration data were provided on request by Carsten Egevang [20], collected from August 2007 to June 2008. Tracks were obtained using leg-mounted light loggers (Mk14 geolocators, mass 1.4 g; British Antarctic Survey). The dataset included full migration paths for 9 individual birds, 8 from the eastern edge of Greenland (Sand Island; 74° 43′ N, 20° 27′ W) and 1 from Iceland (Flatey Island; 65° 22′ N; 22° 55′ W). For this study, I divided each migratory route into a post-breeding southward component, from the breeding region in Greenland and Iceland to the overwintering region in Antarctica (August-December, Fig. 1, Table 1), and a post-wintering northward component, from Antarctica back towards Greenland and Iceland (April-June, Fig. 1, Table 1).
In each case, birds were trapped as adults in breeding areas and light-level geolocators were attached to their legs. Geolocators documented and stored light curves which were translated into latitude-longitude coordinates at midday and midnight, resulting in a roughly 12-hour temporal resolution in sampling. These data were processed following methods by Philips et al. [20,31]. Light logging geolocators provide coarse estimates of location (~185 km resolution) becoming increasingly inaccurate in variable periods around the equinoxes when day-night lengths are approximately the same across all latitudes [31]. Despite geolocation coarseness, these devices can be used to broadly describe the entire cycle of migration and are the only tracking devices currently usable on small birds (<100 g) on a continental scale. Additionally, while the relationships between consecutive points are greatly influenced by geolocation error, it is not expected that this error will introduce a consistent directional bias [32,33].
I removed tracking data from the overwintering period (December 2007 to April 2008), leaving only migration periods. As light-logging geolocators determine latitude position from the midpoint of the light curve, points recorded in the weeks surrounding equinoxes (from approximately September 11 th -October 7 th , 2007) show large error in latitude measurements, and were removed [20]. This excludes a portion of the mid-Atlantic southward migration, which limits the applicability of the analysis for this portion of post-breeding migration.
For each locality on each track, I calculated travel angles for the direction of travel to the next point and for the shortest path to destination (Fig. 2). I defined the direction to destination for each locality as the direction

Environmental Data
I collected multiple sources of environmental data covering the time period of migration. I used wind-speed data from the Cross-Calibrated Multi-Platform Ocean Surface Wind Velocity dataset containing interpolated wind speed measurements (in m/s) for u-wind (longitudinal component) and v-wind (latitudinal component) at 10 m above sea level [34]. This data product extends globally over oceans at 0.25° spatial resolution and 5 day temporal resolution. I used daytime sea surface temperature (SST) data (°C) from the NOAA Optimum Interpolation 1/4 Degree Daily Sea Surface Temperature Analysis dataset [35], which comprises spatially interpolated measurements for global daytime sea-surface temperatures at 0.25° spatial resolution. I also used a post-processed dataset of ocean net primary productivity (NPP, mg C/m 2 /day, http://www. science.oregonstate.edu/ocean.productivity/index.php) as a proximate measure of food availability. This dataset is derived from a Vertically Generalized Production Model (VGPM) estimating chlorophyll based photosynthetic capacity [36]. This dataset has a spatial resolution of 0.083°, and 8-day temporal resolution. I used each environmental dataset in its native resolution. For each sample point on each migratory trajectory, I used 2 fan-shaped sampling neighborhoods of equal area to extract environmental values, one in the direction of the shortest path to the final destination (nearest land-edge of Antarctica for post-breeding tracks, nearest land-edge Greenland and Iceland for post-wintering tracks) and one in the observed direction to the next point in sequence. All fans radiated from sampling points with 10° interior angles. Because portions of the southward migration tracks were removed because of location error, the size of each set of sampling neighborhoods varied depending on speed; I defined the side length of each sampling neighborhood pair as the distance traveled in 12 hours at the speed calculated for the sample point, equivalent to the temporal resolution of the tracking data (Fig. 2). Varying the size of sampling neighborhoods limited overlap of environmental sampling for points close together and prevented under-sampling of environments for points further apart. I sampled mean environmental variables using a polygon extract operation in the 'geosphere' package in R. I log transformed the extracted NPP values to reduce skewness present in extracted values. I sampled wind as mean u-wind and mean v-wind values. Because wind is a directional variable, I used the scalar projection of the wind speed on the direction of the sampling fan, including only headwinds and tailwinds.

Circular dispersion models
To estimate which environmental variables, if any, presented migration resistance, I created a circularlinear regression model [37] of angular dispersion as a function of the difference between environments sampled for each point. Circular-linear regression assumes that some angular response variable θ is a function of a mean direction μ and a concentration parameter κ following a Von-Mises distribution. A circular-normal Von-Mises variable has following probability density function: where I 0 is the modified Bessel function of order 0. When κ = 0, this function becomes the circular uniform distribution; for κ ≥ 2 the density becomes tightly concentrated around μ and can be approximated by a normal distribution with variance 1/κ [37]. For this project the expected angle of travel is that of the shortest direct path to the destination, therefore μ = 0 and θ is the angular deviation of the observed path from the expected. θ is a measure of deviation from the global axis [38] rather than an axis constructed from a local trend [7]. With this framework, smaller vales of κ indicate greater average deviation from the shortest distance path to the destination. If Arctic Terns are following cost-optimal pathways during migration, deviation from the shortest path to destination may indicate a greater net cost to travel the direct route. Regressing differences in environmental neighborhoods onto θ estimates the impact of differences in environment on route decisions.
The log-likelihood for a circular-normal distribution with known mean and unknown concentration was proposed by Fisher and Lee [37]: where the subscript i corresponds to the i th data point in the regression. In this instance, i is a particular sample point on a migratory track. The predictor environmental variables were link-transformed using an exponential function: κ i = exp(− βx i + α) to wrap the linear predictors of the i th environment x i onto a prediction of the concentration parameter κi with range [0, ∞). x i is a vector containing the differences of the environmental values of the observed path and the shortest path for the i th point. β is the vector of regression coefficients where positive values indicate differences in environments associated with deviation (i.e. environments less resistant than the shortest path) and α is the intercept. I maximized the log-likelihood model proposed by Fisher and Lee [37] using the "optim" function in R to determine the maximum-likelihood value of β.
Regression models were evaluated separately for the northward and southward tracks to evaluate differences in observed travel behavior between pre-breeding and post-breeding migration. All 2-way interaction terms were included in initial models, and were retained only when significant.

Linear velocity models
To contrast the circular-linear model of route direction, I constructed a mixed linear model of travel speed as a response to underlying environment with individual bird included as an intercept-only random effect: V = βx + Bird + α where the independent variable V is travel speed, β is the vector of regression coefficients on the matrix of predictor variables x, and α is the intercept. The linear regression used only the environmental neighborhoods between consecutive points rather than the difference between observed and shortest path used in the circular regression because potential velocities in non-traveled environments are unknown. I mean-centered the predictor variables to reduce possible colinearity between main effects and interaction terms [39]. Initial model evaluation included all main effects from predictor variables, all potential 2-way interactions between effects, and random effect of individuals. I estimated model coefficients using maximum-likelihood with the 'nlme' package in R (version 3.1-120, J. Pinheiro et al., http://CRAN.R-project. org/package=nlme). I used backwards selection against uninformative model parameters by comparing model AIC scores of the full model to subset models. When an included parameter failed to significantly decrease AIC (ΔAIC<2) compared to simpler a model with fewer parameters, it was eliminated. Parameters were sequentially discarded until reducing the number of parameters lead to an increase in AIC [40,41].

Results
In the circular-linear dispersion model, more favorable winds and higher values of NPP were significantly and positively related to increased dispersion on the southward migration route. For the northward migration route, only decreased sea surface temperature was significantly associated with increased dispersion (Table 2). For both northward and southward dispersion models, no interaction effects were found to be significant and were not included in final models.
In the southward linear regression model regarding travel speed, the full model was rejected as the final model despite performing best because it was indistinguishable from intermediate models including fewer parameters (ΔAIC < 2, Table 3). The final model for the southward linear regression included all main effects (wind, SST, NPP) and the interaction between wind and SST. For the northward route, all interaction terms were included in final models; exclusion of any parameters from the full model resulted in an increase in AIC (ΔAIC ≥ 3.263, Table 3). For both southward and northward migration, high values of NPP were negatively associated with travel speed. Increased sea surface temperature was a marginally insignificant variable for both southward (p = 0.081) and northward models (p = 0.057). Wind speed was not a significant predictor of travel speed for either direction. Several interaction effects were significant and dissimilar comparing southward and northward models. For the southward model, only the negative interaction between wind and SST was significant. For the northward model, all interaction effects were significant and positive (Table 4).

Discussion
For the southward migration, circular-linear regression supports the hypothesis that Arctic Terns show preferential route choice based upon available environments. Southward migration paths divert from the shortest paths towards  Table 3. Summary of models included in model selection. ΔAIC is the difference between the lowest AIC and the given model. AIC weight shows the relative evidence for each model given model. Final models are in bold. more favorable winds, and regions of higher NPP. These results support the popular hypothesis that food availability significantly influences the choice of travel route. For northward migration, only colder sea surface temperature was significantly associated with dispersion in route choice, and only marginally. Despite being predicted as a major influence, wind was not significantly associated with route dispersion on the post-wintering northward migration in this analysis. This negative result may be caused by the fact that the shortest direction of travel aligned well with the direction of favorable winds, eliminating any observable signal of selectivity. Additionally, NPP was not a significant predictor of route choice on the northward route, which indicates that migration strategy may differ post-wintering. Arctic Terns travel more quickly on the northward, spring route (mean travel time of ~40 days) than on the southward, autumn route (mean travel time of ~97 days) which is a general trend for migratory birds [42].

Route
Linear regressions of travel speed on environmental characteristics showed a different and highly complex picture. On both southward and northward routes, areas of high NPP were associated with slower travel speeds despite the fact that no association between food availability and route choice could be shown on the northward route using dispersion models. This may be a signal of opportunistic feeding in regions of high productivity. Northward and southward models differed highly in all other significant variables. Complex differences in interaction terms between northward and southward models are notably more difficult to explain with reasonable hypotheses.
While wind alone was not predictive, at least one wind interaction term was significant for both migration legs. It is possible that winds may be a secondary driver of travel speed, serving as a frictional force that mediates or supports choice of speed indirectly. It is apparent that the underlying decision-making process is highly complex and surprisingly variable.
The possible reliance of Arctic Terns on strong favorable winds is a troubling result for the species. Climate change is likely to reduce the strength of Hadley cells, in turn calming ocean winds [43], which may have detrimental impacts on the migration success of Arctic Terns in the future. Other sea bird species, most notably albatrosses, have shown changes in migratory routes and habits in response to changes in wind patterns [44]. It remains to be seen if the response observed in albatrosses will be a consistent trend amongst other species of sea birds, such as Arctic Terns, which have longer migration routes than albatrosses, relying more on large-scale wind patterns. Climate change may alter the shape, size and magnitude of important wind currents including the North Atlantic Oscillation and the East Atlantic pattern [45]. In other species with shorter migration routes, changes in wind patterns have affected breeding phenology, migration timing, and community composition [44,46,47]. This analysis includes only Arctic Terns over the Atlantic Ocean; application of these techniques onto Pacific Arctic Tern data [24,25] would be necessary to support a consistent species wide trend.
More generally, the directional model of environmental suitability presented here provides a data-  [20]. Additionally, this study highlights the risk of assuming that slower travel speeds indicate increased travel difficulty as per traditional GIS resistance/ conductance framework [48]. Arctic Terns in this analysis showed significant decreases in travel speed in regions of high NPP, a result that is not easily explained by resistance. Circular-linear models show considerable promise as an alternative and complementary tool to analyze animal movements. Several limitations should be considered as regards to this analysis, however. Tracking data used in this study were obtained through light-level geolocation, which has notable limitations. Light-level data are susceptible to shading effects resulting in low spatial accuracy [31,33]. Additionally, light-level data cannot decipher latitude for several weeks in proximity to equinoxes, which excluded several weeks of the southward migration from analysis [20]. While I do not expect that geolocation error introduced a consistent directional bias, it is likely that the signal of environment in the regression analyses was greatly reduced by noise. Results of analyzing lightlogger data is best used at a broad spatial scale in the context of large-scale behavior [27]. The usage of circularlinear regression on travel trajectories is currently an exploratory exercise rather than a fully predictive tool. Circular-linear regression does not account for the fact that animal movements are often auto-correlated spatially and temporally; while alternate approaches can be used to account for autocorrelation, namely state-space or hidden-Markov models, incorporating environmental variables into these models is currently limited in scope [49]. Both the circular-linear and linear regressions only construct monotonic relationships between the underlying environment and travel direction or speed. Intuitively, the relationship would hold only in an optimal range, and the relationship may decay or reverse at extreme values. Response curves in physiology and functional mechanics are generally bell-shaped, such that intermediate environments are favored: e.g., no wind is unfavorable, but extremely fast tailwinds generated from storm systems may be detrimental as well [50]. The suitability curve approach has been investigated extensively as regards to distributions and scenopoetic ecological niches [51], and should be a next goal in movement analyses of this kind.

Conclusions
Calculation of resistance surfaces from regression coefficients for animal migrations is a complicated procedure. Resistance varies simultaneously in direction, time, and space and computations involved in projections are time-consuming and expensive. Future efforts should focus on both velocity and route choice to explore future migration events. Areas of reduced velocity have traditionally been used as an indication of highfriction, but may be the product of much more complex environmental decisions than previously hypothesized.
With increased availability of high-quality tracking data and the advent of new data-driven analytic techniques, it is possible to obtain more direct measurements of migration preference. These tools can estimate impacts of climate change on migration and illustrate suitable areas for migratory pathway conservation. In addition, track analysis may provide a link between migratory species and the environment for cross-taxonomic analyses. Track analyses can be generalized to any species with track data; temporal and taxonomic generalization of the track analysis approach is a sensible next step to study general patterns of migratory animal movements.