Introduction

Animal movements and migrations can provide functional connectivity between areas that are separated in geographical space by transporting biomass, genes, and less mobile organisms. This connectivity has wider ecological implications for the species’ population structure, and can also provide the dispersal opportunities for organisms like flowering plants by moving pollen and seeds, or pathogens (Altizer et al. 2011; Bauer and Hoye 2014). Identifying connectivity networks and understanding the contribution of animal movement to such networks is a prime motive in ecology and is pivotal to our understanding of spatial structuring processes.

Establishing whether, how, and when animal movement provides a functional connection in space, however, is not easily achieved. Capture-mark-recapture techniques have revealed much about dispersal capabalities of individual animals, thereby providing a history of observed connectivity between distant patches. Estimates such as maximum observed dispersal distances can be used to infer connectivity networks where movement has not been observed, yet there are limitations to their application (Calabrese and Fagan 2004), as distance alone can be insufficient to explain patch connectivity. Estimates of effective distance between patches that incorporate barriers and facilitations to animal movement can be used to improve predictions of connectivity. Algorithms like least-cost paths (e.g., Ferreras 2001; Graham 2001) and electrical circuit theory (McRae et al. 2008) can, in combination with spatially explicit predictors of landscape resistance to movement, provide environmentally informed estimates of connectivity between patches (e.g. for population genetics Row et al. 2010). Often, however, the animal location data used to inform models used for predicting such resistance surfaces lack a behavioural context. Consequently, these resistance surfaces might not be representative of how animals move through the environment (Keeley et al. 2017).

More recently, use of remote tracking technology on wild animals has provided great insights into how, when, and where animals move (Hussey et al. 2015; Kays et al. 2015). Such data are not only a rich source of information about the movement and behaviour of individuals, but can also reveal actual connectivity between spatially separated areas in great detail. In combination with environmental information about the utilised habitat, movement data can provide detailed insight into habitat connectivity for the observed individuals (Almpanidou et al. 2014). Connectivity estimates derived from observed movement, as for example in fragmented landscapes, have been shown to outperform predictions derived from resistance surfaces (LaPoint et al. 2013). Yet to estimate connectivity the use of animal movement data is not without constraints (Calabrese and Fagan 2004). While the miniaturisation of tracking technologies permits scientists to follow ever more individuals of ever smaller species, the cost and effort associated with animal tracking limit sample size, as well as the spatial and temporal extent of the data that can be collected. Thus, the number of individuals that scientists are realistically able to track will remain minuscule compared to even the most conservative estimates of the numbers of moving animals on this planet. The goal of an increasing number of studies is to utilise the knowledge from few, well-studied individuals to estimate the behaviour at a population level. However, such generalisations are not straightforward, mainly because the movement behaviour of individuals and the observed variation may not be representative of the population (e.g., Austin et al. 2004). Individual decision-making is not only influenced by general species properties, but also by variation between individuals and their needs and the surrounding environmental conditions (Nathan et al. 2008). Any kind of movement behaviour is thus to some extent unique to the individual, and explicit in time, space, and its environmental conditions as well as its ecological context.

The literature published on models developed for capturing animal movement is extensive, and such models have been shown to provide useful and sensible estimates of the behaviour of observed as well as unobserved individuals (e.g., Morales et al. 2004; Codling et al. 2008; Michelot et al. 2017; Péron et al. 2017). Providing sensible hypotheses of the routes that animals take might require the contextualisation of observed movement, and the understanding of how animals utilise environmental features for route decision-making. Consequently, movement models that incorporate resource-selection functions (step-selection functions, e.g. Fortin et al. 2005; Thurfjell et al. 2014) are becoming increasingly popular. Step-selection functions have been shown to yield functional estimates of how environmental features influence an animal’s movement through the landscape (e.g., Richard and Armstrong 2010), and have been used to estimate connectivity between patches (Squires et al. 2013). Such step-selection functions, representing resource selection during actual movement, can be used to derive behaviour-specific predictions for resistance of a landscape to movement. In combination with least-cost paths or circuit theory, these context-aware resistance surfaces provide the means to predict the movement of individuals through the landscape (e.g., Zeller et al. 2014, 2016).

In many cases, animals use series of different movement strategies that change in response to the surrounding environment, or in response to the different needs an animal has for different behaviour or life-history stages. Currently, however, even context-aware approaches used for predicting the movement of unstudied individuals often make the assumption that animals follow a single, constant decision rule. As shown by Zeller et al. (2016), these decision rules are considered to be independent of the supply needs of the individual. We think that realistic movement simulations should not only take the environmental context of movement behaviour into account, but also acknowledge the different movement strategies expressed by a species (see e.g. Morales et al. 2004). One example of such a multi-state movement behaviour with striking differences between states are the stepping-stone-like migrations as performed by many migratory bird species that predominantly use flapping flight for locomotion. Here, we refer to stepping-stone migrations as performed by large waterbirds like ducks and geese covering large distances in fast and non-stop flight and using stopover locations for extended staging periods to replenish their fat reserves. Context-aware, multi-state approaches for simulating animal trajectories are uncommon. An additional difficulty for the simulation of stepping-stone migratory movements, is that detailed knowledge about available stopover sites for staging migrants might be necessary.

Here, we introduce a novel approach that allows for inferring environmentally informed migratory trajectories from a multi-state discrete movement model. Using a conditional movement model specifically designed for generating random trajectories from template empirical trajectories (Technitis et al. 2016), we developed this approach with stepping-stone migrations and similar movement strategies in mind. We extend this movement model to represent the two major states of stepping-stone migrations, the non-stop migratory flights and the staging periods, using a stochastic switch informed by empirical estimates of typical duration of both behaviours. Our multi-state movement model can simulate migratory trajectories that realisticically represent empirically-collected migratory movements by exclusively sampling from empirical distribution functions. We develop a measure of route viability that integrates properties of the simulated trajectory and its environmental context to assess the joint suitability of the simulated migratory route and timing strategy. For stepping-stone migrations, we assume that the quality of stopover sites between the breeding grounds and wintering areas predominantly determines how preferable a certain route might be (Green et al. 2002; Drent et al. 2007). While the migration simulation model and the measure of route viability we introduce here are tailored for our study system, the approach in general is flexible and could be applicable to many other study systems and strategies.

Specifically, we apply the conditional movement model on a pronounced long-distance migrant, the bar-headed goose (Anser indicus, Latham 1790). This species of waterbird occurs in Central Asia and is well known for its incredible performance of crossing the Himalayas during migration. The distribution range of bar-headed geese is characterised by four distinct breeding areas which are mirrored by four distinct wintering areas south of the Himalayas. Previous tracking studies have revealed that large parts of the respective populations migrate from their breeding grounds in Mongolia, northern China and the Tibetan Plateau over the Himalayas to their wintering grounds on the Indian subcontinent (e.g., Bishop et al. 1997; Takekawa et al. 2009; Guo-Gang et al. 2011; Hawkes et al. 2011; Prosser et al. 2011). But while the crossing of the Himalayas has been studied in great detail (Hawkes et al. 2011, 2013; Bishop et al. 2015), less is known about the connectivity between range fragments both within the wintering and breeding range (Takekawa et al. 2009). The bar-headed goose thus provides a suitable study species for our approach. We establish a model for bar-headed goose migrations from previously published tracking data, and simulate migrations of unobserved individuals between all fragments of the species’ distribution range. We assess the viability of these trajectories during several times of year using a segmented habitat suitability model to derive a dynamic migratory connectivity network. To assess whether this migratory connectivity network could serve as a quantitative null hypothesis for bar-headed goose migration, we test our predictions against two very simple hypotheses generated from previously published studies.

Stable isotope analyses suggested that the connectivity within the breeding range of bar-headed geese is relatively high (Bridge et al. 2015), a notion that has been supported by tracking data as well (Cui et al. 2010). In the wintering range, however, relatively few movements have been observed (Kalra et al. 2011). Based on these findings (Cui et al. 2010; Kalra et al. 2011; Bridge et al. 2015), we expect to find a higher overall viability of trajectories between the fragments of the breeding range than within the wintering range. We further predict that on average, the temporal variation in viability of simulated migratory routes within the breeding grounds should be higher than within the wintering grounds. Overall, we would like to introduce a new approach for deriving environmentally informed quantitative null hypotheses for animal movement which can be utilised for estimating migratory connectivity based on limited observations (summarised in Fig. 1).

Fig. 1
figure 1

General concept for our approach of environmentally informing simulated stepping-stone migrations: (I) Empirical tracking data are (IIa) used to derive an informed eRTG to simulate conditional movement between sites of interest, and (IIb) combined with environmental correlates to derive predictions of relevant measurements of landscape permeability (here: suitability of stopover sites). (III) Finally, the simulated conditional trajectories are evaluated based on characteristics of the trajectory and permeability using an informed measure of route viability

Full size image

Methods

Tracking data and movement model

Tracking data of bar-headed geese were available to us from a broader disease and migration ecology study implemented by the Food & Agriculture Organization of the United Nations (FAO) and U.S. Geological Survey (USGS). In total, 91 individuals were captured during the years 2007–2009 in several locations: Lake Qinghai in China (hereafter termed “Lake Qinghai”), Chilika Lake and Koonthankulum bird sanctuary in India (hereafter termed “India”), and Terkhiin Tsagaan Lake, Mongolia (hereafter termed “West Mongolia”). All individuals were equipped with ARGOS-GPS tags which were programmed to record the animals’ location every 2 h (ARGOS PTT-100; Microwave Telemetry, Columbia, Maryland, USA). Eighty of the deployed tags collected and transmitted data for \(241 \pm 253\,(\text {mean} \pm \text {SD})\) days. In total, 169, 887 fixes could be acquired over the course of the tracking period (Table 1 and Takekawa et al. 2009; Hawkes et al. 2011). Individuals that were tracked for less than a complete year were excluded from the subsequent analyses, which left a total of 66 individuals (Lake Qinghai: 20, India: 20, West Mongolia: 26). We pooled data from all capture sites for the analyses.

Table 1 A summary of the catching sites and corresponding sample sizes
Full size table

We used the recently developed the empirical Random Trajectory Generator (eRTG, Technitis et al. 2016) to simulate the migrations of unobserved individuals of bar-headed geese. This movement model is conditional, i.e. simulates the movement between two end locations with a fixed number of steps based on a dynamic drift derived from a step-wise joint probability surface. One main advantage of the eRTG is that the trajectories it simulates retain the geometric characteristics of the empirical tracking data (step length, turning angle, as well as covariance and auto-correlation of step length and turning angle), as it relies entirely on empirical distribution functions. Consequently, if a destination cannot be reached within the realms of the empirical distributions of e.g. step lengths and turning angles, the simulation fails rather than forcing the last step towards the destination.

We extended this movement model by incorporating a stochastic switch between the two main states of bar-headed goose migration, non-stop migratory flights (“migratory state”) and movements during staging periods at stopover locations (“stopover state”). We classified the entire tracking data according to the individuals’ movement behaviour to identify these states prior to extracting the empirical distributions functions for the eRTG. First, we clustered the locations in the tracking data using an expectation-maximisation binary clustering algorithm designed for annotating animal movement data (EMbC, Garriga et al. 2016). The EMbC divided the trajectories of bar-headed geese into four behavioural classes (slow speed & low turning angles, slow speed & high tuning angles, high speed & low turning angles, and high speed & high turning angles), which we then re-classified into two behavioural classes, namely high-speed movements (combining the two high speed classes) and low-speed movements (combining the two low speed classes). Within the high-speed behavioural cluster, the average speed between locations was \(8.4 \pm 6.7 \frac{m}{s}\) (mean ± SD) whereas the average speed for the low-speed behavioural cluster was \(0.3 \pm 1.0 \frac{m}{s}\) (mean ± SD). As estimates of speed and turning angle are highly dependent on the sampling rate of the data, we removed those parts of the trajectories that exceeded the average sampling interval of 2 h. Subsequently, we used the low-speed locations for the empirical distribution functions for the stopover state of the two-state eRTG, and the locations classified as high-speed for the empirical distribution functions for the migratory state of the eRTG (see Figure S2). Finally, we derived the step lengths and turning angles from each coherent stretch of data (i.e. only subsequent fixes with a sampling rate of 2 hours). Following this, we calculated the changes in step length and turning angle at a lag of one observation, as well as the covariance between contemporary observations of step length and turning angle. We derived the corresponding empirical distribution functions for both movement states and prepared them for use in the eRTG functions.

Finally, we determined the duration of staging periods, and the duration and cumulative distance of individual migratory legs from the tracking data. We first identified seasonal migration events between breeding and wintering grounds (and vice versa) in the empirical trajectories using the behavioural annotation. We then determined migratory legs (sequential locations classified as migratory state) as well as stopovers (sequential locations classified as stopover state, with a duration \(>12h\)). We used two main proxies to characterise migratory legs, namely cumulative migratory distance as well as duration, and one proxy to characterise staging periods, namely stopover duration. We calculated these proxies for all individuals and migrations, and determined the maximum observed distance (\(\text {dm}_{\text {max}}\)) and duration (\(\text {Tm}_{\text {max}}\)) of a migratory leg. As we did not distinguish between extended staging (e.g. during moult, or after unsuccessful breeding attempts) from use of stopover locations during migration, we calculated the \(95\%\) quantile of the observed stopover durations (\(\text {Ts}_{\text {max}}\)) rather than the maximum.

Simulating a bar-headed goose migration with the two-state eRTG

When simulating a conditional random trajectory between two arbitrary locations a and z, the two-state eRTG initially draws from the distribution functions for the migratory state, producing a fast, directed trajectory. To determine the time available for moving from a to z, we assumed the mean empirical flight speed derived for the migratory state, and calculated the number of required steps accordingly. While simulating the trajectory, after each step modelled by the eRTG, the cumulative distance of the trajectory as well as the duration since the start of the migratory leg were calculated. By using cumulative distance and duration as well as the empirically derived \(\text {dm}_{\text {max}}\) and \(\text {Tm}_{\text {max}}\), our two-state eRTG was based on a binomial experiment with two possible outcomes: switching to the stopover state with a probability of \(p_{ms}\), or resuming migration with a probability of \(1-p_{ms}\). We defined \(p_{ms}\), the transition probability to switch from migratory state to stopover state, as

$$\begin{aligned} p_{ms}(t) = \frac{\sum _{i=0}^{t}(dm)}{dm_{\text {max}}} \times \frac{\sum _{i=0}^{t}(Tm)}{Tm_{\text {max}}} \end{aligned}$$
(1)

where \(\text {dm}\) and \(\text {Tm}\) represent the distance and duration between two consecutive locations during a migratory leg. At step t, the simulation of the migratory movement can switch to the unconditional stopover state, corresponding to a correlated random walk, with a probability of \(p_{ms}(t)\). Likewise, the simulation can switch back from stopover state to migratory state with the probability \(p_{sm}(t)\), which we defined as as

$$\begin{aligned} p_{sm}(t) = \left( \frac{\sum _{i=0}^{t}(Ts)}{Ts_{\text {max}}}\right) ^2 \end{aligned}$$
(2)

where \(\text {Ts}\) represents the duration between two consecutive locations during a stopover. This process is then repeated until the simulation terminates because: either the trajectory reached its destination, or the step-wise joint probability surface did not allow for reaching the destination with the remaining number of steps (resulting in a dead end or zero probability).

Evaluating the plausibility of simulated migrations

We estimated the plausibility of each simulated trajectory, representing a unique migratory route, using a measure we called route viability \(\Phi\) aimed to integrate the ecological context into the movement simulations. We developed this measure specifically with the stepping-stone migratory strategy of bar-headed geese or similar species in mind, and it is defined by the time spent in migratory mode, the time spent at stopover sites, and the habitat suitability of the respective utilised stopover sites. For this specific measure of route viability, we made two main assumptions: (1), it is desirable to reach the destinations quickly, i.e. staging at a stopover site comes at the cost of delaying migration, and (2), the cost imposed by delaying migration is inversely-proportional to the quality of the stopover site, i.e. the use of superior stopover sites can counterbalance the delay. Our argument for these assumptions is that during spring migration, the arrival at the breeding grounds needs to be well-timed with the phenology of their major food resources (Bauer et al. 2008). Furthermore, the quality of stopover sites has been shown to be of crucial importance for other species of geese with similar migratory strategies (Green et al. 2002; Drent et al. 2007).

Each simulated multi-state trajectory between two arbitrary locations a and z can be characterised by a total migration duration \(\tau _{a,z}\), which consists of the total flight time \(\tau _{M,a,z}\) and the total staging time at stopover sites \(\tau _{S,a,z}\). The total flight time \(\tau _{M,a,z}\) is the sum of the time spent flying during each migratory leg l, and is thus \(\tau _{M,a,z} \, = \sum _{l=0}^{n}t_M(l)\), with \(t_M(l)\) corresponding to the time spent flying during migratory leg l. Similarly, the total staging time \(\tau _{S,a,z}\) consists of the staging times at all visited stopover sites, corresponding to \(\tau _{S,a,z} \, = \sum _{k=0}^{n}t_S(k)\), where \(t_S(k)\) amounts to the staging time at stopover site k. For our metric of route viability, we will consider the time spent staging at stopover locations \(\tau _{S,a,z}\) as a delay compared to the time spent in flight. This delay is, however, mediated by the benefit b an individual gains at the stopover site from replenishing its fat reserves. We define this benefit gained by staying at stopover site k, b(k), as proportional to the time spent at site k, \(t_S(k)\), and the habitat suitability of site k, S(k). This habitat suitability S should range between [0, 1], which allows our measure of route viability to range between [0, 1] as well. We further assume the effects of several sequential stopovers to be cumulative, and thus define the total benefit of a migratory trajectory between locations a and z with n stopovers as \(B_{a,z} \, = \sum _{k=0}^{n}S(k) \times t_{S}(k)\). Finally, we define the route viability \(\Phi _{a,z}\) of any trajectory between a and z as:

$$\begin{aligned} \Phi _{a,z} \, = \frac{\tau _{M,a,z}}{\tau _{M,a,z} + \tau _{S,a,z} - B_{a,z}} \, = \frac{\tau _{M,a,z}}{\tau _{a,z} - B_{a,z}} \end{aligned}$$
(3)

Thus, the viability of a trajectory with no stopovers and a trajectory with stopovers of the highest possible quality (\(S(k)=1\)) will be equal, and is defined solely by the time the individual spent in migratory state (\(\Phi _{a,z}=1\)). For trajectories with stopovers in less than optimal sites, however, the viability of trajectories is relative to both the staging duration and quality of stopover sites, and should take values of \(\frac{\tau _{M,a,z}}{\tau _{a,z}}< \Phi _{a,z} < 1\). Using this metric, we assessed simulated trajectories in a way that is biologically meaningful for bar-headed geese. In the next section, we detail how we calculated the route viability \(\Phi\) for each simulated migration.

A migratory connectivity network for bar-headed geese

We simulated migrations of bar-headed geese within the native range of the species which naturally occurs in Central Asia (68–107\(^{\circ }\)N , 9–52\(^{\circ }\)E). According to BirdLife International and NatureServe (2013), both the breeding and wintering range are separated into four distinct range fragments (see also Figure S1), with minimum distances between range fragments ranging from 79 km to 2884 km. For this study, we investigated how well, in terms of an environmentally informed measure of route viability and the number of stopovers required to reach a range fragment, these range fragments can be connected by simulated migrations of bar-headed geese.

To choose start- and endpoints for the simulated migrations, we sampled ten random locations from each of the range fragments indicated in the distribution data provided by BirdLife International and NatureServe (2013). We simulated 1000 trajectories for all pairs of range fragments (100 trajectories per location pair) and counted the number of successes (trajectories reach the destination) and failures (trajectories terminate in a dead end). We proceeded to calculate the viability of simulated routes in the following way: Initially, we determined the total duration of the trajectory between locations a and z, \(\tau _{a,z}\), the number of stopover sites used, \(n_{a,z}\), as well as the time spent at each stopover site, \(t_{S}(k)\), for each of the total \(n_{a,z}\) stopovers (corresponding to the number of steps multiplied with the location interval of 2 h). We determined the habitat suitability of stopover locations S(k) using habitat suitability landscapes for bar-headed geese during five periods of the year (see Figure S3): winter/early spring (mid-November–February), mid-spring (mid March–mid April), late spring/summer (mid April–mid August), early autumn (mid August–mid September), and late autumn (mid September–mid November). We identified these periods using a segmentation by habitat use (van Toor et al. 2016, for details see Section A in the Electronic Supplementary Material (ESM)). The segmentation-by-habitat-use procedure uses animal location data and associated environmental information to identify time periods for which habitat use is consistent. Habitat suitability models derived for these time periods should thus reflect differences in habitat use by bar-headed geese throughout the year. We used time series of remotely sensed environmental information and Random Forest models (Breiman 2001) to derive habitat suitability models corresponding to these five time periods, and predicted the corresponding habitat suitability landscapes (Section A in the ESM). Following the prediction of habitat suitability landscapes for winter/early spring, mid-spring, late spring/summer, early autumn, and late autumn, we annotated all stopover state locations of the simulated trajectories with the corresponding habitat suitability. We then calculated the benefit b gained by using a stopover location k using the mean suitability for each of the stopover locations, S(k), and the duration spent at stopover locations, \(\tau _S(k)\).

To calculate the route viability \(\Phi _{a,z}\), we also required an estimate for duration of migration if a simulation were exclusively using the migratory state \(\tau _{M,a,z}\), without the utilisation of stopover sites. We used a simple linear model to predict flight time as a function of geographic distance which we trained on the empirical data derived from the migratory legs (see Section B in the ESM for details). By basing the linear model on the empirical migratory legs rather than mean flight speed, the estimate for \(\tau _M\) retains the inherent tortuosity of waterbird migrations. For each simulated trajectory, we then calculated the geographic distance between its start- and endpoint, and predicted the expected flight time \(\tau _{M,a,z}\). Finally, we calculated route viability \(\Phi _{a,z}\) for all trajectories using Eq. 3, repeating the process for each of the five suitability landscapes derived from the segmentation by habitat use. This resulted in five different values of \(\Phi _{a,z}\) for every simulated trajectory, corresponding to winter/early spring, mid-spring, late spring/summer, early autumn, and late autumn, respectively.

Calculating migratory connectivity as average route viability

We calculated migratory connectivity between range fragments as the average route viability \(\Phi _{\text {avg.}}\) of all trajectories connecting two range fragments. We calculated this average by using non-parametric bootstrapping on the median route viability \(\Phi _{\text {avg.}}\) (using 1000 replicates), and also computed the corresponding \(95\%\) confidence intervals (CI) of the median route viability \(\Phi _{\text {avg.}}\). We did this for each of the five time periods represented in the suitability landscapes, and also computed an overall migratory connectivity by averaging all five habitat suitability values for each stopover site prior to calculating \(\Phi\).

We wanted to compare migratory connectivity within the breeding range and migratory connectivity in the wintering range to test our first hypothesis stating that migratory connectivity should be higher within the breeding range. To do so, we differentiated between route viability among breeding range fragments (\(\Phi_{\text{breeding}}\)), among the wintering areas (\(\Phi_{\text{wintering}}\)), and between breeding and wintering range fragments (\(\Phi_{\text{mixed}}\)). We computed the median and 95% CIs of route viability with non-parametric bootstrapping with 1000 replicates, using the average habitat suitability of all five suitability landscapes for all trajectories within the breeding range, all trajectories in the wintering range, and all trajectories connecting breeding range fragments with wintering range fragments.

To test our second hypothesis, stating that variation in migratory connectivity throughout the year should be higher in the breeding range than in the wintering range, we calculated the standard deviation of route viability for the five suitability landscapes in the breeding range and in the wintering range. We did this by again differentiating trajectories in the wintering range, trajectories in the breeding range, and trajectories connecting breeding range fragments with wintering range fragments. We computed route viability \(\Phi\) for each of the five suitability landscapes for all trajectories, and pooled the corresponding values for \(\Phi _{\text {late \, winter/early\, spring}}\), \(\Phi _{\text {mid-spring}}\), \(\Phi _{\text {late \, spring/summer}}\), \(\Phi _{\text {early \,autumn}}\), and \(\Phi _{\text {late\, autumn}}\) for the wintering range, for the breeding range, and for trajectories connecting breeding range fragments with wintering range fragments separately. We then used a non-parametric bootstrapping (1,000 replicates) on the standard deviation over the five time periods, and determined the corresponding 95% CIs on the standard deviation.

Calculating route viability for empirical migrations

Following the above described procedure, we annotated the stopover locations of empirical migrations with the habitat suitability of the corresponding time period, and calculated the route viability for these migratory trajectories in the same way as described above. We then used non-parametric bootstrapping on the median route viability for all empirical migrations (\(\Phi _{\text {emp., total}}\)), only spring migrations (\(\Phi _{\text {emp., spring}}\)) and only autumn migrations (\(\Phi _{\text {emp.,autumn}}\)), and computed 95% CIs for the median of \(\Phi _{\text {emp., total}}\), \(\Phi _{\text {emp., spring}}\), and \(\Phi _{\text {emp., autumn}}\).

Results

Route viability of empirical and simulated migrations

The simulations resulted in a total of 30,730 simulated trajectories, of which 8945 trajectories connected breeding range fragments (simulation success rate: \(74.5\%\)), 5393 trajectories connected wintering range fragments (simulation success rate: \(44.9\%\)), and the remaining 16,392 trajectories connected breeding and wintering range fragments (simulation success rate: \(51.2\%\); see Figure S4). While all these trajectories were successful in connecting origin and destination (i.e. did not result in a dead end), they differed profoundly in their route viability \(\Phi _{\text {simulated}}\), which ranged between 0.014 and 0.59. We found that simulated migrations had a higher route viability for late spring and summer than for autumn (Fig. 2).

The range of route viability for simulated migrations was comparable to that of the empirical migrations (\(\Phi _{\text {emp., total}}\): 0.01–0.38). Overall, we found that route viability of empirical migrations was higher for spring migrations (\(\Phi _{\text {emp., spring}}\): [0.0614; 0.1070]; 95% CIs on the median) than for autumn migrations (\(\Phi _{\text {emp., autumn}}\): [0.0270; 0.0514]; 95% CIs on the median). This was caused both by differences in the habitat suitability of utilised stopover locations and by differences in migration duration between spring and autumn migrations. We found that bar-headed geese on average stayed longer at stopover locations during autumn than during spring migrations (spring: \(6.8 \pm 14.2\) days, autumn: \(11.8 \pm 12.2\) days; mean ± SD).

Fig. 2
figure 2

The route viability \(\Phi\) of empirical and simulated migrations. Here we show \(\Phi\) for spring and autumn migrations, as well as the \(\Phi\) for the simulated trajectories across all five suitability landscapes (Seg. 1: winter/early spring, Seg. 2: mid-spring, Seg. 3: late spring/summer, Seg. 4: early autumn, Seg. 5: late autumn). The black bars show the 95% CIs for the respective medians, and the grey dots and violin plots show the observed (empirical trajectories) and route viability densities (simulated trajectories)

Full size image

Migratory connectivity network informed by route viability

We separated the simulated trajectories into movements within the breeding range, movements within the wintering range, and movements resembling seasonal migrations between the breeding and wintering range. Here, we found that viability of trajectories was highest within the breeding range (95% CIs for median \(\Phi _{\text {breeding}}: [0.0676;0.0684]\); 95%-quantiles of median \(\Phi _{\text {breeding}}:[0.1469;0.1546]\)), and lowest within the wintering range (95% CIs for median \(\Phi _{\text {wintering}}:[0.0590; 0.0596]\); 95%-quantile of median \(\Phi _{\text {wintering}}:[0.1090;0.1147]\)), predicting that movements between range fragments should occur more often within the breeding than in the wintering areas. The median route viability for migrations between breeding and wintering range fragments was intermediate (95% CIs for median \(\Phi _{\text {mixed}}:[0.0618; 0.0622]\); 95%-quantile of median \(\Phi _{\text {mixed}}:[0.1224;0.1296]\)). These patterns are reflected in the simplified network of average migratory connectivity \(\Phi _{\text {avg.}}\) (Fig. 3). We also identified the single trajectory with the maximum route viability between range fragments rather than the median (Figure S5). This network of maximum migratory connectivity shows that migrations that connect the breeding and wintering ranges have the highest route viability. Finally, the number of stopover locations of movements was proportional to the geographic distance between range fragments (Figure S6).

Fig. 3
figure 3

The median route viability \(\Phi\) between range fragments of bar-headed geese. We summarised \(\Phi\) for all pairwise range fragment trajectories using the median route viability. The thickness of edges represents the sample size. Blue polygons show the native breeding area of the species. Green polygons show the native wintering range. Long edges are curved for sake of visibility

Full size image

Temporal variability of migratory connectivity

We found that the spatial patterns of migratory connectivity varied across the five habitat suitability landscapes representing five periods of consistent habitat suitability (Fig. 4; see also Figure S3 for details on the temporal correspondence of the time periods). For the suitability landscapes derived for winter/early spring, mid spring, and late spring/summer, the estimated connectivity predicted that bar-headed goose migrations are most likely to occur between the wintering and breeding range, and within the breeding range. For early autumn, connectivity patterns predicted that movement is most likely between breeding and wintering areas. For late autumn, we observed connectivity within the wintering range of the species. We calculated the 95% CIs for the overall migratory connectivity values for each time period (Fig. 2), which predicted the highest median route viability for the periods from winter/early spring (mid-November–February) as well as from late spring/summer (mid April–mid August). We also compared the standard deviation of route viability across suitability landscapes and found the highest variation for the breeding range (95% CIs for s.d. of \(\Phi _{\text {breeding}}\): [0.0124; 0.0133]) and the lowest variation for the wintering range (95% CIs for s.d. of \(\Phi _{\text {wintering}}: [0.0041;0.0046]\)). Again, the trajectories between breeding and wintering range fragments showed intermediate values (95% CIs for s.d. of \(\Phi _{\text {mixed}}:[0.0084;0.0089]\)).

Fig. 4
figure 4

Temporal dynamics of the route viability \(\Phi\). Here we show the predicted movements for each of the five suitability landscapes separately. The visible edges of the network have a median route viability \(\Phi\) that is higher than 75% of the route viability for the complete network. The respective time periods associated to these networks is displayed in Figure S3

Full size image

Discussion and conclusions

Using tracking data of bar-headed geese and the empirical Random Trajectory Generator (eRTG), we were able to successfully develop a model that simulates the migratory movements of bar-headed geese. Our extension of the eRTG with a stochastic switch between a migratory state and a stopover state was sufficient to capture the overall migratory strategy of this species. With this model for bar-headed goose migrations, we inferred the migrations of unobserved individuals between all fragments of the species’ distribution range, and used an environmentally-informed measure of route viability to derive average estimates of migratory connectivity between range fragments. We put this simplified predictive network of migratory connectivity to a simple test using predictions derived from the literature. Indeed, we found that the average route viability, as an indicator of migratory connectivity, was higher within the species’ breeding range (\(\Phi _{\text {breeding}}\)) than in the wintering areas (\(\Phi _{\text {wintering}}\)), confirming the expectations from the literature (Cui et al. 2010; Kalra et al. 2011; Bridge et al. 2015). While bar-headed geese are thought to be philopatric to their breeding grounds (Takekawa et al. 2009), the post-breeding period seems to be a time of great individual variability and extensive movements (Cui et al. 2010). This flexibility in long-distance movements has also been observed for other Anatidae species (e.g., Gehrold et al. 2014), and due to the temporary flightlessness during moult the choice of suitable moulting sites is critical to many waterfowl species. As the average route viability within the breeding range and during the summer months is high, unsuccessful breeders and individuals in the post-breeding period may not be not limited by sufficiently suitable stopover locations when moving between breeding range fragments. Furthermore, our results confirmed that the temporal variability of migratory connectivity was higher in the breeding areas north of the Himalayas than in the subtropical wintering areas.

Simulating trajectories with multiple movement states and an element of randomness can be useful to infer the movements of unobserved individuals. Here, we simulated migratory trajectories under the assumption that the movements of the tracked individuals are similar to those of other individuals. While this limits the informative value of an estimate of migratory connectivity, through repeated simulations it may be possible to explore the routes of unobserved individuals according to the movement behaviour observed in empirical data on the same temporal scale. Indeed, average route viability corroborated previous studies on the within-range movements for bar-headed geese even without additional filtering. Consequently, in combination with relevant environmental and ecological information, the simulation of unobserved migrations using a model like our two-state eRTG may provide a sensible and quantitative null hypothesis for the migrations of bar-headed geese or species with similar strategies. While we determined the route viability using only the habitat suitability of the stopover locations and measures of migratory duration, other correlates such as wind support or altitude profile may easily be incorporated for the migratory state. Similarly, the transition probabilities that mediate the switch between movement modes may be extended to include environmental conditions. In general, our stochastic switch performed reasonably well in replicating the movement behaviour observed from recorded tracks. We used simple functions to determine transition probabilities due to the long fix interval (2 h) and the amount of missed fixes in the data. If a larger sample size were available, the functions we used (see Eqs. 1, 2) could be replaced by a probability distribution function that more adequately represents the decision-making of bar-headed geese. Alternatively, algorithms such as state-space models could be integrated to simulate animal movement with a more complex configuration of movement states (Morales et al. 2004; Patterson et al. 2008). With a few modifications specific to the species of interest, the approach described in this study could be adapted for other scenarios of animal movement. One important application for our approach could be to support capture-mark-recapture data, especially when tracking data for multiple individuals are hard to acquire. Simulations from a multi-state movement model informed by the movements of a few representative individuals could be used to infer alternative routes connecting the re-sightings of individually marked animals. A corresponding relevant measure of route viability could then be used to explore alternative strategies from an ecologically informed perspective. In such a study, it could also be of interest to use Bayesian approaches to approximate ideal routes using the environmental context.

Furthermore, our results highlight the importance of integrating temporal changes in habitat use of moving animals into measures of landscape connectivity. Zeigler and Fagan (2014) argue that the ecological function of landscape connectivity through animal movement is not only determined by where, but also when the environment provides the conditions that allow an individual to move from a to z. In our study, estimates of migratory connectivity were affected by changes in the predicted habitat suitability of stopover locations, whereas in other cases, changes in vegetation density throughout the year or the temporary freezing of waterbodies can be imagined to change connectivity between distant sites. Using time series of environmental information in combination with an approach that segments a species utilisation of the environment for moving, as shown here, could help with the identification of temporal patterns of landscape connectivity. Accounting for such temporal changes in connectivity could also help to better understand how for example, diseases can spread through through a metapopulation (such as white-nose syndrome, Blehert et al. 2009; Turner et al. 2011, or Influenza A viruses in birds, Gaidet et al. 2010; Newman et al. 2012).

Overall, models that incorporate a species’ movement behaviour and its utilisation of the environment can provide sensible estimates for landscape connectivity. This approach possibly provides the basis for a wider range of applications, for example the estimation of seed or pathogen dispersal on a population level. Our approach provides a starting point for complementing tracking efforts with ecologically relevant estimates of a species’ potential to migrate through a landscape and act as a link between patches, populations, and ecosystems.