Consequences of the genetic threshold model for observing partial migration under climate change scenarios

Abstract Migration is a widespread phenomenon across the animal kingdom as a response to seasonality in environmental conditions. Partially migratory populations are populations that consist of both migratory and residential individuals. Such populations are very common, yet their stability has long been debated. The inheritance of migratory activity is currently best described by the threshold model of quantitative genetics. The inclusion of such a genetic threshold model for migratory behavior leads to a stable zone in time and space of partially migratory populations under a wide range of demographic parameter values, when assuming stable environmental conditions and unlimited genetic diversity. Migratory species are expected to be particularly sensitive to global warming, as arrival at the breeding grounds might be increasingly mistimed as a result of the uncoupling of long‐used cues and actual environmental conditions, with decreasing reproduction as a consequence. Here, we investigate the consequences for migratory behavior and the stability of partially migratory populations under five climate change scenarios and the assumption of a genetic threshold value for migratory behavior in an individual‐based model. The results show a spatially and temporally stable zone of partially migratory populations after different lengths of time in all scenarios. In the scenarios in which the species expands its range from a particular set of starting populations, the genetic diversity and location at initialization determine the species’ colonization speed across the zone of partial migration and therefore across the entire landscape. Abruptly changing environmental conditions after model initialization never caused a qualitative change in phenotype distributions, or complete extinction. This suggests that climate change‐induced shifts in species’ ranges as well as changes in survival probabilities and reproductive success can be met with flexibility in migratory behavior at the species level, which will reduce the risk of extinction.

The inclusion of such a genetic threshold model for migratory behavior leads to a stable zone in time and space of partially migratory populations under a wide range of demographic parameter values, when assuming stable environmental conditions and unlimited genetic diversity. Migratory species are expected to be particularly sensitive to global warming, as arrival at the breeding grounds might be increasingly mistimed as a result of the uncoupling of long-used cues and actual environmental conditions, with decreasing reproduction as a consequence. Here, we investigate the consequences for migratory behavior and the stability of partially migratory populations under five climate change scenarios and the assumption of a genetic threshold value for migratory behavior in an individual-based model. The results show a spatially and temporally stable zone of partially migratory populations after different lengths of time in all scenarios. In the scenarios in which the species expands its range from a particular set of starting populations, the genetic diversity and location at initialization determine the species' colonization speed across the zone of partial migration and therefore across the entire landscape. Abruptly changing environmental conditions after model initialization never caused a qualitative change in phenotype distributions, or complete extinction. This suggests that climate change-induced shifts in species' ranges as well as changes in survival probabilities and reproductive success can be met with flexibility in migratory behavior at the species level, which will reduce the risk of extinction.

K E Y W O R D S
evolution, genetic diversity, individual-based model, passerines, selection landscape, spatially explicit

| INTRODUCTION
Migration is a widespread phenomenon across the animal kingdom as a response to seasonality in environmental conditions, with a wide variety in migratory strategies across and within species (Dingle & Drake, 2007). In quite some species of birds, the young do not migrate in groups or families, clearly demonstrating that migratory behavior in space and time is endogenous (Liedvogel, Åkesson, & Bensch, 2011).
Many of these migratory species actually show a geographical cline in migratoriness, meaning that the species' area of distribution ranges from populations with migratory individuals to populations with residential individuals, following an environmental gradient (Newton, 2008;Sahashi & Morita, 2013). Partially migratory populations, in which both phenotypes co-exist (Chapman, Brönmark, Nilsson, & Hansson, 2011b), can in such species occur between these two kinds of populations. The stability of partially migratory populations has long been debated (Chapman, Brönmark, Nilsson, & Hansson, 2011a).
Clearly, under such conditions, all individuals carry the genetic migration program, yet it is differentially expressed between the migratory and residential phenotypes (Franchini et al., 2017). As reasons for the differential expression of the genes related to migratory behavior, most empirical studies point toward fluctuating differences in survival or fecundity between phenotypes (Rolandsen et al., 2017) and an important role for density dependence in the maintenance of the two phenotypes at the same location (Holt & Fryxell, 2011;Kaitala, Kaitala, & Lundberg, 1993;Kokko, 2011;Lundberg, 2013;Shaw & Levin, 2011).
In passerines, studies looking into the heritability of migratory behavior have concluded that the inheritance of migratory activity is currently best described by the threshold model of quantitative genetics (Pulido, 2007(Pulido, , 2011Pulido & Berthold, 2003Pulido, Berthold, & Van Noordwijk, 1996;Van Noordwijk et al., 2006). In a recent study by Cobben and Van Noordwijk (2016), a parsimonious model was used to investigate whether the inclusion of such a genetic threshold model for migratory behavior can explain the existence of partially migratory populations. They conclude that assuming a threshold model leads to a spatially and temporally stable zone of partially migratory populations under a wide range of demographic parameter values. While the specific location and width of this zone of partial migration vary with dispersal distances and the strength of density dependence, structural differences between phenotypes and density dependence are no prerequisites for obtaining such a zone. Cobben and Van Noordwijk (2016) investigate the equilibrium selection pressures on the threshold value across the species' range under model initializations with the number of individuals and the level of genetic diversity at carrying capacity.
Global warming has profound effects on geographical species distributions, showing range shifts for many taxa as a result of changing environmental conditions (Chen, Hill, Ohlemueller, Roy, & Thomas, 2011;Parmesan & Yohe, 2003). Particularly, for migratory species, it is expected that arrival at the breeding grounds might be increasingly mistimed under global warming as a result of the uncoupling of long-used cues and actual environmental conditions, with decreasing reproduction as a consequence (Both & Visser, 2001;Lameris et al., 2017;Miller-Rushing, Høye, Inouye, & Post, 2010;Møller, Rubolini, & Lehikoinen, 2008). In this study, we investigate the consequences for migratory behavior and the stability of the zone of partial migration under the assumption of a threshold model for migratory behavior, under a wide range of climate change-related scenarios. For this, we use the model of Cobben and Van Noordwijk (2016) with the purpose to provide testable hypotheses for empirically investigating the support for this threshold model in the field. In addition, we relax their initialization assumptions regarding the level of genetic diversity to investigate the consequences of a reduced availability of genetic variation for the model outcomes. Specifically, we look into these five scenarios: (i) initialization with zero genetic diversity, both with a universally beneficial threshold value and a locally beneficial threshold value; (ii) initialization of the species at only the extreme ends of the landscape, so with migrant populations or resident populations only, to investigate the speed of range expansion after initialization; (iii) improving winter survival probability from the viewpoint that an increased winter temperature will relax stressful environmental conditions; (iv) deteriorating winter survival probability from the viewpoint that species interactions can mismatch and predictability of the conditions is lowered; (v) deteriorating reproductive success for migratory individuals from the viewpoint that these will suffer most from unpredictability of the environment, resulting in mistimed arrival at the breeding grounds. We are interested in the relationship between the level of genetic diversity and flexibility of migratory behavior at the species level, and specifically explore the speed of changes in the frequency of migratory phenotypes across generations following changes in environmental conditions, and the speed of colonization across the partial migration zone under initializations with different levels of genetic variation.

| METHODS
We use the same model as in Cobben and Van Noordwijk (2016), a spatially explicit individual-based metapopulation model of a haploid species with discrete generations, and the rules for determining the population dynamics are inspired by passerine ecology. The species has a single trait that determines the threshold for migration, which is allowed to evolve.

| Landscape
The simulated landscape consists of 100 columns (x-dimension) of 25 breeding patches each (y-dimension), all with carrying capacity K. We assume bouncing borders in both x and y directions.
Hence, an individual cannot leave the landscape by dispersal, but migratory individuals are assumed to not spend their winters in the landscape that signifies only the species' breeding area. There are two gradients in the x-dimension, of increasing winter survival and decreasing reproduction, simulating average conditions on a trajectory from the pole to the equator. The experienced nominal average winter survival for residents is 0.01 times the x-location, so changes linearly from 0.01 at x = 1 to 1 at x = 100 (Equation 1).
Reproduction is corrected for position along the x-dimension through a reproduction factor that decreases linearly from 1 at x = 1 to 0.8 at x = 100 (Equation 2), incorporating the process that these species prefer more polar breeding areas where a highquality food peak can sustain larger broods. Both winter survival and reproduction are reduced by density-dependent factors (see below for further details).

| Winter survival and population dynamics
Local populations are composed of haploid individuals, each of which is characterized by a single trait, its migration threshold. If the local expected winter survival is lower than the genetically determined threshold T, the individual migrates (Able & Belthoff, 1998). This expectation includes all processes determining winter survival, so also population density as explained below, and thus models the process that an individual makes an estimate, ahead of time, of whether local winter conditions will be good enough to reach its own endogenous threshold in survival probability, based on environmental cues. A migrant is subjected to a constant probability of winter (thus migration) survival s m . In contrast, a resident can survive the winter locally with probability s r that is determined by its x-location. Resident winter survival is in addition controlled by the local resident density, so in all determined according to the following equation: where x is the x-coordinate of the patch [1,100]; Nr x,y,t is the number of resident individuals in patch x, y at time t; K is the carrying capacity; c_dens s governs the strength of density dependence, with no density dependence at c_dens s = 1, and s r = 0 under c_dens s = 0.6 and density = 2.5 (Fig. S1a in Appendix A and Table 1). In our simulations, the winter survival of migrants s m was 0.5. In our model, the migrant survival is not determined by local winter population densities, as individuals are known to move further south if the local carrying capacity is reached. In addition, migrant survival is thought to be mostly determined by mortality during migration.
If an individual survives the winter season, either as a resident or a migrant, it can disperse from its native patch with probability d to a neighboring cell with a maximum dispersal distance of δ cells. It then produces a number of offspring. This number is randomly drawn from a discrete uniform distribution with a minimum value of zero and a maximum of six, and then corrected for density and location, leading to the following equation: where No is drawn from U {0,6}; N x,y,t is the number of individuals in patch x,y at time t (i.e., including both the migrants and residents); c_loc r governs the sensitivity of reproduction with location x, with no dependence for c_loc r = 1, and maximum dependence of c_loc r = 0; c_dens s governs the strength of density dependence, with no dependence at c_dens s = 1; c tot is then the total reproduction correction factor (Fig. S1b in Appendix A). After reproduction, the individual dies.
Although carrying capacity K is equal for both Equations 1 and 2, the population sizes and thus the densities change over the seasons.
In addition, the parameters governing the strength of density dependence are different for both equations (c_dens s and c_dens r , Table 1), which results in different effects of density on resident winter survival and reproduction, respectively.

| Genetics
The individuals have one genetic trait, the migration threshold, which is located at a single gene. The threshold value T is a continuous number in the range [0,1]. The individuals are haploid and reproduce asexually, thus the offspring inherit the threshold values from their mothers. To allow for some variation, we apply a probability of mutation m. When such a mutation occurs, the threshold value is changed into a new randomly chosen value [0,1].

| Simulation experiments
With our simulation experiments, we want to test how sensitive the model outcomes in Cobben and Van Noordwijk (2016) are to different model starting conditions (initializations) and abrupt changes in environmental conditions during the model run, inspired by  Table S1).
Changing survival and reproduction conditions were implemented 300 years after the start of the model run, when demographic equilibrium was reached. The only exceptions to this are the scenarios implemented with populations initialized across the landscape in combination with zero genetic diversity, in order to prevent genetic diversity from increasing before changing conditions were applied.
For each scenario, we performed 100 replicate simulations.

| Analysis
The individual phenotypes, that is, migrancy versus residency, were documented in time and space throughout the simulations and summed per x-location. The migration threshold values were averaged per x-location. We defined the polar border of the partial migration zone as the smallest x-location where the total number of residents was at least 1% of the total number of individuals. The equatorial border is equivalently defined as the largest x-location where at least 1% of the total number of individuals was a migrant. The difference between both borders is the width of the partial migration zone, while the location of the zone was defined by its middle x-location.

| RESULTS
All runs of all scenarios converge to an equilibrium situation with fully migratory populations at the lower x-values, fully residential populations at the higher x-values, and a spatially and temporally stable zone of partially migratory populations between these (see Figure 1 for a typical example).

| Range expansion scenarios
In the range expansion scenarios, the genetic diversity and location at initialization determine the species' colonization speed across the zone of partial migration and therefore across the entire landscape T-values (and thus migratory phenotype) can later replace residents at x-values larger than 70 (Fig. S3).

| DISCUSSION
Cobben and Van Noordwijk (2016)  ing areas, and also population shifts either toward more migratory or more sedentary and even residential behavior (Able & Belthoff, 1998;Berthold, Helbig, Mohr, & Querner, 1992;Fiedler, 2003;Knudsen et al., 2011;Partecke & Gwinner, 2007;Pulido & Berthold, 2010). This flexibility is often attributed to plasticity as opposed to evolutionary change in migratory behavior (Knudsen et al., 2011). However, the finding that the inheritance of migratory activity is best described by the threshold model of quantitative genetics (Pulido, 2011) reconciles these viewpoints regarding the issue of migratory versus residential behavior for species where migration is controlled endogenously. In this study, we investigate the evolutionary change in local threshold values that determine migratory behavior in combination with local survival conditions, which is thus a plastic trait. In migratory fish, this was studied by Sahashi and Morita (2013), who show that selection on the threshold value (in their case a threshold in size at maturity), resulted in a decreasing threshold value with increasing distance from the sea, leading to an increased number of residential individuals as migration costs increase. Also here, migratory behavior is plastic, as a fish of a certain size at maturity is migratory under downstream conditions, and residential in upstream populations, which is comparable to our trait, where a specific T-value will give different migratory phenotypes under different environmental conditions.
Despite these observed rapid changes in migratory behavior in response to changing environmental conditions, it is widely acknowledged that migratory species in general seem to suffer more from population declines than species that do not migrate, with climate change as one of its four drivers (Wilcove & Wikelski, 2008).
It is predicted that communities of migratory bird species in Europe will be altered through changes in migratory behavior under current projected climatic changes, rather than through reassembly of bird communities, the latter of which would lead to lower proportions of migratory species (Schaefer, Jetz, & Böhning-Gaese, 2008 (Matsumura, 1996;Saulich & Musolin, 1996;Söderlind & Nylin, 2011). Some of these traits, for example, diapause, also fulfill the condition that local densities for one phenotype are lower due to the existence and inherent absence of the other phenotype. Although density dependence is not essential for the maintenance of populations consisting of both phenotypes, it does substantially increase the area where mixed populations are observed (Cobben & Van Noordwijk, 2016). A major difference with the life history studied here is that in insects the reproduction is much higher. From this, one can expect an increase in the number of mutations which will likely cause an even wider zone of mixed populations. The modest maximum dispersal distance employed in this study is fairly agreeable to insect life history, in all implying that our model is also applicable to insect life histories and well-known threshold traits in these.
Our model is artificial in a number of ways: We assume discrete generations and haploid inheritance. The results thus mainly provide a view of the selection landscape, in interaction with the limited mutation and dispersal rates. The result that shifts in the selection pressure lead to quick responses while maintaining stable mixed populations were obtained under a wide range of scenarios and parameter values. Another simplification of reality was made when assuming that individuals can perfectly predict what their local survival probability is, that is, they "know perfectly" where they are. This aspect has been thoroughly tested by Cobben and Van Noordwijk (2016), who show that decreased predictability of survival conditions has a small effect on the mean and standard deviation of the location of the partial migration zone, but does not cause a qualitatively different result.
In this study, we investigated the flexibility of migratory behavior at the landscape and species level, and stability of partially migratory populations under a set of different initialization conditions and F I G U R E 3 Scenarios with changing environmental conditions. The average x-location of the zone of partial migration over 100 replicate runs, with the error bars expressing standard deviations, in time after the change in survival probability was abruptly implemented in the model simulation for eight different scenarios (a, b). Change was applied 300 generations after model implementation except for the "compl-zero" scenario. Decreased reproduction for migrants (c) was simulated in two different scenarios after demographic equilibrium was reached at generation 300. "Compl" = initialized with individuals 1 ≤ x ≤ 100, "migr" = 1 ≤ x ≤ 10, "res" = 90 ≤ x ≤ 100, "full" = initialized with 0 ≤ T ≤ 1, "zero" = T = 0.5, "zerofit" = T = 0.9 (migrants) and T = 0.1 (residents). a: resident winter survival increase of 20%. b: resident winter survival decreased of 20%. c: migrant reproductive success decrease of 20% climate change-related scenarios. The results indicate that assuming a genetic threshold model for migratory behavior always leads to a spatial zone consisting of partially migratory populations and supports rapid changes of local observed migratory phenotypes, which is in agreement with empirical data.

ACKNOWLEDGMENTS
MC was funded by the Open Program of the Netherlands Organization of Scientific Research (NWO). We thank Bart Nolet and Thomas Lameris for fruitful discussions.