Introduction

Many animal populations benefit from seasonal shifts in resource distribution by migrating between seasonally distinct ranges (Fryxell and Sinclair 1988; Avgar et al. 2013). Theoretical studies demonstrate that migration is an evolutionarily stable strategy in environments where areas with the highest population net-growth shift markedly between seasons (Holt and Fryxell 2011; Fryxell and Holt 2013), which results in partial or complete migratory populations (Chapman et al. 2011). Migration is thus an adaptive and spatially extensive utilization of resources that enhances individual fitness (Mariani et al. 2016). The large-scale redistribution of individuals following these long-distance migrations results in an ideal-free distribution within each season (Fretwell and Lucas 1969, 1972), but only when migration has negligible cost. Currently, snow deposition and seasonal resource pulses driving migrations are rapidly shifting due to climate change (Middleton et al. 2013; Brown and Mote 2009). In addition, increasing anthropogenic movement barriers are turning terrestrial migration into a globally threatened phenomenon (Berger 2004; Bolger et al. 2008; Harris et al. 2009; Wilcove and Wikelski 2008).

Anthropogenic infrastructure development that has blocked traditional migration corridors causes declines in abundance of the affected migratory populations (Berger 2004; Bolger et al. 2008; Harris et al. 2009; Wilcove and Wikelski 2008). Barriers to movement can affect the proportion of migratory individuals in a partial migratory population mainly in two ways (Shepard et al. 2008). First, barriers can increase the cost of migration by increasing mortality risk or energy expenditure associated with crossing the barrier. For instance, collisions with vehicles during road crossings are a common source of mortality for terrestrial animals (Trombulak and Frissell 2000; Shepard et al. 2008). Second, a structure can also act as a barrier by altering animal behavior, i.e., an individual may avoid a barrier resulting in a decreased crossing probability without necessarily a direct effect on its survival (Beyer et al. 2016; Shepard et al. 2008; Jaeger and Fahrig 2004). The latter can be a non-adaptive behavioral response, when mortality risk is low and benefits of migration are high.

At northern latitudes, seasonal contraction of ungulate ranges is driven by snow deposition in fall increasing mortality risk in their highland ranges, while range expansion in the spring is driven by “green waves” of newly emerging, high-quality forage enhancing growth and recruitment of migratory individuals (Bischof et al. 2012; Aikens et al. 2017). Although geographic variation (e.g., altitudinal) in snow depth often drives migratory behavior in northern latitudes, other sources of geographic variation may motivate animals to migrate (Fryxell and Sinclair 1988). Snow deposition is very sensitive to global warming; even moderate emission scenarios alter snow cover markedly, which alter future habitat suitability of migratory ungulates (Rivrud et al. 2019). A reduction in snow cover would expand suitable year-round ranges and lead to a decreasing benefit of migration in northern latitudes. Such changing migratory patterns have far-reaching implications, as migratory animals transport nutrients, energy, and other organisms between disparate locations (Bauer and Hoye 2014).

Many migratory animal populations provide key resources for indigenous people and rural economies often through harvest and have been called “mobile agent–based ecosystem services” (Kremen et al. 2007; Semmens et al. 2011). Optimal harvest strategies of fluctuating populations has received extensive focus both from theoretical (Lande and Engen 1994; Lande et al. 2003) and empirical studies (Forrest and Walters 2009; Sæther et al. 2001); however, these models typically assume a single resident population (Ling and Milner-Gulland 2008). Recent research efforts have extended such approaches to include animal movements into their harvest models, for instance, dispersal movements in a metapopulation (e.g., Brooks and Lebreton 2001; Costello and Polasky 2008; Neubert 2003; Salinas et al. 2005; Sanchirico et al. 2001) or dispersal movements with protected areas (e.g., Křivan and Jana 2015). Nilsen et al. (2009) found that high levels of migration resulted in a coupling of seasonal ranges, which led to an economic loss, when harvest on both ranges is managed independently.

We investigate how population dynamics and sustainable harvest of (partial) migratory populations are affected by barriers to migration and changes in seasonality. A changing seasonality will affect the relative benefits of being migrant versus resident, while barriers reduce the migration probability or inflict a cost of migration that leads to a violation of the ideal-free distribution across seasons. We build upon the migratory population model used by Fryxell and Holt (2013), but replaced the migratory morphs in that model with behaviorally plastic ideal-free migration (Mariani et al. 2016) and allowed for population harvest. This allowed us to investigate the effects of changing seasonality and movement barriers through both fitness cost and non-adaptive barrier avoidance on population dynamics and harvest strategies.

The model

Partial migratory population dynamics

We used the population model presented by Fryxell and Holt (2013), where the dynamics of population size follow a Ricker model (Ricker 1954), with scaled densities (i.e., HCode \(N = N^{\prime }/K^{\prime }\), where \(N^{\prime }\) is the unscaled population size in spring just prior to the reproduction season and \(K^{\prime }\) is the population size at which on average each individual replaces itself during summer. Table 1 provides an overview of all symbols). Following Fryxell and Holt (2013), we assumed recruitment (r) during summer to be density dependent and mortality during winter (μ) to be density independent. Both assumptions are supported by many empirical studies of ungulates (Sæther 1997; Gaillard et al. 2000).

Table 1 Symbols
Full size table

Assuming non-overlapping generations and using the Ricker formula to represent episodes of summer reproduction, the multiplicative growth rate equals \(\exp (r[1-N^{\prime }/K^{\prime }])=\exp (r[1-N])\), where er is the maximum per capita recruitment during summer. Assuming density-independent winter survival probability e−μ, the number of animals after 1 year at the end of winter is calculated as follows:

$$ N (t+1) = N (t) \exp\Big(r \big[ 1 - N (t) \big] - \mu \Big) $$
(1)

The (scaled) equilibrium population size (N∗) for this model is (i.e., where r[1 − N∗] − μ = 0, from Eq. 1):

$$ N^{*} = 1-\frac{\mu}{r} = 1 - \phi $$
(2)

The equilibrium population size is thus less than the summer carrying capacity; the difference is given by the ratio (ϕ) between mortality (μ) and recruitment (r); when ϕ ≥ 1, a population is not viable year-round.

The seasonal fitness in an area will fluctuate between \(\exp (r[1 - N])\) in summer and \(\exp (- \mu )\) in winter (Fryxell and Holt 2013). Similar to Mariani et al. (2016), we define the seasonality (𝜃) of an area as the difference between these seasonal fitnesses. Due to the density dependence of the summer fitness, this definition of seasonality is also density dependent. At equilibrium, the seasonality will depend solely on the winter fitness, which can be approximated using a first order Taylor expansion with a very simple expression:

$$ \theta = \exp(r[1-N^{*}])-\exp(-\mu) = e^{\mu}-e^{-\mu} \approx 2 \mu $$
(3)

Although at high density the amplitude of the variation in seasonal fitness depends only on the winter mortality rate, it is not intuitive to define seasonality in an area exclusively based on the demographic parameters from that area during a single season. We therefore defined seasonality based on the amplitude of the fitness fluctuation at low population density:

$$ \theta = \exp(r[1-0])-\exp(-\mu) = e^{r}-e^{-\mu} $$
(4)

Thus, we define seasonality as the difference between the intrinsic summer recruitment and the winter survival.

Following Fryxell and Holt (2013), we linked two population ranges through the movement of migratory animals, which move with migration probability, m, and cost, c, after the winter season from one range to the other (see Fig. 1 and in Appendix Fig. 7). For convenience, we will refer to these ranges as either lowland (L) or highland (H), since many ungulate migrations follow strong elevation gradients (Mysterud et al. 2011). However, depending on the ecological context other labels could be more appropriate, e.g., southern and northern range for latitudinal migration; our results and conclusions are not limited to the case of altitudinal migration only. The lowland range is characterized by lower winter mortality (i.e., μL ≤ μH). Therefore, we assumed during fall only migration from high- to lowland; the so-called “perverse” migrants (sensu Fryxell and Holt 2013) from low- to highland during fall were not considered. All migrants that moved during spring to the highland range move back in the fall together with their offspring to the lowland range for the winter season.

Fig. 1
figure 1

Illustration of migratory events throughout the year between the two ranges: highlands (left) and lowlands (right). The migration barrier is represented by the road in between. During winter (upper row), only highland residents can be found in the highlands, while both residents and potential migrants share the lowland range. During summer (lower row), the highlands are occupied by highland residents together with migrants from the lowland area, while only lowland residents remain in their range. Note that such resident or migrant individuals may not be present for certain parameter combinations, see main text

Full size image

The number of individuals migrating from the lowland range is \(mN^{\prime }_{L}\) so that the number of individuals contributing to population regulation in the highland range is \(N^{\prime }_{H} + mN^{\prime }_{L}\). Following the population model in Eq. 1, the multiplicative growth rate in the highland range for both residents and migrants equals: \(\exp (r_{H}[1-(N^{\prime }_{H}+mN^{\prime }_{L})/K^{\prime }_{H}])=\exp (r_{H}[1-N_{H}-\frac {K^{\prime }_{L}}{K^{\prime }_{H}}mN_{L}])\). The number of animals after 1 year at the end of winter on the highland range is calculated as follows:

$$ N_{H} (t+1) = N_{H} (t) \exp\Big(r_{H} \big[ 1 - N_{H} (t) - \frac{K^{\prime}_{L}}{K^{\prime}_{H}}mN_{L} (t) \big] - \mu_{H} - {h_{H}^{S}} - {h_{H}^{W}}\Big) $$
(5)

where \(e^{-{h_{H}^{S}}}\) and \(e^{-{h_{H}^{W}}}\) are the probabilities of surviving respectively the summer and winter harvest season in range H. We use subscripts to indicate the harvested range and superscripts to indicate the timing of the harvest (see details in the section on “population harvest” below).

The number of animals at the end of winter in the lowland range, NL, after 1 year is the sum of residents in L and migrants, which returned to L after the summer on H:

$$ \begin{array}{@{}rcl@{}} N_{L} (t+1) = (1-m)N_{L} (t) \exp\Big(r_{L} \big[ 1 - (1-m)N_{L} (t) \big] - \mu_{L}- {h_{L}^{S}} - {h_{L}^{W}}\Big)\\ +mN_{L} (t) \exp\Big(r_{H} \big[ 1 - N_{H} (t) - \frac{K^{\prime}_{L}}{K^{\prime}_{H}}mN_{L} (t) \big] - \mu_{L} - c- {h_{H}^{S}} - {h_{L}^{W}}\Big) \end{array} $$
(6)

where \(e^{-{h_{L}^{S}}}\) and \(e^{-{h_{L}^{W}}}\) are the probabilities of surviving respectively the summer and winter harvest season in the lowland range. For simplicity, we assumed that migration cost (c) is incurred after summer, so that the fraction of migrants surviving winter is \(e^{-\mu _{L}}e^{-c} = e^{-\mu _{L}-c}\). Density-dependent recruitment (r) takes place during summer, which is for migrants in the highland range (i.e., rH) and in the lowland range for those residents (i.e., rL), whereas the density-independent winter mortality (μ) occurs when migrants and residents share their common lowland range.

Population harvest

We defined harvest on each range as the proportion of the population harvested in Eqs. 5 and 6 during either the summer (\({H^{S}_{H}}\) and \({H^{S}_{L}}\)) or the winter season (\({H^{W}_{H}}\) and \({H^{W}_{L}}\)), we assumed that harvest never occurs during both seasons (i.e., if \({H^{S}_{H}} > 0\) or \({H^{S}_{L}} > 0\), then \({H^{W}_{H}} = {H^{W}_{L}} = 0\), or if \({H^{W}_{H}} > 0\) or \({H^{W}_{L}} > 0\), then \({H^{S}_{H}} = {H^{S}_{L}} = 0\)). The proportion of the population in a range (H or L) surviving the harvest season (S or W ) is then \(1 - H^{\text {season}}_{\text {range}} = \exp (-h^{\text {season}}_{\text {range}})\). Whether harvest takes place before or after migration does not affect residents directly (i.e., residents in H experience either \({H^{S}_{H}}\) or \({H^{W}_{H}}\), and residents in L either \({H^{S}_{L}}\) or \({H^{W}_{L}}\)); however, for migrants, the timing of harvest will affect whether they experience the same harvest pressure as highland residents during summer (\({H^{S}_{H}}\)) or as lowland residents during winter (\({H^{W}_{L}}\)).

The long-term equilibrium dynamics of this system, based on Eqs. 5 and 6, is defined in terms of the equilibrium population sizes on the lowland (\(N_{L}^{*}\)) and highland range (\(N_{H}^{*}\)), which are calculated as \(N_{L}^{*} = N_{L}(t) = N_{L}(t+1)\) and \(N_{H}^{*} = N_{H}(t) = N_{H}(t+1)\), respectively. The overall maximum sustainable yield is then given by the harvest proportions (\({H^{S}_{H}}\) and \({H^{S}_{L}}\), or \({H^{W}_{H}}\) and \({H^{W}_{L}}\)) that realize the largest yield (for summer harvest in Eq. 7, and winter harvest in Eq. 8) at equilibrium summed over both ranges. The total long-term yield from summer harvest YS is the summed yield from residents and migrants on the highland range and residents on the lowland range:

$$ Y^{S} = {H^{S}_{H}} (N_{H}^{*} + \frac{K^{\prime}_{L}}{K^{\prime}_{H}}mN_{L}^{*}) + {H^{S}_{L}} \frac{K^{\prime}_{L}}{K^{\prime}_{H}}(1-m)N_{L}^{*} $$
(7)

with the yield scaled by the carrying capacity of range S (\(K^{\prime }_{H}\)). Similarly, the total yield from winter harvest YW is the summed yield from the highland range and the lowland range, with the migrants returned to their lowland range for the winter:

$$ Y^{W} = {H^{W}_{H}} N_{H}^{*} + {H^{W}_{L}} \frac{K^{\prime}_{L}}{K^{\prime}_{H}}N_{L}^{*} $$
(8)

Density-dependent migration with a semi-permeable barrier

The population dynamics on the lowland and highland range are linked through partial seasonal migration. As more animals migrate (i.e., increasing migration probability, m) from the lowland range, the density of animals in the lowland range decreases during summer and the summer fitness of residents in that range increases, whereas the density of animals in the highland range increases during summer and the summer fitness of migrants in the highland range decreases (Fig. 2). We assumed an ideal-free migration strategy (Mariani et al. 2016), where the ideal-free migration probability (\(\hat {m}\)) equalizes the summer fitness for lowland range residents and migrants, calculated from Eq. 6:

$$ \exp \Big(r_{L} (1 - (1-m)N_{L}) - {h_{L}^{S}} \Big) = \exp \Big(r_{H} (1 - N_{H} - \frac{K^{\prime}_{L}}{K^{\prime}_{H}}mN_{L}) - c - {h_{H}^{S}} \Big) $$
(9)

Since this equation is linear in m, the ideal-free migration probability (\(\hat {m}\)) is calculated as follows:

$$ \hat{m} = \frac{r_{H}[1-N_{H}] - r_{L}[1-N_{L}] - ({h_{H}^{S}} - {h_{L}^{S}}) - c}{[r_{L} + r_{H}\frac{K^{\prime}_{L}}{K^{\prime}_{H}}] N_{L}} $$
(10)
Fig. 2
figure 2

The summer fitness of the lowland residents (dashed line) and the migrants (dot-dashed line) as a function of the migration probability. The solid line shows the average summer fitness for both movement tactics. The vertical gray lines indicate the ideal-free migration probability (\(\hat {m}\)) – where the fitness of both tactics is equal – at 0.59, and the migration probability maximizing the average population summer fitness or growth (m∗) at 0.49. The parameters for Eq. 9 were rH = 0.47, rL = 0.33, KL/KH = 1, c = 0, \({h^{S}_{H}}={h^{S}_{L}}\), NH = 0, and NL = 1. It is interesting to note that \(m^{*} < \hat {m}\); thus, maximum growth of the population at the lowland range is realized at a migration probability below the ideal-free migration probability

Full size image

Thus, the ideal-free migration probability increases as the difference in saturation between the high- and lowland ranges increases, weighted by their respective levels of summer recruitment. Moreover, it decreases when the cost of migration increases, and when harvest pressure during summer is higher on the highland range than on the lowland range. We assumed omniscient or “ideal” individuals; in reality, however, an individual’s knowledge of these parameters will be incomplete. The effects of incomplete knowledge on migration behavior and population dynamics are an interesting avenue for future research.

Interestingly, the ideal-free migration probability is not the migration probability that maximizes the average fitness for the lowland range population (Fig. 2). The migration probability that maximizes the average fitness will be lower than the ideal-free migration probability when rH > rL. Instead, the ideal-free migration probability (\(\hat {m}\)) is the one that equalizes the summer fitness on the high- and lowland range (Fig. 2). This solution is therefore a Nash-equilibrium (Osborne 2009), since an individual changing its movement tactic would increase the density in the range it selects and experiences a lower fitness than if it would have stayed (Křivan et al. 2008).

The ideal-free migration probability accounts for the fact that migration can be costly (c > 0); however, it does assume that animals are “free” to move. In addition to increasing the cost of migration (c), barriers can also affect migration behavior through behavioral avoidance of the barrier (i.e., beyond an adaptive response to an increasing migration cost), which further reduces the probability of migration. This non-adaptive behavioral avoidance effect of a barrier (b) will lead to a lower realized migration probability (m in Eqs. 5 and 6) compared to the ideal-free migration probability in the absence of such a barrier (\(\hat {m}\), in Eq. 10):

$$ m = (1-b) \times \hat{m} $$
(11)

Thus, a barrier can affect migration probability either through behavioral avoidance, b, in Eq. 11, or through the migration cost, c, in Eq. 10.

Simulations

We used 𝜃L = 0.5 and 𝜃H = 0.7 (see Eq. 4) as reference values for seasonality and ϕL = ϕH = 0.5 (see Eq. 2) for year-round range performance. The vital rates corresponding to these reference values are rL = 0.31, rH = 0.41, μL = 0.15, and μH = 0.21, similar to those used by Fryxell and Holt (2013). Table 2 shows a summary of the different parameter values used. In addition to this reference with a relatively high year-round range performance, we also provide in Supplementary Figures the results from a scenario with the same seasonality, but with low year-round range performance: ϕL = ϕH = 0.95 (vital rates: rL = 0.25, rH = 0.35, μL = 0.24, and μH = 0.33). We initialized population sizes in both ranges at 0.01 for all simulations (i.e., NL(t = 0) = NH(t = 0) = 0.01); we tested the sensitivity to the initial values by comparing with both NH(t = 0) = 0.9 and NL(t = 0) = 0.9. Our results were not affected by these starting values, unless stated otherwise.

Table 2 Simulation parameter values
Full size table

First, we investigated the role of seasonality on ideal-free migration by varying the seasonality difference between the high- and lowland range (𝜃L = 0.5 and 0.5 ≤ 𝜃H ≤ 0.7) for different year-round population performance (ϕL = ϕH = 0.5, ϕL = 0.5 with ϕH = 0.7, and ϕL = ϕH = 0.95) in the absence of a barrier (c = b = 0). Secondly, we looked into the yield obtained from the total population using the reference values for different harvest probabilities, with harvest taking place either during summer or winter, in the absence of a barrier (c = b = 0). Finally, we investigated how both migration cost (0 ≤ c ≤ 0.15) and behavioral avoidance (0 ≤ b ≤ 1) of barriers influence the population size and the maximum sustainable yield using the reference values.

We used the optim-function with the constrained quasi-Newton algorithm (i.e., “L-BFGS-B,” Byrd et al.1995) in R (R Core Team 2016) to find the maximum yield at equilibrium (i.e., the maximum sustainable yield).

Results

Ideal-free migration in a seasonal environment

The reference parameter values result in 56% of the population migrating and 44% remaining resident year-round in the common lowland range (Fig. 3). However, when the year-round carrying capacity is low (1 − ϕ = 0.05), then 87% of the population migrates and only 13% remains resident. Under both scenarios, residents in the highland range cannot persist because they are out-competed by the migrants, who benefit from the lower mortality on the lowland range (μL < μH).

Fig. 3
figure 3

The population dynamics of the highland residents (dotted line), lowland residents (dashed line), and migrants (dot-dashed line) over 100 years. Both population ranges were initiated at 0.01. The panel on the left shows results with the reference values (Table 2, and on the right reference values with low year-round carrying capacity, ϕL = ϕH = 0.95). Both populations show partial migration. However, as the year-round carrying capacity of each range decreases towards 0, the migration behavior of the population moves from partial towards obligatory

Full size image

More generally, the proportion of migrants increases with increasing differences in seasonality or decreasing year-round carrying capacity of high- and lowland ranges (Fig. 4). The model is capable of producing three main states. First, in the absence of a difference in seasonality (𝜃L = 𝜃H) and year-round carrying capacity (ϕL = ϕH) the proportion of migrants is zero and the population consists simply of two resident populations. Second, when differences exist between ranges in seasonality and/or year-round carrying capacity for a migration tactic to exploit, animals from the lowland range migrate to the highland range and out-compete highland residents (as μL < μH). This results in a partially migratory population on the lowland range, without residents on the highland range (as is the case with the reference parameter values shown in Fig. 3). Third, when the lowland range would not be viable year-round (i.e.. ϕL ≥ 1), the animals on the lowland range become obligatory migrants without residents in either range. Thus, ideal-free migration (\(0 < \hat {m} \leq 1\)) will exist whenever the ranges differ in winter mortality (i.e., μH > μL). The proportion of migrants increases in the population as the year-round carrying capacity on the lowland range decreases (i.e., decrease in 1 − ϕL), or as the recruitment on the highland range increases (i.e., increase in rH). When μH = μL, the proportion of migrants was dependent upon the initial population sizes, because the fitness of migrants and highland range residents is equal when migration is cost-free.

Fig. 4
figure 4

Proportion of migrants at equilibrium as a function of the seasonality difference between both ranges (d𝜃 = 𝜃L − 𝜃H) and year-round carrying capacity of each range (1 − ϕL and 1 − ϕH). The seasonality of the lowland range was fixed (𝜃L = 0.5). For migration to take place, the two ranges need to differ in seasonality and/or carrying capacity, in the absence of both the proportion of migrants will be zero. The proportion of migrants increases as the difference in seasonality between both ranges increases, and as the year-round viability of the lowland range decreases

Full size image

Harvest of a migratory population

The yields from harvesting a partial migratory population are dependent upon the timing of harvest (Fig. 5), and unsurprisingly on the year-round carrying capacity (i.e., 1 − ϕ; Supplementary Fig. 8). Harvest during the summer season, when migrants are on the highland range, gives a higher maximum sustainable yield (MSY), than harvest during the winter season, when migrants are on the lowland range. The winter harvest rate with MSY makes lowland residency not viable (i.e., \(\frac {\mu _{L} + {h_{L}^{W}}}{r_{L}} \geq 1\), with the reference values from Table 2), thereby restricting summer recruitment to the highland range, which results in a lower MSY for winter than summer harvest.

Fig. 5
figure 5

Yield as a function of harvest fraction on the highland range (HH) and lowland range (HL) for summer (panel on the left, HS) and winter harvest (panel on the right, HW). See Table 2 for all parameter values

Full size image

Changes in the yield from summer harvest are most affected by the highland range harvest probability (\({H_{H}^{S}}\); Fig. 5). Beyond intermediate harvest pressures during summer on the lowland range, residency on that range is no longer viable and the population is dominated by migrants; therefore, the yield comes from the harvest of migratory individuals, which during summer are on the highland range. For the same reason, the yield from winter harvest is most affected by the lowland range harvest probability (\({H_{L}^{W}}\)), when migrants are back on their range. Interestingly, when the winter harvest pressure increases on the lowland range, a population of highland range residents can grow due to reduced densities of competing immigrants from the lowlands.

Not surprisingly, both yield and harvest probability at the maximum sustainable harvest are lower when year-round carrying capacity (i.e., 1 − ϕ) on both ranges is lower (Supplementary Fig. 8). At low year-round carrying capacity, the population is a virtually obligatory migrant (Fig. 3), which is reflected in the fact that the yield is largely dependent upon harvest from the range where the migrants are present, i.e., yield from the highlands for summer harvest and yield from the lowlands for winter harvest.

Effects of migration barriers on population size and harvest

Both migration cost (c) and behavioral avoidance (b) of a barrier lead to a reduction in total population size (Fig. 6). The migration cost and behavioral avoidance of a barrier interact negatively (Supplementary Figs. 9 and 10): with an increasing migration cost (c) the behavioral avoidance effect (b) of the barrier decreases, due to the ceiling on the total barrier effect (when the barrier is maximal, migration no longer occurs). Although a behavioral avoidance generally leads to a decrease in the total population size, a small avoidance effect actually increases the total population size (i.e., \(b = 1- \frac {m^{*}}{\hat {m}} \approx 0.2\); Fig. 6d) as it decreases the ideal-free migration probability towards the migration probability that maximizes average fitness (Fig. 2).

Fig. 6
figure 6

Population size (panels on the left ) and harvest yield (panels in the center and on the right) as a function of migration cost (c, top row) and non-adaptive barrier avoidance (b, bottom row). Population size in the absence of harvest is shown for residents in the lowland range (dashed line), residents in the highland range (dotted line), migrants (dot-dashed line), and in total (solid line). The yield is shown on the lowland range (dashed line), highland range (dotted line), and in total (solid line) from harvest with the maximum sustainable yield for summer (center panels) and winter harvest (right hand panels). See Table 2 for parameter values

Full size image

The overall reduction in population size due to barriers is partially buffered by an increase in resident population size, when the year-round carrying capacity is relatively high (i.e., ϕ = 0.5), especially in the highland range. This increase in the highland resident population follows the reduced competition from migrant individuals. When the year-round carrying capacity in each range is low (i.e., ϕ = 0.95), then little buffering will occur and the negative barrier effects on total population size are stronger (Supplementary Figs. 9 and 10).

Also, the maximum sustainable yield (MSY) is negatively affected by both migration cost and behavioral avoidance of barriers (Fig. 6). This decrease is particularly dramatic when the year-round carrying capacity is low (Supplementary Figs. 11 and 12), and little buffering from resident populations takes place.

Paradoxically, barriers can lead to a local increase in yield, despite the overall decrease in the MSY due to barriers. When the year-round carrying capacity is sufficiently high, barriers lead to a local increase in yield from the highland range. Therefore, this barrier-mediated increase in harvest yield from the highland residents results in a more equitable distribution of the yield across both ranges (Fig. 6c and f).

Discussion

We presented a simple model for the dynamics of a density-dependent migratory population to investigate the effects of changes in seasonality and migration barriers on population size and harvest. Our model predicts partial migration in environments where ranges differ in seasonality and/or year-round carrying capacity. Migration barriers—irrespective of whether they increase the cost of migration or are simply avoided—lead to a reduction in population size and harvest yield. However, when year-round carrying capacity is sufficiently high, this reduction is buffered by an increase in the resident populations. This increase in the resident population can result also in a local increase in harvest yield at the expense of a global decrease in yield due to migration barriers.

Our model simplifies the demography of ungulates by omitting sex-, age-, state-dependent migration, etc. Although these different elements will probably affect the finer details of our model’s behavior, our main results stem from the migrants experiencing (a) higher reproduction during summer compared to lowland residents and (b) reduced winter mortality compared to highland residents. The higher summer reproduction resulted from a release of density-dependent competition with lowland residents, whereas the lower winter mortality was independent of density. This seasonal structure was inspired by our choice of ungulates in the northern hemisphere as our model species group (Sæther 1997; Gaillard et al. 2000), being the economically most important terrestrial migratory resource subject to harvest. However, we do expect the general dynamics of our model to be applicable to a broader range of migrating species in environments with a different seasonal structure as long as migrants can have an advantage over the residents on both seasonal ranges. By assuming density dependence only in one season, Eq. 9 becomes linear in m and thus mathematically convenient to derive. However, species may experience density dependence on both ranges. Introducing density-dependent mortality on the lowland range into our model could reduce the benefit that migrants have over highland range residents from returning to the lowland range; thus, under certain parameter values, highland range residents can be expected (see Holt and Fryxell 2011). However, a full exploration of the added complexity of this double density dependence on ideal-free migration was beyond the scope of the present paper.

Global warming and changing seasonality

In an empirical setting, we would only observe winter mortality in the shared range of residents and migrants, while our theoretical setting allows us to explore the effects of the non-observed (expected) mortality of remaining on the highland (recruitment) range. Similar to Mariani et al. (2016), we find density-dependent migration (either partial or obligatory) when seasonal differences exist between both ranges (i.e., rH > rL and/or μH > μL). These results confirm the findings by Holt and Fryxell (2011) and Fryxell and Holt (2013) when considering the evolution of a migratory morph. Thus, in seasonal environments, we would expect migration behavior independently of whether it is an ideal-free behavioral choice or genetically fixed.

A primary driver of seasonal migration of ungulates at northern latitudes is snow deposition leading to range contraction in autumn (Nelson 1995; LeResche 1974; Brazda 1953), since snow depth is a major driver of overwinter mortality (Cederlund and Lindström 1983). The response of snow deposition patterns to global warming is affected by the interaction between precipitation and the temperature threshold around freezing (Brown and Mote 2009); therefore, even moderate climate change scenarios can markedly affect future habitat suitability of ungulates at northern distribution ranges (Rivrud et al. 2019). However, whether climate changes will reduce or increase snow cover depends on whether temperature or precipitation is limiting snow accumulation (resp. coastal vs. continental regions, Brown and Mote 2009). Thus, the effects of climate change on migratory populations and their harvest can be opposite in different regions.

Moreover, an increase in weather extremes and variation is expected due to climate change (Coumou and Rahmstorf 2012), which can lead to years when migration becomes a particularly important tactic for a population to deal with such variation spatially (e.g., Stien et al. 2010). In general, the increased variability in living conditions due to climate change will require an increased plasticity in migratory behavior (Moore 2011).

Migration mortality and behavioral avoidance of barriers

Our approach and results enhance the understanding of how the reported negative effects of migration barriers on population dynamics may arise (reviewed in Bolger et al. 2008). In our analysis, we considered effects from both migration cost and behavioral avoidance of barriers on animal migration (Shepard et al. 2008). The effects of an increased cost to migration or behavioral avoidance of a barrier were similar as both led to a reduction in the migration probability. A migratory population becomes fully resident through either complete behavioral avoidance of the barrier (i.e., b = 1), or when the cost associated to migration outweighs the benefits (i.e., \(c > r_{H}[1-N_{H}] - r_{L}[1-N_{L}] - ({h_{H}^{S}} - {h_{L}^{S}})\)). In our model, we assumed that animals responded to the increased cost of migration. It is, however, possible that certain species are not responding adequately to the increased mortality associated with a barrier. In such cases, a barrier could have even larger effects on the dynamics of a population.

The increase in population size we found for small behavioral barriers is a consequence of the difference between ideal-free and average-fitness-maximizing migration (Fig. 2); this difference results from the density regulation. With the Ricker-type density regulation in our model, no positive effect would occur from a barrier on the ideal-free migration probability when rH ≤ rL, the ideal-free migration probability would then equal or be lower than the average-fitness-maximizing one.

In the most strongly seasonal environments (i.e., where each seasonal range cannot sustain a population year-round: 1 − ϕ ≤ 0), we expect that the introduction of a strong migration costs and/or behavioral avoidance of barrier will lead to a population collapse. However, even when one or both of the seasonal ranges can sustain a resident population a migration barrier often results in a decreased population growth (e.g., Skogland and Mølmen1980).

Migration, harvest, and management

Migratory animals create a spatial dependency between areas, which results in spatial subsidies (Semmens et al. 2011) when the societal costs and benefits from a species are season dependent (Skonhoft and Olaussen 2005; Nilsen et al. 2009). Such spatial subsidies are often a challenge for managers, when administrative boundaries are crossed (Boyce 1991). In environments with high year-round carrying capacity, barriers can locally increase populations and harvest yields resulting in a more equitable distribution of costs and benefits thereby reducing societal conflicts. However, such positive barrier effects are likely countered by the decrease in overall population size and harvest yield, when the year-round carrying capacity of the ranges is low (i.e., none or only small resident populations are supported).

In addition, changing harvest and mortality can shift a population from being resident to migratory and vice versa. Empirical studies have documented increased residency following increased mortality during migration due to non-human predators (Middleton et al. 2013; Hebblewhite and Merrill 2011). Although several studies have demonstrated the effects of human hunting on the timing of migration (Bechet et al. 2003; Rivrud et al. 2016), we are not aware of studies showing that human harvest has altered the migratory behavior of a population altogether.

The timing of harvest is crucial for the spatial distribution of the yield from a migratory population. Summer harvest will lead to a high yield in the highland range and results in the maximum sustainable yield for partially migratory populations. On the other hand, winter harvest provides the benefits from harvest to the lowland range. For example, the harvest season has been extended to allow for winter harvest of migratory moose to redistribute benefits among landowners in some municipalities of Norway (Skonhoft and Olaussen 2005). Also in red deer, a change in the harvesting season resulted in a marked shift in the spatial distribution of the harvest (Loe et al. 2016). Thus, the timing of the harvest has important consequences for the spatial subsidies that areas provide to each other (Semmens et al. 2011).