Migration pulsedness alters patterns of allele fixation and local adaptation in a mainland-island model

Abstract

Geneflow across populations is a critical determinant of population genetic structure, divergence, and local adaptation. While evolutionary theory typically envisions geneflow as a continuous connection among populations, many processes make it fluctuating and intermittent. We analyze a mainland-island model where migration occurs as recurrent “pulses.” We derive mathematical predictions regarding how the level of migration pulsedness affects the effective migration rate, for neutral and selected mainland alleles. We find that migration pulsedness can either decrease or increase geneflow, depending on the selection regime. Pulsedness increases geneflow for sufficiently (counter)selected alleles (⁠$s<s_1$⁠), but reduces it otherwise. We provide a mathematical approximation of the threshold selection strength $s_1$⁠, which is verified in stochastic simulations. Migration pulsedness thus affects the fixation rate at different loci in opposite directions, in a way that cannot be described as a change in effective population size. We show that migration pulsedness would generally reduce the level of local adaptation and introduce an additional genetic load: the “pulsedness load.” This is detrimental to the adaptation and persistence of small peripheral populations, with implications in management and conservation. These results indicate temporal variability in migration patterns may be an important, yet understudied, controller of geneflow and local adaptation.

temporally variable migration, effective migration rate, geneflow, migration load, selection

Introduction

Geneflow between populations, as a major determinant of evolutionary dynamics, has received considerable interest since almost a century. Depending on its intensity and interactions with other evolutionary forces, geneflow has a range of contrasting effects (e.g., Bürger, 2014; Felsenstein, 1976; Lenormand, 2002; Tigano & Friesen, 2016). In a focal population, it can enhance genetic diversity, prevent inbreeding, or on the contrary hamper local adaptation (Bürger & Akerman, 2011; Garant et al., 2007; Gomulkiewicz et al., 1999). Across populations, it controls the spatial spread of novel mutations, the maintenance of polymorphisms, the level of population divergence and, eventually, the possibility of speciation (e.g., Feder et al., 2019; Johnson et al., 2000; Mailund et al., 2012; Maynard-Smith, 1966; Rousset, 2013; Yeaman & Otto, 2011).

The consequences of geneflow have been studied in a range of spatial configurations, such as two interconnected populations (Maynard-Smith, 1966; Yeaman & Otto, 2011), mainland-island systems (Bürger & Akerman, 2011; Johnson et al., 2000) or metapopulations (Feder et al., 2019; Rousset, 2013; Slatkin, 1981). However, theoretical studies usually consider the process of migration (here used in the sense of dispersal) as constant and only rarely consider temporal variability in the flow of propagules. This is unfortunate, as migration is often governed by time-fluctuating and potentially highly unsteady phenomena that would make it temporally variable (Peniston et al., 2019).

Causes of temporal variability in migration rate are almost endless. Among the most frequently cited phenomena, a first category relates to environmental variation. Geographical barriers can change depending on land bridges (Keyse et al., 2018; Morris-Pocock et al., 2016), sea levels fluctuations (Hewitt, 2000), or habitat fragmentation (Peacock & Smith, 1997). Dispersal can also be affected by variation in oceanic and atmospheric currents (Renner, 2004; White et al., 2010; Smith et al., 2018; Benestan et al., 2021; see also Catalano et al., 2020 for a recent attempt to quantify dispersal variability). Last, but not least, dispersal is affected by extreme meteorological or climatic events such as floods and storms (Boedeltje et al., 2004; Reed et al., 1988) that can cause rafting events (Carlton et al., 2017; Garden et al., 2014). These extreme phenomena are bound to become more prevalent with climate change (Masson-Delmotte et al., 2018). A second category includes variation caused by species life-history traits, such as mast seeding, ballooning (Bishop, 1990) or various forms of group or clump dispersal (Soubeyrand et al., 2015). And finally, a third category encompasses variation related to dispersal by animal vectors (Martin & Turner, 2018; Yamazaki et al., 2016) or human activities (through ballast waters for instance, see Carlton & Cohen, 2003). All three categories are known to result in variable dispersal rates, in the form of random fluctuations (Yamazaki et al., 2016), intermittent flows (Hewitt, 2000; Keyse et al., 2018), or episodic bursts of migration (Carlton & Cohen, 2003; Morris-Pocock et al., 2016; Peacock & Smith, 1997; Reed et al., 1988; Reiners & Driese, 2004). The latter form is particularly common and is often referred to as “pulsed migration” (Bobadilla & Santelices, 2005; Boedeltje et al., 2004; Folinsbee & Brooks, 2007; Martin & Turner, 2018; Smith et al., 2018).

The relatively few existing theoretical studies suggest that migration variability can have important implications for population genetics and adaptation. Nagylaki (1979), Latter and Sved (1981), and Whitlock (1992) have all used discrete-time island models and considered neutral genetic variation only (infinite alleles model). All three concur in predicting that temporal variability in migration rates should decrease effective geneflow, and thus increase differentiation among populations (see also Rousset, 2013). In a two-population model, Yamaguchi and Iwasa (2013), and following papers, studied the fixation of incompatible mutations and the progress to allopatric speciation, and found that speciation occurred faster with variable migration. Some studies have investigated non-neutral cases in spatially heterogeneous environments, in particular in the context of source-sink population systems. Gaggiotti and Smouse (1996) found that spacing out migration events, while keeping the mean number of migrants constant, decreases the level of genetic variation in the sink population. Rice and Papadopoulos (2009) studied a mainland-island model and suggested that neglecting migration stochasticity would generally lead to overestimating the impacts of migration on adaptation. More recently, Peniston et al. (2019) extended the results of Gaggiotti and Smouse (1996), investigating the impact of temporally pulsed migration on the level of local adaptation in the sink population. They found that spacing out migration events (with a constant mean number of migrants) can either hamper or facilitate adaptation in a harsh sink environment, depending on genetic scenarios. In a different context, Matias et al. (2013) studied species diversity in a metacommunity model, and found that with randomly fluctuating migration rates, larger mean dispersal values were needed to produce the same local species richness. Overall, these studies seemingly converge on a negative impact of dispersal variability on effective migration rate. However, they remain few and limited in their scope. Most of them underscore the need to pay more attention to the consequences of migration variability (Peniston et al., 2019).

The present study aims at providing a more comprehensive appraisal of the consequences of temporal variability in migration, considering both neutral, beneficial, and deleterious alleles, under various forms of selection and levels of dominance. We will consider a mainland-island system where a small island population receives migrants from a large mainland population (Bürger, 2014; Felsenstein, 1976; Wright, 1931). In the island population, we will model the dynamics of allele frequencies in continuous time, at one locus, in diploid individuals. We will focus on the case of “pulsed” migration patterns, i.e., geneflow from the mainland occurs as bursts of migration of variable size and frequency (e.g., Peniston et al., 2019; Rice & Papadopoulos, 2009). We explore the range of migration “pulsedness” from continuous migration (independent migration of individuals) to very pulsed migration (rare migration of groups of individuals).

Mathematical predictions for how the level of migration pulsedness should affect the effective migration rate are derived under a low-migration limit. These predictions are validated in stochastic simulations with large migration rates and logistic population dynamics on the island. We show that for neutral alleles, migration pulsedness reduces the effective migration rate, thus decreasing the rate of allele fixation, which is in accordance with earlier theoretical results. However, under general scenarios of selection, this conclusion can change quantitatively, and even reverse. We find that if the selection coefficient falls below some negative threshold value $s_1$⁠, for which we provide an analytical expression, the qualitative impact of migration pulsedness switches to positive, so that the effective migration rate increases for those deleterious alleles. Moreover, for sufficiently recessive deleterious alleles, migration pulsedness may have a non-monotonous effect on the effective migration rate. The interplay of selection, dominance and drift is thus found to play an essential role in determining the impact of migration variability on effective migration rate.

Our results show that the effect of migration pulsedness is not uniform across loci subject to different selection regimes and dominance levels. Over an entire genome, migration pulsedness effectively homogenizes the fixation rates across loci, in a way that cannot be simply described as a change in effective population size. This has consequences for the dynamics of mean fitness and local adaptation, by creating an additional migration load that we call the “pulsedness load.” Overall, our results highlight migration variability as an important property to be taken into account when studying geneflow and local adaptation. Specifically, pulsed migration patterns could overall be detrimental to the local adaptation and persistence of small peripheral populations.

Methods

Mainland-island model

A large mainland population is connected to a small island population of finite size $N$⁠, and migrant individuals flow from the former into the latter. We consider arbitrary loci at which some allele $A$ is fixed in the mainland whereas a different allele $a$ was initially fixed in the island population (Figure 1). Individuals are diploid (there are $2N$ genes at any locus) and mate randomly within the island population. The three diploid genotypes $AA$⁠, $Aa$⁠, and $aa$ have fitness values $1+s$⁠, $1+hs$⁠, and $1$⁠, with $h$ the degree of dominance and $s$ the selection coefficient, as usually defined (Thurman & Barrett, 2016). Results directly extend to haploids if one considers genic selection (⁠$h=0.5$⁠) and takes the number of genes to be $N$⁠. For simplicity, mutation is neglected. The evolution of allele frequencies in the island population is modeled in continuous time (overlapping generations, Moran, 1962).

Figure 1.

Number of migrants as a function of time in the continuous case (⁠$n=1$⁠, independent migration of individuals) and in a pulsed migration case (⁠$n=15$⁠). The left panel averages the number of migrants over one generation. The right panel is a close-up over five generations, showing the arrival of individual migrants at each migration event. Illustrative parameters: $m=10$⁠, $T=1$⁠.

Modeling of pulsed migration

Migration occurs as discrete migration events (pulses of migrants) corresponding to the arrival of $n$ individuals into the island population (see Yamaguchi & Iwasa, 2013, for a similar description of migration). The overall intensity of migration is controlled by the migration rate $m$⁠, the mean number of migrants per unit of time. To keep this mean number of migrants constant, we impose that the more individuals arrive per migration event (the larger $n$⁠), the less frequent migration events are (Peniston et al., 2019). Specifically, the rate at which migration events occur is taken to be $m/n$⁠. If $n=1$⁠, individuals migrate independently, corresponding to the classical continuous form of migration. When $n$ exceeds one, migrants arrive more sporadically, but in larger packs. We will thus increase $n$ to make migration more and more pulsed, and evaluate how this affects predictions compared to continuous migration (⁠$n=1$⁠).

Parameter $n$ controls the degree of migration pulsedness. More precisely, the number of migration events that occur in the time interval corresponding to one generation (⁠$T$⁠) follows a Poisson distribution with mean $mT/n$⁠. Therefore, $mT/n$ migration events occur on average, with a variance equal to $mT/n$⁠. It follows that the variance in the number of migrants (⁠$n\times mT/n$⁠) per generation is $n^{2}mT/n=nmT$⁠. Parameter $n$ thus quantifies the degree of overdispersion (variance inflation) in the number of migrants per generation, relative to the continuous case. Migration events can range from very frequent and of small intensity, with minimal temporal variance in the number of migrants per generation (low $n$⁠), to more infrequent and intense, with large temporal variance (large $n$⁠; Figure 1).

Metrics of interest

In such a mainland-island model, fixation of mainland alleles in the island eventually occurs with probability one, but we are interested in how long fixation takes, depending on the pulsedness level $n$⁠. In other words, we seek to understand whether migration pulsedness, all else equal, is favorable or unfavorable to geneflow. We want to know this for mainland alleles with different selective values (⁠$h$ and $s$⁠), and across different population sizes (⁠$N$⁠). To compare pulsed migration with continuous migration, we will compute the effective migration rate $m_e$⁠, defined as the migration rate that would be required to produce the same time to allele fixation, all else equal, under continuous migration (i.e., with $n=1$⁠). This definition is a variant of others (Kobayashi et al., 2008; Rice & Papadopoulos, 2009; Wang & Whitlock, 2003), adapted to the study of migration pulsedness. The value of $m_e$ is by definition $m$ when $n=1$⁠, and can deviate from $m$ for larger $n$⁠: if $m_e>m$⁠, migration pulsedness promotes geneflow, whereas if $m_e<m$⁠, it reduces geneflow. Following Kobayashi et al. (2008), we will further compute the geneflow factor, $m_e/m$⁠: if the geneflow factor is greater than one, migration mulsedness promotes geneflow, whereas if it is less than one, it reduces geneflow.

Note that all comparisons are done for a given locus with a given selective regime, not across different selective regimes. Obviously, beneficial alleles will on average go to fixation before neutral or maladaptive ones, for any kind of migration. Hence our effective migration rate and geneflow factor are not “relative to a neutral allele,” as sometimes used, but rather “relative to continuous migration, for an identical allele” (Kobayashi et al., 2008).

Stochastic simulations

We simulated the island population as a continuous time birth-death-immigration process, using a Gillespie algorithm. The population dynamics of the island population followed a stochastic logistic model with carrying capacity $K$ (Goel & Richter-Dyn, 1974). Each diploid individual has some basal death rate $d$ and birth rate $b=d$⁠, and density dependence acts through the death rate: total death rate is $d\frac{N^2}{K}$ while total birth rate is $bN$⁠. Selection occurs at reproduction (fertility selection): reproducing individuals are picked at random with odds proportional to their selective values. The selective values of genotypes otherwise do not affect the overall population dynamics. For computational efficiency, we used an optimized Gillespie algorithm (see Supplementary Information; Section 1 for a full description of the simulation algorithm).

The island is initially fixed for allele $a$ and at carrying capacity (⁠$N=K$⁠). At the beginning of a simulation, it starts receiving $AA$ migrants from the mainland, under the stochastic migration process described above (Figure 1). The simulation ends when allele $A$ gets fixed in the island. We conducted 10,000 replicates per parameter set to capture the stochasticity in migration times, population size fluctuations, and genetic drift. Effective migration rates were determined from the observed time it takes for the mainland allele to get fixed in the island, using a pre-computed abacus (Supplementary Information; Section 1.3).

Mathematical approximation

To allow analytical progress, we make the standard assumption that population regulation is fast, so that population size is effectively constant at $N$ in the island. However, we cannot employ the equally standard method of diffusion approximations. Indeed, the latter would require $n$ to be smallish relative to $N$⁠, de facto preventing the description of pulsed migration patterns (see Kimura, 1962; Yamaguchi & Iwasa, 2013). We rather use an alternative approximation, assuming that migration events are rare enough, so that genetic drift and selection result in allele fixation before the next migration event occurs.

With these assumptions, at each migration event, $n$ homozygous $AA$ migrants arrive in the island population containing either $N$$aa$ individuals (if fixation did not occur yet) or $N$$AA$ individuals (if it did). Following every migration event, population regulation quickly takes population size back to $N$⁠. Until fixation of the mainland allele, each migration event brings the latter to initial frequency $f=\frac{2n}{2n+2N}=\frac{n}{N+n}$⁠, and the outcome of the particular migration event is either fixation (success), or disappearance (failure) (see also Yamaguchi & Iwasa, 2013). The probability of success can be well approximated by the fixation probability of an allele in initial frequency $f$ in a closed population of size $N_e$⁠, $u(f,N_e,s,h)$⁠. Expressions of $u$ are are well known (Kimura, 1962; Whitlock, 2003). For simplicity, we will typically consider that $N_e$ is equal to $N$ in computations, even though $N_e$ could often be smaller than population size, including in our simulations (owing to demographic fluctuations). To simplify notations, we’ll write interchangeably $u\left(f\right)$ or $u(f,N_e,s,h)$⁠.

From this, we can derive the probability that the mainland allele has not yet swamped the island after a given time, and the mean time before it does, from an exponential distribution with rate $-u\left(f\right)m/n$ (Supplementary Information; section 2.2).

Parameter values

The island population should be sufficiently isolated, otherwise it will immediately be swamped by the mainland and geneflow is a trivial issue. A famous rule-of-thumb is the so-called “one migrant per generation rule” (Mills & Allendorf, 1996, see also Blanquart et al., 2013), stating that migration is strong if more than one migrant enters the population every generation. Without loss of generality, we will set the basal death rate as $d=N$ in our simulations, so that the generation time (defined as the average time for $N$ deaths to occur, see Moran, 1962) is $T\approx 1$⁠. Since the mean number of migrants per generation is $mT$⁠, this means that, roughly speaking, $m$ should not exceed $1$⁠. In turn, our mathematical approximation requires the average time between migration events (⁠$n/m$⁠) to be much larger than the average time to fixation after one event. Using classical theory (Kimura & Ohta, 1969; Otto & Whitlock, 2013; Whitlock, 2003), we can show that this always holds if $m\ll \frac{1}{4N}$⁠. This is a sufficient condition: the approximation should hold well even for larger migration rates (Supplementary Information; section 2.1).

We will consider four scenarios regarding the intensity of migration in simulations: (a) very low migration (⁠$m=0.001$⁠); (b) low migration (⁠$m=0.1$⁠); (c) medium migration (⁠$m=1$⁠), and (d) strong migration (⁠$m=10$⁠). The first scenario will be used to check that simulations do converge to our mathematical approximation. The three other scenarios span the one migrant per generation rule by one order of magnitude in both directions. The last scenario (strong migration) is the one illustrated in Figure 1.

The level of migration pulsedness $n$ will be systematically varied between $1$ and $N$⁠, the latter being an extreme case such that migration pulses temporarily double the local population size and brings migrants to $50\mathrm{\%}$ frequency. Selection parameters will also be systematically varied in plausible ranges: $s$ between $-0.1$ and $0.1$ (Martin & Lenormand, 2015; Thurman & Barrett, 2016), and $h$ between $0$ and $1$⁠. The population size $N$ will be varied between $50$ and $200$ (Palstra & Ruzzante, 2008; Peniston et al., 2019). With such parameters, the average proportion of immigrant individuals per generation (⁠$m/(N+m)$⁠) varies between $0.0005$ and $0.2$⁠.

Results

Mathematical predictions

General criterion for the impact of migration pulsedness

Under our mathematical approximation, the rate of fixation of a mainland allele is the product of $m/n$ and the probability of success $u\left(\frac{n}{N+n}\right)$ of a particular migration event (see above), and the mean time to fixation is the inverse of this rate (Supplementary Information; section 2.3). By definition of the effective migration rate, we must have $m_{e}u\left(\frac{1}{N+1}\right)=\frac{m}{n}u\left(\frac{n}{N+n}\right)$⁠. It follows that in order to have a geneflow factor greater than one (⁠$m_e>m$⁠, i.e., migration pulsedness increases geneflow), we must have:

$$\begin{array}{l}u\left(\frac{n}{N+n}\right)>nu\left(\frac{1}{N+1}\right)\end{array}$$

(1)

This condition yields a simple graphical criterion for determining whether migration pulsedness promotes or decreases geneflow, based on the shape of the $u$ function corresponding to any particular selective regime (Supplementary Figure S3). It also provides an operational way to derive quantitative mathematical predictions for specific types of genetic variation, as we will proceed to do now.

Neutral variation

At a locus where the mainland allele is neutral, we can apply criterion (1) using the fixation probability of a neutral allele, which is well known to be its initial frequency:

$$\begin{array}{l}u\left(\frac{n}{N+n}\right)=\frac{n}{N+n}\end{array}$$

(2)

In that case, it is straightforward to see that criterion (1) is never met, for any $n$⁠. Note that in the above equation, $N$ is the actual census population size, not $N_e$⁠, even if they differ. We thus predict that migration pulsedness should invariably decrease the effective migration rate, i.e., decrease geneflow. The more pulsed migration is, the larger the effect is.

Selected variation

In the general case of selection with dominance, an expression for function $u$ is:

$$\begin{array}{l}u\left(f\right)=\frac{{\int }_0^{f}G\left(x\right)dx}{{\int }_0^{1}G\left(x\right)dx}\end{array}$$

(3)

with $G(x)=\exp [2N_{e}s((2h-1)x^2-2hx)]$ (Kimura, 1962).

Note that in the expression of $G\left(x\right)$⁠, the effective population size $N_e$ matters if different from the census population size $N$⁠.

This expression allows much greater flexibility in the shape of the function. Function $u\left(f\right)$ is well known to be either strictly concave, strictly convex, or convex-concave with an inflexion point, depending on $s$ and $h$⁠. A particular case is the additive case (frequency-independent selection in haploids, or codominance in diploids, i.e., $h=1/2$⁠), for which the explicit solution $(1-\exp (-2f{N}_{e}s))/(1-\exp (-2N_{e}s))$ is well known (Kimura, 1962; Maruyama, 1970; Whitlock, 2003).

Predictions based on numerical analyses for all $s$ and $h$ values are presented in Figure 2. Even though equation (3) is much less mathematically tractable, we can still make some analytical progress. Criterion (1) can be viewed as a comparison to linearity for function $u\left[f\left(n\right)\right]$ (Supplementary Information; section 2.4): if initially (for low levels of pulsedness, i.e., $n$ close to 1), the function is concave (has negative second derivative), then the criterion will not be met. If on the contrary it is convex, then the criterion will be met. We can therefore determine whether making migration pulsed increases or decreases geneflow by studying the curvature of $u\left[f\left(n\right)\right]$ at $n=1$⁠. Since $f$ is a concave function of $n$⁠, we know that $u\left[f\left(n\right)\right]$ will be more prone to being concave than $u\left(f\right)$⁠. This graphical interpretation of equation (1) is illustrated in Supplementary Figure S3.

Figure 2.

Mathematical predictions. For all combinations of selection coefficient $s$ and dominance $h$⁠, the plane is partitioned according to the impact of migration pulsedness on effective migration rate (geneflow). Predictions were computed numerically (from equations (1) and (3)). Triangles show the observed result in stochastic simulations under the very low migration scenario (⁠$m=0.001$⁠), to check consistency with mathematical analysis. The dashed curve is the mathematical expression for $s_1$ (equation (4)). Subpanels are cross-sections of the main panel at four levels of dominance (⁠$h=0.8,0.3,0.1,\mathrm{a}\mathrm{n}\mathrm{d}0$⁠), showing $n$ values that increase (white) or decrease (gray) the effective migration rate. Note the square-root y-scale. The blue curves show the $n$ values that maximize (solid) or minimize (dotted) geneflow. Other parameters: $N=N_e=100$⁠.

Differentiating $u\left[f\left(n\right)\right]$ twice at $n=1$⁠, using equation (3) and $f=n/(N+n)$⁠, we find after some calculations (Supplementary Information; section 2.4) that migration pulsedness increases the effective migration rate (promotes geneflow) if:

$$\begin{array}{l}s<s_1=-\frac{{(N+1)}^2}{2N_{e}N(hN+1-h)}\end{array}$$

(4)

The threshold selection coefficient $s_1$ is always negative. In the special case of additivity (⁠$h=1/2$⁠), if population size is large and no different from the effective population size, we obtain the approximation $s_1\approx -1/N$⁠. This becomes $s_1\approx -1/\left(2N\right)$ for fully dominant alleles, whereas for fully recessive alleles, we get the quite extreme value $s_1\approx -1/2$⁠. These results are in close agreement with the numerics and with stochastic simulations under a very low migration scenario (Figure 2).

We remark that the values of $n$ at which an initially positive effect of pulsedness would revert to negative are usually unrealistically large (thick dark-green curves in subpanels of Figure 2; see also Supplementary Information; section 2.7). Therefore, equation (4) in practice means that for sufficiently deleterious alleles (⁠$s<s_1$⁠), migration pulsedness increases the effective migration rate, instead of reducing it as it does for neutral alleles. It reduces geneflow, as for neutral alleles, for slightly deleterious alleles (⁠$s_1<s<0$⁠), and for all beneficial alleles (⁠$s>0$⁠).

The condition $s<s_1$ is sufficient, but there still is the possibility that function $u\left(n\right)$ is first convex, then switches to concave to such an extent that criterion (1) is met for some larger $n$ value, even though it was not for small values. In such circumstances, the consequences of migration pulsedness further depend on the value of $n$⁠: it would generally decrease the effective migration rate, except for some intermediate range of $n$ values, for which it would on the contrary increase it. This situation requires function $u\left(f\right)$ to have a pronounced sigmoidal shape, which is known to occur only for recessive deleterious alleles. This is illustrated in Supplementary Figure S3c.

Accordingly, as can be seen in Figure 2, for sufficiently recessive alleles, increasing $s$ above $s_1$ may not immediately turn the effects of migration pulsedness from positive to negative, but rather there can be an intermediate range $(s_1,s_2)$ of values for which migration pulsedness has mixed effects depending on its intensity. One can determine the level of dominance at which this bifurcation occurs (⁠$h_c$⁠), by remarking it is the point $(h,s_1\left(h\right))$ at which the second inflexion point of function $u\left(n\right)$ collides with the first one (located at $n=1$⁠). In other words, it is the point at which function $u\left(n\right)$ has a null third derivative at $n=1$ while verifying the condition $s=s_1$ (Supplementary Information; section 2.5). After some calculations, we get the following condition:

$$\begin{array}{l}h<h_c=\frac{N-1}{3N-1}\end{array}$$

(5)

Note that this criterion depends on the census population size only, not on the effective population size. If population size is large, we get the approximation $h\approx 1/3$⁠, in close agreement with the numerics (Figure 2).

We remark that, for $h<h_c$ and $s<s_2$⁠, the minimal level of pulsedness required to get an increased effective migration rate is typically quite large (⁠$40--100$⁠; thin light-green curves in the subpanels of Figure 2). Therefore, in practice, this situation is essentially similar to that of neutral alleles (migration pulsedness reduces the effective migration rate).

To summarize, migration pulsedness should increase the effective migration rate for sufficiently maladapted and not-too-recessive alleles (⁠$s<s_1$⁠), whereas it should reduce the effective migration rate for slightly maladapted alleles, neutral alleles and beneficial alleles.

Confronting predictions with stochastic simulations

The above mathematical predictions were well verified in stochastic simulations in which the overall rate of migration $m$ was allowed to take larger values (Figure 3). To facilitate comparisons across parameter sets, we follow the usual practice of plotting results in terms of $N_{e}s$ rather than $s$⁠, to reflect the effective strength of selection relative to drift (Lande, 1994; Wright, 1931).

Figure 3.

Results of stochastic simulations. Left: Effect of migration pulsedness on effective migration rate, for $n$ between 1 and $N$⁠, and for a range of $Ns$ values. The solid black lines represent the mathematical prediction, obtained numerically from equations (1) and (3). Results are provided for the tree scenarios of migration intensity (columns), and for three population sizes (rows). Right: Geneflow factor (⁠$m_e/m$⁠) as a function of $n$⁠, for different selection coefficients $s$ and for the three migration scenarios. The $s$ values used are also positioned in the left panels as vertical lines. Symbols are the values obtained in stochastic simulations, and solid lines are mathematical predictions, obtained as before. Other parameters: $h=1/2$⁠.

Qualitatively, mathematical predictions hold almost perfectly in the low and medium migration scenarios. In other words, migration pulsedness could either increase or decrease the effective migration rate, depending on the type of selective regime, in the direction predicted. In the strong migration scenario, predictions also held except that the space of selective regimes for which migration pulsedness increases the effective migration rate was shifted to lower $s$ values (i.e., $s_1$ gets more negative).

Results did not depend importantly on population size, provided one scales $n$ with $N$⁠. This is again consistent with mathematical predictions: the left-hand side of equation 1 is unchanged if $N_{e}s$ and $n/N$ (and thus $f$⁠) are both kept constant (see equation 3). The right-hand side changes very little: its relative variation with $N$ is of order less than $1/N^2$ (see Supplementary Information; section 2.8). Therefore, mathematical predictions are almost invariant for given values of $N_{e}s$ and $n/N$⁠.

In quantitative terms, the impact of migration pulsedness on the effective migration rate (the geneflow factor) was very close to the one predicted by our mathematical approximation in all cases, except for deleterious alleles under the strong migration scenario (Figure 3). In the latter case, the observed increase in geneflow could be much smaller than the one predicted. In the most extreme example, a predicted six-fold increase in effective migration rate turned out to be slightly more than a two-fold increase with $m=10$ (Figure 3). Overall, large migration rates (i.e., beyond the one-migrant-per-generation rule) shift predictions toward a decrease in effective migration rate: more negative $s$ values are required to obtain the same impact as for low migration rates. They otherwise have little qualitative impact. Remark that these quantitative differences are most important for large pulsedness levels (⁠$n$⁠); results for low $n$ values, and thus the value of $s_1$⁠, are in much closer agreement with mathematical predictions (Figure 3). Similar conclusions hold for other values of the dominance level $h$ (Supplementary Figure S5).

Genome-wide consequences: mean fitness and the pulsedness load

The above results have established that migration pulsedness should affect different loci in different ways, depending on their selective regime. Overall, it will promote geneflow (promote fixation) of deleterious or maladapted alleles, especially if not too-recessive, whereas it will reduce geneflow (slowdown fixation) of slightly deleterious, neutral and beneficial alleles, especially recessive ones.

As pulsedness promotes the fixation of maladapted alleles and opposes the fixation of beneficial alleles, its overall effect is to homogenize geneflow (fixation rates) across loci with different selective regimes in a genome. Such a homogeneizing effect could at first be regarded as similar to a decrease in effective population size: demographic fluctuations or other sources of variability that decrease $N_e$ would reduce the strength of selection and also homogenize fixation rates across different loci.

However, simply decreasing population size, in a model with continuous migration, cannot reproduce the patterns observed with migration pulsedness (Figure 4). Even if adjusting potentially unknown parameters (⁠$m$ or the time since onset of migration), in order to match exactly the fixation rate of neutral alleles, a reduction in $N_e$ can either fit the fixation patterns observed for maladapted alleles, or those of beneficial alleles, but never fits both (Figure 4). In other words, the change in $N_e$ that is required to mimic the consequences of migration pulsedness is not identical across loci differing in their selection regimes. Therefore, the consequences of migration pulsedness cannot be interpreted as some reduction in effective population size.

Figure 4.

Probability of fixation of different alleles varying in their selection coefficient, under pulsed migration (⁠$n=30$⁠; red curve), and under continuous migration with different effective population sizes (gray dashed curves). (A) $N=100$ and (B) $N=200$⁠. The time since the onset of migration (⁠$t$⁠) was set at an intermediate value so that about $50\mathrm{\%}$ of neutral alleles have gone to fixation in the pulsed case: $t=150$⁠. To minimize the difference between pulsed and continuous curves, the time was adjusted for the continuous cases, so as to match the neutral fixation rate: (A) $t=120$ and (B) $t=140$⁠. If using the actual times, the discrepancies between the two migration scenarios are even more pronounced. The curves are from the mathematical model with parameters $m=1$ and $h=1/2$⁠. Note the logistic scale for the y-axis.

One can extrapolate the fitness consequences of migration pulsedness over the entire genome, assuming a large number of independent loci, each with its particular selection regime. The distribution of $h$ and $s$ values over loci would vary greatly across organisms and contexts, but for illustration we can assume some plausible distribution of fitness effects (DFE) with $h=1/2$⁠, such as the one proposed in Martin and Lenormand (2015). This DFE, illustrated in Figure 5A (see Supplementary Information; section 4.1), features an excess of deleterious (maladapted) alleles, representative of situations where migration from the mainland mostly decreases fitness on the island. At any time, the mean fitness $w$ in the island population is the product of the fitness contributions of all loci, depending on their fixation status (see Supplementary Information; section 4.2).

Figure 5.

The pulsedness load. (A) An illustrative distribution of fitness effects (DFE) for mainland alleles across loci. Most alleles are slightly deleterious or neutral, some are beneficial (Martin & Lenormand, 2015). Three important classes of alleles are distinguished: $s<s_1$⁠, $s_1<s<0$ and $s>0$⁠. (B) The instantaneous rate of fixation (fixation flux) of mainland alleles as a function of the time since onset of migration, for continuous migration and two levels of migration pulsedness. The two inserts show the DFE of alleles that are currently going to fixation, at times $50$ and $500$⁠, for continuous migration and pulsed (⁠$n=20$ and $n=30$⁠) migration. (C) Genetic load (⁠$1-$ mean fitness) as a function of time since onset of migration, for continuous and pulsed migration. The excess of genetic load caused by migration pulsedness (the pulsed load) is highlighted as the shaded area. (D) Value of the pulsedness load (relative to the genetic load with continuous migration) as a function of time since onset of migration, for three different levels of migration pulsedness. Results are from the mathematical model with parameters $h=1/2$⁠, $m=1$⁠, and $N=100$⁠.

Following the onset of immigration, the island population mean fitness evolves through time, asymptotically settling at its equilibrium value when all mainland alleles have fixed, and thus controlled by the DFE. Under recurrent migration, regardless of the level of pulsedness, beneficial (adapted) alleles tend to go to fixation first, followed by neutral alleles, and ultimately deleterious (maladapted) alleles. These dynamics of allele fixation are shown in Figure 5B, with the DFE of currently fixing alleles, at two different times, shown as inserts. Introducing migration pulsedness slows down the fixation of beneficial and neutral variation, which decreases mean fitness in the island population, compared to what it should be under continuous migration. Conversely, it hastens the fixation of deleterious (maladapted) alleles with $s<s_1$⁠, which is also detrimental to mean fitness (see Figure 5A). As a result, migration pulsedness causes mean fitness to be lower than under continuous migration (Figure 5C), and we call this fitness deficit (i.e., additional genetic load) the “pulsedness load.” If most fitness variation is driven by adaptation to local habitat conditions, the pulsedness load translates into lower local adaptation. It must be stressed that if one simply extrapolated the fact that migration pulsedness decreases the effective migration rate, as is observed for neutral alleles, to the entire genome, one would reach the opposite prediction: migration pulsedness would increase mean fitness at any time.

As shown in Figure 5D, the pulsedness load can readily cause a $10\mathrm{\%}$ increase in the genetic load, compared to what is predicted under continuous migration. Interestingly, slightly deleterious alleles with selective values in the range $(s_1,0)$ have reduced fixation rates under pulsed migration, which has the opposite effect of boosting mean fitness: this contributes negatively to the pulsedness load. As a consequence, the pulsedness load might transiently decline (and possibly become negative), at the time when most alleles getting to fixation fall into this category (Figure 5D). Generically, the dynamics of the pulsedness load thus follows three consecutive phases over time, driven by the fixation of different categories of alleles (Figure 5A). The relative importance and timing of the three phases of course depends on the exact DFE.

Discussion

Temporal variability in migration is probably pervasive in nature. Plants, fungi, birds, mollusks and other marine invertebrates are most represented in the literature about variable dispersal, with many other taxa presumably concerned, including hominoids (Folinsbee & Brooks, 2007). A common form of migration temporal variability is migration pulsedness, i.e., the positively correlated migration of individuals, producing pulses of migration interspersed with periods of low or absent migration. However, classical evolutionary theory largely rests on the assumption of constant migration rates (e.g., Johnson et al., 2000; Mailund et al., 2012; Peniston et al., 2019; Rousset, 2013; Yeaman & Otto, 2011). Considering the growing evidence for non-constant migration processes, it is important to gain a general theoretical understanding of their evolutionary consequences.

We here proposed a novel approach and derived mathematical predictions regarding how the level of migration pulsedness should impact geneflow (the effective migration rate) at loci subject to an arbitrary selection regime. This is an important difference from earlier existing studies, which usually focused on neutral variation alone, or on some specific forms of selection. Though we here addressed pulsed migration patterns (Peniston et al., 2019; Yamaguchi & Iwasa, 2013), we can reasonably expect our predictions to apply more broadly to other forms of variable migration patterns that cause a temporal overdispersion in the number of migrants per generation. The latter is quantified by parameter $n$ in our model.

Our main finding is that the effect of migration pulsedness on geneflow depends not only quantitatively, but also qualitatively, on the selective regime considered. We find that migration pulsedness should decrease the effective migration rate (reduce geneflow) for alleles that are neutral or beneficial. However, for sufficiently deleterious (maladapted) alleles, the effect can be opposite. We found that there exists a threshold selection coefficient $s_1$⁠, of which we derived a mathematical expression, below which migration pulsedness on the contrary increases the effective migration rate (increases geneflow). In the additive case $h=1/2$⁠, the value of $s_1$ coincides with the well-known threshold for selection to be stronger than drift, corresponding to the divide between “truly deleterious” and “nearly neutral” alleles. The value of $s_1$ increases with the dominance level $h$⁠, so that an increase in geneflow with pulsedness is more likely for dominant alleles, and much less likely for recessive alleles.

The effect of migration pulsedness is therefore not homogeneous across different loci over the genome that experience different selection regimes. Migration pulsedness increases geneflow, and thus increases the speed and probability of fixation, for deleterious (maladapted) alleles, but it does the opposite for beneficial alleles. This homogeneizes fixation rates over the genome. Such a homogeneization of fixation rates across loci experiencing contrasting selection regimes could also be brought about by a reduced effective population size, i.e., a greater importance of drift relative to selection. However, we showed that a simple change in effective population size, retaining a continuous mode of migration, cannot adequately reproduce the fixation patterns created by pulsed migration. The action of migration pulsedness therefore leaves a signature that, at least in principle, could be distinguished from other forms of variability that reduce effective population size (random population fluctuations, random inter-individual fecundity variations, biased sex-ratios). The signature is more subtle. It is more comparable to that of deterministic processes such as directional changes of population size through time (Otto & Whitlock, 1997). In particular, an increased fixation of deleterious alleles and concomitant decreased fixation of beneficial alleles is analogous to what is predicted for shrinking populations (Otto & Whitlock, 1997). Our results also suggest that migration pulsedness would generate a relative excess of dominant deleterious alleles, and a relative deficit of recessive deleterious alleles (see Figure 2), which adds an extra dimension to its genomic signature.

Another consequence of migration pulsedness is that, by promoting the fixation of deleterious (maladapted) alleles, while slowing down the fixation of beneficial (adaptive) alleles, it overall reduces mean fitness and the level of local adaptation. We showed that migration pulsedness creates an additional genetic load, on top of the migration load expected under a continuous migration of similar intensity. This so-called “pulsedness load” can be non-negligible, and may easily represent a $10\mathrm{\%}$ increase in the genetic load, or even more, for plausible distributions of fitness effects and pulsedness values.

Our results generalize some existing results and are compatible with them. Most importantly, our predictions for neutral alleles are consistent with the few earlier population genetics studies, even though they used entirely different modelling approaches (Latter & Sved, 1981; Nagylaki, 1979; Rousset, 2013; Whitlock, 1992). We predict, as they did, that temporal variation in migration rates always decreases the effective migration rate for neutral alleles.

Even fewer studies have so far considered non-neutral variation. Peniston et al. (2019) investigated the consequences of pulsedness of migration in a mainland island model, like us, but they specifically considered source-sink situations: the island population was initially not viable on its own, and had to adapt in the face of maladaptive geneflow. This is quite a different scenario, as we here assumed that the island population was viable, not a sink. Our results thus complement the latter study, and we find, as they did, that migration pulsedness can importantly affect the dynamics of local adaptation. Even though direct comparisons are difficult, we remark that among all the scenarios they considered, for the one closest to ours, Peniston et al. (2019) report that more pulsedness migration hampers local adaptation, which is consistent with our predictions.

Yamaguchi and Iwasa (2013) and follow-up papers considered parapatric speciation through the accumulation of neutral mutations between two isolated populations, until prezygotic isolation develops. They used a formalism similar to ours, but in general they did not vary migration pulsedness and kept the number of migrants per migration event (⁠$n$⁠) at very low values. In one instance they did not, and found migration pulsedness to reduce geneflow, thus promoting parapatric speciation. This is in full agreement with our predictions for neutral variation. To go beyond the assumption of neutral mutations, we can apply our results to to the build-up of reproductive isolation under the effect of incompatible mutations. We can consider Dobzhansky-Muller–DM–mutations (Gavrilets, 2004) by seeing those as the limit, where $s\to 0^{\text{~}\text{-}}$⁠, and $h\to \mathrm{\infty }$ (strong overdominance). For instance, $s=-0.001$ and $h=100$ describes a mutation that causes a $-0.1$ fitness disadvantage only in heterozygotes. Extrapolating Figure 2 to those values, we would conclude that migration pulsedness should slightly decrease geneflow for them. This may seem paradoxical, as such mutations obviously benefit from arriving at higher frequency. However, their probability to go to fixation by chance is almost zero, unless their frequency reaches or exceeds $50\mathrm{\%}$⁠, a value that would require unrealistic large values of $n$⁠. We can tentatively posit that migration pulsedness should not much affect the accumulation of DM incompatibilities, even though establishing this would require a dedicated study.

Most of our predictions can be understood in terms of the trade-off between fewer migration events (rate $m/n$⁠), which reduces geneflow, and higher initial frequency after an event (⁠$n/(N+n)$⁠), which promotes geneflow. Increasing the pulsedness level $n$ decreases the number of migration events in inverse proportion, but the initial frequency at each event increases less than linearly (⁠$n/(N+n)$⁠). In general, and for neutral alleles in particular, the net effect of increasing $n$ is therefore detrimental to geneflow. Intuitively, this is because concentrating the arrival of alleles in time causes those mainland alleles to compete more among themselves, which does not help their fixation over local alleles. It is only when the fact of arriving in larger initial frequency is very beneficial to fixation that this effect can overcome that of a decreased number of migration events. This is the case for sufficiently deleterious alleles, whose chance of establishing is smallish from a low initial frequency, but much better if in higher frequency, owing to genetic drift. We have shown that sufficiently deleterious here means $s<s_1$⁠. The beneficial effect of a larger initial frequency is of course strengthened by positive frequency-dependent selection: this is why the value of $s_1$ increases with the level of dominance $h$ in a diploid context, as it generates positive frequency dependence at the allele level. An increase in geneflow with migration pulsedness would similarly be promoted by other sources of positive frequency dependence, such as aposematic color signaling (Endler, 1988): migration pulsedness could probably promote the successful migration of foreign aposematic morphs.

Our mathematical predictions held very well in stochastic simulations, for different population sizes, and even when the average migration rate (⁠$m$⁠) was quite large. However, in the latter case the possibility for migration pulsedness to increase geneflow is shifted toward more deleterious alleles, i.e., requires more stringent conditions than at lower migration rates. This can be understood intuitively in terms of mass-effects: with large migration rates, the influx of mainland alleles is strong enough to increase their frequency in the island on its own, regardless of the local selection-drift dynamics. In other words, deleterious alleles “pile-up” over successive migration events, and this can be an important driver of their eventual fixation. The latter effect is virtually unaffected by migration pulsedness; it solely depends on the average rate $m$⁠. As a consequence, the positive impact of migration pulsedness that we have identified, and explained just above, is in some sense diluted and gets relatively weaker. This is not the case for beneficial alleles, which tend to get fixed after a small number of migration events, so that the “piling-up” effect does not significantly contribute to their fixation. Further mathematical investigations allowing for large migration would be needed to validate this interpretation.

The finding that migration pulsedness all else equal favors the fixation of deleterious alleles and reduces local adaptation has direct implications for biological conservation, population management, and island biology. For instance, population reinforcement programs usually take the form of periodic releases of groups of individuals, whose frequency and intensity must both be optimized. Our results suggest that rare introductions of relatively large numbers of individuals, i.e., pulsed introduction patterns, would be more detrimental to population fitness and local adaptation, relative to more continuous fluxes of individuals, possibly compromising population viability and persistence. Extending our reasoning to island communities (Cowie & Holland, 2006), this also suggests that sporadic but potentially intense bouts of immigration, as brought about by rafting and human-driven invasions, could favor the establishment of relatively maladapted mainland species within communities competing for similar resources.

Of course, these impacts would combine with other effects, in particular demographical ones. It has been shown for instance that in small populations subject to Allee effects, pulsed migration patterns can be favorable to population establishment and persistence (Bajeux et al., 2019; Rajakaruna et al., 2013). Similarly, Peniston et al. (2019) found that migration pulsedness could sometimes favor local adaptation and persistence in populations that are demographic sinks (see also Gaggiotti & Smouse, 1996). These examples underline how complex the consequences of dispersal can be when demo-genetics effects are taken into account (Garant et al., 2007). It would be interesting to extend our model considering a broader range of demographic scenarios, and see how predictions are affected. For instance, the decrease in mean fitness through time as mainland alleles fix in the island population could have a negative impact on local population size. A decrease in $N$ could have various impacts on the probabilities of fixation. At the simplest level, it would increase the $n/N$ ratio, which would make the effects of pulsedness more and more pronounced. It would also increase the importance of drift, i.e., decrease $N_{e}s$⁠, which would shift the value of $s_1$⁠. How exactly this would impact the overall dynamics of the pulsedness load remains to be determined. At a more genetic level, we assumed for simplicity that alleles at different loci were unlinked and behaved independently. This will not always be the case, and it would be interesting to study the role of migration pulsedness for pairs of linked mutations that are in strong linkage disequilibrium among migrant individuals.

We conclude, along the lines of Peniston et al. (2019), by advocating for tests of predictions on migration pulsedness in empirical systems with well-designed experimental set-ups. We add that available data on migration flows, e.g., from monitoring of oceanic rafting, individual tracking of dispersal routes, or marine and aerial current models (Lagomarsino Oneto et al., 2020; Ser-Giacomi et al., 2015), combined with genomic data, may also provide opportunities for testing such predictions in the wild. Understanding the consequences of ecological variability is increasingly important in the present context of climate and habitat changes, and the spatio-temporal variability in migration patterns probably deserves to receive more attention.

Data availability

All simulation files and Java, R and Mathematica scripts are available on Zenodo (DOI 10.5281/zenodo.7133443).

Author contributions

V.C. and F.A. designed the research. B.L. wrote a first version of the codes and made preliminary analyses. F.A. and V.C. conducted the research and wrote the manuscript with inputs from L.M.

Conflict of interest: The authors declare no conflict of interest.

Acknowledgments

The authors are grateful to the OPAL infrastructure from Université Côte d’Azur for providing computing resources and support. Special thanks to cluster NEF (Inria SAM) and cluster BPI (ISA). The authors thank Thomas Guillemaud for regular discussions and comments on earlier versions of the manuscript. F.A. benefited from a PhD fellowship from Université Côte d’Azur (IDEX UCAJEDI, ANR-15-IDEX-01). V.C. acknowledges funding from ANR, project MUSCADO (ANR-20-CE02-0019).