Genetic load and extinction in peripheral populations: the roles of migration, drift and demographic stochasticity

We analyse how migration from a large mainland influences genetic load and population numbers on an island, in a scenario where fitness-affecting variants are unconditionally deleterious, and where numbers decline with increasing load. Our analysis shows that migration can have qualitatively different effects, depending on the total mutation target and fitness effects of deleterious variants. In particular, we find that populations exhibit a genetic Allee effect across a wide range of parameter combinations, when variants are partially recessive, cycling between low-load (large-population) and high-load (sink) states. Increased migration reduces load in the sink state (by increasing heterozygosity) but further inflates load in the large-population state (by hindering purging). We identify various critical parameter thresholds at which one or other stable state collapses, and discuss how these thresholds are influenced by the genetic versus demographic effects of migration. Our analysis is based on a ‘semi-deterministic’ analysis, which accounts for genetic drift but neglects demographic stochasticity. We also compare against simulations which account for both demographic stochasticity and drift. Our results clarify the importance of gene flow as a key determinant of extinction risk in peripheral populations, even in the absence of ecological gradients. This article is part of the theme issue ‘Species’ ranges in the face of changing environments (part I)’.


Introduction
Most outcrossing populations carry a substantial masked mutation load owing to recessive variants, which can contribute significantly to inbreeding depression in peripheral isolates or after a bottleneck. The extent to which the increased segregation (or fixation) of deleterious mutations due to drift (drift load) exacerbates extinction risk in isolated populations has been a subject of long-standing interest [1][2][3][4]. Theory predicts that moderately deleterious mutations contribute the most to genetic load and extinction in small populations [2,5]; however, the prevalence of such deleterious variants of mild or moderate effect and their dominance values remain poorly characterized, except for a few model organisms [6,7].
The relative risks posed by mutation accumulation and demographic stochasticity to a population depend crucially on its size, with theory suggesting that these may be comparable for populations in their thousands [1]. Additionally, environmental stochasticity-catastrophic events, as well as fluctuations in growth rates and carrying capacities, may dramatically lower extinction times [8]. Both demographic and environmental fluctuations, in turn, reduce the effective size of a population, making it more prone to fix deleterious alleles; the consequent reduction in fitness further depresses size, pushing populations into an 'extinction vortex', which is often characterized by a complex interaction between the effects of genetic drift, demographic stochasticity and environmental fluctuations [9].
Peripheral populations at the edges of species' ranges receive dispersal from the core to an extent which varies over space and time. Moreover, ranges may be fragmented owing to habitat loss and individual sub-populations connected to each other via low and possibly declining levels of migration. Under what conditions are such extinction vortices arrested by migration, and what are the genetic and demographic underpinnings of this effect, when it occurs?
Migration boosts numbers, mitigating extinction risk owing to demographic and environmental stochasticity, or at the very least, allows populations to regenerate after chance extinction. The demographic consequences of migration are especially important in fragmented populations with many small patches [10]: above a critical level of migration, the population may survive as a whole over long timescales even if individual patches frequently go extinct [11,12].
Migration also influences extinction risk by shifting the frequencies of fitness-affecting variants: the resultant changes in fitness may decrease or increase population size, thus further boosting or depressing the relative contribution of migration to allele frequency changes within a population, setting in motion a positive feedback which may culminate in extinction (when gene flow is largely maladaptive; e.g. [13,14]) or evolutionary rescue (if gene flow supplies variation necessary for local adaptation, or reduces inbreeding load; e.g. [5,15,16]).
The maladaptive consequences of migration have largely been explored for extended populations under spatially varying environments: here, gene flow typically hinders local adaptation, especially at range limits, leading to 'swamping' and extinction [17]. However, the consequences of gene flow for fitness, and consequently survival, are not always intuitive when the fitness effects of genetic variants are uniformly deleterious (or beneficial) across populations. For example, while gene flow may alleviate inbreeding load by preventing the fixation of deleterious alleles in small populations, it may also render selection against recessive mutations less effective by increasing heterozygosity. A striking consequence is that under a range of conditions, the fitness of metapopulations is maximized at intermediate levels of migration [18] and more generally, at intermediate levels of population structure [19].
A key consideration is whether or not gene flow is symmetric, i.e. whether some sub-populations are merely influenced by the inflow of genes from the rest of the habitat or if all sub-populations influence the genetic composition of the population as a whole [20]. Asymmetric dispersal is common at the geographical peripheries of species' ranges or on islands. Moreover, populations occupying small patches within a larger metapopulation with a wide distribution of patch sizes, or sub-populations with lower-thanaverage fitness (and consequently, atypically low numbers) may also experience predominantly asymmetric inflow of genes. Asymmetric gene flow allows for allele frequency differences across the range of a population even in the absence of environmental heterogeneity, e.g. when population sizes (and hence the efficacy of selection relative to drift) vary across the habitat. This, in turn, may generate heterosis or outbreeding depression across multiple loci, when individuals from different regions hybridize.
From a conceptual viewpoint, the consequences of asymmetric gene flow are typically simpler to analyse as we can focus on a single population, while taking the state of the rest of the larger habitat as 'fixed'. Such analyses are key to understanding more general scenarios where genotype frequencies and population sizes across different regions co-evolve.
Here, we analyse the eco-evolutionary dynamics of a single island subject to migration from a larger mainland in a scenario with uniform selection across the two populations, i.e. where fitness is affected by a large number of variants that are unconditionally deleterious. We ask: under what conditions can migration from the mainland alleviate inbreeding load, thus preventing 'mutational meltdown' and extinction of the island population? Further, how are the effects of migration mediated by the genetic architecture of load, i.e. by the genome-wide mutation target and fitness effects of deleterious variants? A key focus is to understand the coupled evolution of allele frequencies (across multiple loci) and population size: to this end, we consider an explicit model of population growth with logistic regulation, where growth is reduced by an amount equal to the genetic load.
While the effects of maladaptive gene flow on marginal populations have been studied under various models [21,22], there has been little work (under genetically realistic assumptions) on the ( possibly) beneficial effects of migration on inbreeding load and survival. In particular, modelling the polygenic nature of fitness variation is crucial, as changes in load (e.g. owing to migration) at any locus can affect all other loci by effecting changes in population size, which in turn influences the efficacy of selection across the genome.

Model and methods
Consider a peripheral island population subject to one-way migration from a large mainland. Individuals are diploid and carry L biallelic loci that undergo bidirectional mutation between the wild-type and deleterious state at rate u per generation per individual per locus in either direction. Mutations have the same fitness effects on the mainland and island, i.e. there is no environment-dependent fitness component.
Deleterious variants across loci affect fitness multiplicatively (no epistasis): individual fitness is e À P L j¼1 s j ðX j þh j Y j Þ , where X j (Y j ) equals 1 if the jth locus is homozygous (heterozygous) for the deleterious allele, and is zero otherwise. Here, s j is the homozygous selective effect and h j the dominance coefficient for the deleterious allele at locus j. We assume 0 ≤ h j ≤ 1/2, so that deleterious alleles are (partially) recessive. In individual-based simulation (see the electronic supplementary material, appendix A), we use the form (1 − hs) n (1 − s) m 1 m 0 (which is equivalent to the above fitness function for small s) where n represents the number of heterozygous loci, m represents the number of loci homozygous for the deleterious allele and m 0 represents the number of homozygous loci for the wild-type allele with n + m + m 0 = L.
In each generation, a Poisson-distributed number of individuals (with mean m 0 ) migrate from the mainland to the island. For simplicity, we assume that the mainland population is large enough that deleterious allele frequencies among migrants (denoted by fp ðmÞ j g) are close to the deterministic predictions for a single locus under mutation-selection equilibrium.
We assume density-independent selection on individuals on the island and density-dependent population growth with logistic regulation, where the baseline growth rate is reduced by the genetic load: the population size n t in any generation t is then Poisson-distributed with mean equal to n tÀ1 exp½r 0 ð1 À n tÀ1 =KÞW, where n t−1 is the population size in the previous generation, r 0 the baseline growth rate, K the carrying capacity, and W the mean genetic fitness. The n t individuals in the tth generation are formed by randomly sampling 2n t parents (with replacement) from the n t−1 individuals in the previous generation with probabilities equal to their relative fitnesses, followed by free recombination between parental haplotypes to create gametes. Diploid offspring are then formed by randomly pairing gametes.

(a) Simulations
We carry out two kinds of simulations: individual-based simulations that explicitly track multi-locus genotypes of all individuals on the island, and simulations that assume linkage equilibrium (LE) and neglect inbreeding. The latter kind of simulations are computationally less intensive as they only track allele frequencies at the L loci and the size of the population. However, they make two simplifying assumptions: first, that genotypes at any locus are in Hardy-Weinberg proportions (i.e. no inbreeding); second, that any statistical associations, e.g. between the allelic states of different loci (linkage disequilibria; LD) or between the probability of identity by descent, and consequently homozygosity, at different loci (identity disequilibria; ID)are negligible. Then, individual genotypes are simply random assortments of deleterious and wild-type alleles, and can be generated, e.g. in a simulation, by independently assigning alternative allelic states to different loci with probabilities equal to the allele frequencies. Details of the simulations are provided in the electronic supplementary material, appendix A.
Because selection pressures are identical on the mainland and island, systematic differences in allele frequencies or homozygosity between the two populations across multiple loci (which would generate LD and ID, respectively) must arise solely owing to differences in population size, which would cause the efficacy of selection to be different on the mainland and island. In general, we expect LD and ID to be negligible when all ecological and evolutionary processes are slower than recombination [14]: this may not hold, however, when populations are small (and drift significant), making it necessary to evaluate how associations between deleterious variants affect extinction thresholds. In the rest of the paper, we will only show results of the allele frequency simulations (assuming LE and zero inbreeding); we compare these with individualbased simulations and discuss how LD and ID affect population outcomes in the electronic supplementary material, appendix B.

(b) Joint evolution of allele frequencies and population size
Assuming LE and no inbreeding, the evolution of the island population is fully specified by how allele frequencies {p j } and the population size n co-evolve in time t. If all evolutionary and ecological processes (except recombination) are slow, we can describe this co-evolution in continuous time. We rescale all rates by the baseline growth rate r 0 and population sizes by the carrying capacity K. This yields the following dimensionless parameters: τ = r 0 t, S = s/r 0 , M 0 = m 0 /(r 0 K), U = u/r 0 , N t = n t /K and an additional parameter ζ = r 0 K which governs the strength of demographic fluctuations. The joint evolution of allele frequencies {p j } and the scaled population size N in continuous time is described by the following equations (see also [14]): S j p j ð2h j q j þ p j Þ: ð2:1Þ The random processes λ N and l p j satisfy E½l p j ¼ E½l The four terms in the first equation correspond (in order of appearance) to changes in allele frequency due to selection, mutation, migration and drift. Note that the strength of migration is inversely proportional to population size N, reflecting the stronger (relative) effect of migration on the genetic composition of smaller, as opposed to larger, island populations. The second equation describes the evolution of population size N: the first term describes changes in N under logistic growth, where the growth rate is reduced by a factor proportional to the log mean fitness (i.e. the genetic load); the second term captures the effect of migration; the third term corresponds to demographic fluctuations (whose variance is proportional to N, the size of the population). These equations capture a key feature of polygenic eco-evolutionary dynamics-namely, that the evolution of allele frequencies at different loci is coupled via their dependence on a common N, which in turn is influenced by the degree of maladaptation at all loci via R g . Thus, allele frequencies do not evolve independently, even if allelic states at different loci are statistically independent at any instant (under LE).
For fixed N and under LE and IE, the joint distribution of allele frequencies at mutation-selection-migration-drift equilibrium is a product of the single-locus distributions. This was given by Wright [23], and for a given locus j is (in terms of scaled parameters): e À2zNSp j ½p j þ2hð1Àp j Þ : ð2:2Þ Integrating over this distribution yields the expected allele frequency Eðp j jNÞ and the expected heterozygosity Eð2p j q j jNÞ at any locus, and thence the expected total load EðR g jNÞ ¼ P j S j ½Eðp j jNÞ À ð1=2 À hÞEð2p j q j jNÞÞ (scaled by the baseline growth rate r 0 ) for fixed N.
However, in reality, N is not fixed and will fluctuate-both owing to random fluctuations in fitness (due to the underlying stochastic fluctuations of {p j } at equilibrium) as well as in the reproductive output of individuals (demographic stochasticity). Thus, N itself follows a distribution. While we can write down an equation for the stochastic co-evolution of N and {p j }, no explicit solution for the joint equilibrium distribution is possible unless mutation rates are strictly zero (as assumed by Szép et al. [14]). Thus, we must employ various approximations to describe the coupled dynamics of population sizes and allele frequencies.

(c) Approximate semi-deterministic analysis
We expect the population size distribution to be sharply peaked around one or more values {N * } if demographic fluctuations are weak (i.e. for ζ = r 0 K ≫ 1) and if fluctuations in royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210010 mean fitness are also weak (i.e. for large L). Such sharply peaked distributions correspond to populations that transition only rarely between alternative peaks, i.e. where the alternative peaks {N * } of the size distribution represent 'metastable' states, allowing allele frequencies sufficient time to equilibrate at any N * . Then, the genetic load would be close to the expected value EðR g jN Ã Þ under mutation-selectiondrift-migration balance (given N * ) when populations are in one of the alternative metastable states, though not necessarily while they transition between states. In order to determine {N * }, we postulate that these must represent stable equilibria of the population size dynamics, neglecting demographic stochasticity and assuming that mean fitness (given N * ) is close to the expectation under mutation-selection-drift-migration balance for that N * . Then we have: The equality follows from equation (2.1), by setting the third (noise) term to zero and assuming that R g is close to its equilibrium expectation EðR g jN Ã Þ, given N * . The second inequality is the condition for N * to be a stable equilibrium: this means that populations starting at an arbitrary N in the vicinity of N * would evolve towards this equilibrium size. Equation (2.3) can be solved numerically to obtain the equilibria {N * }, which, under the above assumptions, will be close to the peaks of the population size distribution. As we see below, depending on the parameter regime, there may be one or two stable equilibria of equation (2.3), corresponding to population size distributions that are unimodal or bimodal. Bifurcations (i.e. critical parameter thresholds where one of the two stable equilibria vanishes) thus correspond to qualitative transitions in the state of the population.
We will refer to the analysis above as a 'semi-deterministic analysis' as it accounts for the stochastic effects of genetic drift on allele frequencies and load (via the expectation E½R g jN Ã ) but neglects demographic stochasticity. In general, we expect the semi-deterministic analysis to become more accurate for larger ζ = r 0 K (see also the electronic supplementary material, appendix B), since increasing ζ, which implies higher number of births per generation, results in weaker demographic fluctuations [14].

(a) Metastable populations and extinction times in the absence of migration
We first consider peripheral populations in the absence of migration. Such populations necessarily become extinct in the long run: however, depending on the mutation target for deleterious mutations (relative to the baseline growth rate r 0 ), the fitness effects of mutations and the carrying capacity of the island, extinction times may be very long and populations metastable. We can use equation (2.3) to gain intuition for the conditions for metastability under zero migration. Setting M 0 = 0, it follows that there is always an equilibrium at N = 0 (corresponding to extinction): this equilibrium is stable for L(S/2) > 1 where S = s/r 0 , and unstable otherwise (see the electronic supplementary material, appendix B for details). There may exist a second equilibrium at if the equilibrium genetic load is lower than the baseline growth rate.
Because the equilibrium load E½R g jN Ã , M 0 ¼ 0 depends on four independent parameters (which we choose as ζS = Ks, ζU = Ku, h and 2LU = 2L(u/r 0 )), mapping the conditions for metastability boils down to asking: in the absence of demographic fluctuations, where in this four-dimensional parameter space, can we find populations with a non-zero equilibrium size or sufficiently low load (figure 1a)? In reality, there is a fifth parameter ζ = r 0 K, which governs demographic royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210010 stochasticity: thus, all other parameters being equal, extinction times will be longer when demographic stochasticity is weaker, i.e. r 0 K larger (figure 1b). For simplicity, we assume that deleterious alleles anywhere on the genome have the same selective effects and dominance coefficients; this assumption is relaxed in the electronic supplementary material, appendix D. We use Ks c to denote the critical selection strength per homozygous deleterious allele (scaled by the carrying capacity K), such that a non-zero equilibrium N * exists for Ks > Ks c but not for Ks < Ks c . Figure 1a shows Ks c as a function of 2LU = 2L(u/r 0 ), the total mutation rate relative to the baseline growth rate, for different Ku (different colours) for nearly recessive (h = 0.02; solid lines) and additive (h = 0.5; dashed lines) alleles.
For populations to be metastable, the total genetic load must be less than the baseline growth rate. The total load scales with the mutation target L; the load per locus is approximately 2u for strongly deleterious variants (u/hs ≪ 1 and Ks ≫ 1). For smaller Ks, drift will typically inflate load above this deterministic expectation when deleterious alleles are additive (h ∼ 1/2) but may also reduce load (via more efficient purging) when alleles are recessive (h ∼ 0). This is reflected in figure 1a: the threshold Ks c required for metastable populations is higher for additive deleterious alleles (dashed lines) than recessive alleles (solid lines). Moreover, because genetic load in this drift-dominated regime depends only weakly on the mutational input Ku, the threshold Ks c is largely independent of Ku (different colours). Finally, for all parameter combinations, the threshold Ks c increases as the mutation target becomes larger: this simply reflects the fact that for the total load to be less than the baseline growth rate, the load per locus must be lower (requiring stronger selection) if deleterious variants segregate at a greater number of loci. Accordingly, for very large mutation targets 2LU ¼ 2Lðu=r 0 Þ * 1, the total load will exceed r 0 (and populations will fail), irrespective of the strength of selection against deleterious mutations.
The critical selection thresholds for metastability shown in figure 1a are computed by neglecting demographic stochasticity, i.e. by assuming ζ = r 0 K to be very large. However, for moderate ζ, stochastic fluctuations in reproductive output from generation to generation may accelerate extinction: this effect can be especially significant in smaller populations as these tend to fix more deleterious alleles, which further reduces fitness and size, thus rendering populations even more vulnerable to stochastic extinction.
To investigate how demographic stochasticity contributes to extinction, we compare populations residing on islands with different carrying capacities K, but characterized by the same values of Ks, Ku, 2LU = 2L(u/r 0 ) and h. Mathematically, this involves taking the limit s → 0, u → 0, K → ∞, L → ∞, while holding Ks, Ku, 2LU constant: then, increasing K has the sole effect of weakening demographic stochasticity (by increasing r 0 K). Populations are initially perfectly fit, but accumulate deleterious variants over time, eventually becoming extinct owing to the combined effects of genetic load and demographic stochasticity. Figure 1b shows the extinction 'half-time' T 1/2 = r 0 t 1/2 (scaled by the baseline growth rate r 0 )-the time by which precisely half of all 1000 simulation replicates are extinct, as a function of 2LU for various K for nearly recessive (main plot) and additive alleles (inset). All results are from allele frequency simulations (assuming LE and zero inbreeding). Dashed vertical lines indicate the threshold 2LU above which metastable populations (with N * = 1 − E[R g |N * ] > 0) cannot exist (even for large r 0 K).
We find that extinction times increase with increasing K for all parameters. However, for parameter combinations that correspond to extinction in the large r 0 K limit (i.e. to the right of the dashed lines), this increase is approximately linear in K, while for parameters leading to metastability in the large r 0 K limit (left of dashed lines), this increase is faster than linear. In fact, in the metastable regime, even with a carrying capacity K of a few thousand individuals, demographic stochasticity will be sufficiently weak and extinction times large enough that isolated populations can persist over geological timescales: for these parameter regimes, it is environmental fluctuations (which our model ignores), rather than mutation load or demographic stochasticity, that will primarily influence extinction risk.
We also compute the average (scaled) population size N = n/K in the metastable state (figure 1c). This declines with increasing 2LU, i.e. as we approach the threshold for loss of metastability, and is close to the semi-deterministic prediction (dashed black lines) for all values of K.

(b) Effect of migration on equilibrium population sizes and genetic load
We now consider how migration from a large mainland influences population dynamics on the island, for different genetic architectures of load, i.e. given certain selective effects and dominance coefficients of deleterious alleles. As before, we assume that all deleterious alleles have equal selective effects and dominance values, relegating the discussion of more general scenarios, where alleles with different fitness effects segregate, to the electronic supplementary material, appendix D. We first identify critical parameter thresholds associated with qualitative changes in population outcomes for the two extremes: nearly recessive (h = 0.02 in figure 2) and additive (h = 0.5) alleles. We then consider how these thresholds depend on the dominance of deleterious alleles (figure 3). Figure 2 illustrates the effect of migration on population sizes and load on the island for three parameter combinations, corresponding to weakly deleterious (Ks < Ks c ; left column), moderately deleterious (Ks * Ks c ; middle) and strongly deleterious (Ks ≫ Ks c ; right) nearly recessive alleles (h = 0.02). Recall that Ks c is the selection threshold (obtained by neglecting demographic stochasticity) such that populations rapidly go extinct in the absence of migration for Ks < Ks c , but can be metastable (even without migration) for Ks > Ks c and r 0 K → ∞ ( figure 1). For Ks < Ks c , there exists a single stable equilibrium at all migration levels (figure 2a), with the corresponding population size approaching zero as migration becomes rarer. Accordingly, the stochastic distribution of N exhibits a single peak for all m 0 ( figure 2d). An increase in migration causes the genetic load to decrease (inset, figure 2a) and the equilibrium population size to increase. The stochastic distribution of population sizes is approximately exponentially distributed for low m 0 (corresponding to a sink state), but shifts towards higher N as m 0 increases, and is approximately normally distributed about the deterministic equilibrium N * (indicated by empty triangles in figure 2d) for large m 0 .
Increasing migration has two effects in this case: it reduces genetic load, typically by reducing homozygosity-a genetic effect, but also increases population numbers-a demographic effect. The genetic effect of migration, i.e. the dependence of the expected load EðR g jNÞ ¼ P j S j ½Eðp j jNÞ À ð1=2 À hÞ Eð2p j q j jNÞ on m 0 , can be further decomposed into effects on the expected frequencies Eðp j jNÞ and expected heterozygosity Eð2p j q j jNÞ of deleterious alleles. For weakly deleterious recessive alleles, the expected frequency decreases with increasing migration for low m 0 , is minimum at intermediate m 0 , and then increases as m 0 increases further (not shown). However, genetic load still decreases monotonically with increasing migration (inset, figure 2a), owing to the much sharper increase in heterozygosity with migration. 2c) corresponds to a sink state with the associated population size approaching zero as m 0 → 0. The other equilibrium (depicted in red) corresponds to a large population, which can persist at a finite fraction of carrying capacity even as migration becomes exceedingly rare, i.e. has size N * , which approaches 1 À E½R g jN Ã , M 0 ¼ 0 . 0 as m 0 → 0. The existence of two stable equilibria at low m 0 translates into bimodal population size distributions (in black and brown in figure 2f), with one peak close to extinction and the other at the semi-deterministic N * (depicted by filled triangles). Population size N(t) also exhibits characteristic dynamics (see time series in black and brown in figure 2i): populations tend to remain 'stuck' in either the sink state or the large-population state for long periods of time, typically exhibiting rather small fluctuations about the characteristic size associated with either state, and only transitioning rarely between states.
The sink and large-population states can also be thought of as 'high-load' (i.e. genetic load r g greater than the baseline growth rate r 0 ) and 'low-load' (r g < r 0 ) states. Thus, in this parameter regime, populations exhibit a genetic Allee effect, wherein load is sufficiently low and net growth rates positive only above a threshold population size: therefore, starting from an initially empty island, populations cannot grow deterministically (and persist instead as migrationfed sinks), even though large populations can maintain themselves, at least in the absence of demographic fluctuations. Transitions between the sink (high-load) and large-population (low-load) states are thus inherently stochastic, arising owing to demographic fluctuations (which are aided by higher migration) rather than owing to a systematic drive towards the alternative state.
Changes in m 0 have qualitatively different effects on the two states: at low m 0 , increasing migration reduces load and thus increases numbers in the sink state (blue curve), but has the opposite effect on the large-population equilibrium (red curves). This is owing to the non-monotonic dependence of the frequency of deleterious recessive alleles on population size: increasing size causes drift to become weaker relative to selection but also reduces homozygosity, so that fewer deleterious alleles are exposed to selection; frequencies are thus minimum at intermediate population sizes, reflecting the tension between these two opposing effects. In small sink populations, migration from the continent reduces homozygosity as well as deleterious allele frequencies, thus reducing load. However, in larger populations, increasing migration reduces the homozygosity but raises the frequency of deleterious alleles (since deleterious alleles are purged less efficiently in very large mainland populations than in intermediate-sized or large island populations); the overall effect of migration is thus to increase load in the large-population state. Above a critical migration threshold, which we denote by m c,1 , the sink equilibrium vanishes. Thus, for m > m c,1 , populations always have a positive growth rate and reach a finite fraction of carrying capacity, regardless of starting size. As before, this qualitative change in the (semi-deterministic) population size dynamics at a critical migration level has its analogue in a qualitative change in the stochastic distribution of population sizes (figure 2f): increasing migration causes the population size distribution to change from bimodal to unimodal (with a sole peak at the large-population equilibrium N Ã ¼ 1 À E½R g jN Ã ). At higher migration rates, i.e. as we approach the critical migration threshold m c,1 starting from low m 0 , turnover between the sink and large-population states becomes more rapid-note the shorter intervals between transitions in the brown versus black time series in figure 2i.

(e) Moderately deleterious recessive alleles
Consider now the case where selection is stronger than the threshold Ks c but not too strong, i.e. Ks * Ks c (middle column in figure 2). As in the Ks ≫ Ks c case, there exist two alternative stable equilibria at low migration levels, one corresponding to the sink state (blue lines in figure 2b) and the other to the large-population state (red lines). As before, increasing migration alleviates inbreeding load and increases population size in the sink state but elevates load (by hindering purging) and decreases population size in the largepopulation state. Unlike in the Ks ≫ Ks c state, the latter effect is much stronger: thus, above a critical migration threshold, which we denote by m c,2 , it is the large-population equilibrium that vanishes, so that there is only a single (sink) equilibrium for m > m c,2 .
The analogous change in the population size distribution with increasing migration (figure 2e) is somewhat subtle: for the lowest value of m 0 (black), the distribution is clearly bimodal, with most of the weight close to N = 0 (sink state) and a very small peak at the large-population equilibrium. As migration increases, the distribution of sizes associated with the sink state widens, while the peak corresponding to the large-population equilibrium shifts towards lower N (e.g. brown plot). At high migration levels (orange and blue plots), there is no distinct second peak and only the sink state persists. In this state, the genetic load exceeds r 0 on average; however, its distribution and the corresponding distribution of N is quite wide.
If migration increases further, then the average load associated with the sink state continues to fall; at a third threshold (denoted by m c,3 ; here approx. 1.33), the load again becomes lower than the baseline growth rate r 0 or, alternatively, scaled load R g = r g /r 0 less than 1 (depicted by a dashed line in figure 2b). Thus, for m 0 > m c,3 , populations can grow, starting from small numbers, and reach a finite fraction of carrying capacity. However, unlike the large-population state that emerges at higher Ks, in this case, populations are highly dependent on migration and would rapidly collapse (owing to the fixation of deleterious alleles) if cut off from the mainland.
Population size dynamics (figure 2h) are characterized, as in the Ks ≫ Ks c case, by occasional transitions between the sink state and the large-population state for low values of m 0 (black), with transitions becoming more frequent with increasing m 0 (orange). Transition times are much longer than in the Ks ≫ Ks c case (note the values on the x-axis); thus transitions are unlikely to occur over realistic timescales, and populations will typically be observed in the sink state. At higher migration levels, there are no obvious transitions, with population sizes and load fluctuating (with some skew) about a mean value.

(f ) Additive alleles
With additive effects (h ∼ 0.5), any deleterious allele experiences the same selective disadvantage, irrespective of whether it appears in the heterozygous or homozygous state: thus, there is no purging in smaller populations (which have higher homozygosity) and allele frequencies decrease monotonically with increasing population size. As a result, migration from the larger mainland always decreases load (by decreasing the deleterious allele frequency) and increases the size of the island population.
This implies that there are only two qualitatively distinct regimes: with weakly deleterious alleles (Ks less than the corresponding Ks c ; figure 1a), populations tend to extinction as migration declines. There is a single equilibrium of the semideterministic population size dynamics for all m 0 or analogously, a single peak of the stochastic population size distribution, which shifts towards higher sizes as m 0 increases.
For strongly deleterious alleles (Ks > Ks c ), population size N and load R g are largely insensitive to migration, since populations can always grow, starting from small numbers, at least when demographic stochasticity is unimportant (i.e. for r 0 K ≫ 1). With smaller r 0 K, population sizes and load depend weakly on migration, since demographic stochasticity may depress N and inflate the effects of drift (even with Ks > Ks c ). Moreover, in this regime, populations are also prone to stochastic extinction for m 0 < 1/2 [14], such that the distribution of N is inherently bimodal. Thus, where demographic stochasticity is significant, a low level of migration may be needed for stable populations even with strong selection against additive alleles.
The analysis outlined here (based on equation (2.3)) also applies to loci with a distribution of effects: in the electronic supplementary material, appendix D, we consider examples where a fraction of deleterious mutations are additive and the remaining fraction recessive. We also compare the results of allele frequency simulations with those of individual-based simulations with unlinked loci (electronic supplementary material, appendix C). These show fairly close agreement, suggesting that LD and ID do not significantly affect allele frequency dynamics, at least for the typical parameter values considered here: this is also consistent with earlier work, which suggests only a modest effect of disequilibria on background selection in sub-divided populations under soft selection [19].
royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210010 (g) Disentangling genetic versus demographic effects of migration on recessive alleles In summary, given a mutation target 2LU & 1 and assuming equal-effect loci, there is a critical selection threshold Ks c , such that large populations are metastable in the m 0 → 0 limit for Ks > Ks c but not for Ks < Ks c . For Ks < Ks c , there is a genetic Allee effect at low migration levels and with partially recessive alleles: populations cannot grow subsequent to recolonization of an initially empty island, and persist only as demographic sinks until a chance fluctuation increases numbers sufficiently that load can be purged and the alternative (large-population) equilibrium attained; such large populations can then be maintained over long periods of time (figure 2i). For very low dominance values (h & 0:15), there is a second threshold Ks c,2 with Ks c,2 > Ks c (see below), which separates parameter regimes characterized by qualitatively different effects of migration on population outcomes: for Ks > Ks c,2 , increasing migration destabilizes the sink state, so that only the large-population state persists above a migration threshold m c,1 , whereas for Ks c < Ks < Ks c,2 , increasing migration destabilizes the large-population state, so that only the sink state persists above a threshold m c,2 . However, such migration-fed sink populations can be quite large: above a third threshold m c,3 , populations may even have sufficient heterozygosity to again attain a low-load (r g < r 0 ) state.
To what extent can we attribute such qualitative changes in population outcomes to the genetic versus demographic effects of migration? As before, one approach is to compare critical thresholds for populations with the same scaled parameters Ks, Ku, 2LU = 2L(u/r 0 ) and h, while increasing the carrying capacity K (simultaneously increasing L and lowering u, s). Then, as K increases, the demographic effects of migration (which depend on the dimensionless parameter M 0 = m 0 /(r 0 K); see equation (2.1)) can be neglected, while its genetic effects on load (which depend on Ks, Ku and m 0 ) remain important. Figure 3 shows these comparisons for two dominance values-h = 0.02 (figure 3a,b) and h = 0.1 (figure 3c,d ). Note that we only consider predictions of the semi-deterministic analysis (equation (2.3)), which assumes that population sizes are sharply clustered around the peaks of the distribution and transitions between peaks are infrequent. This assumption clearly breaks down close to transition thresholds (e.g. see figure 2e,f ); thus, thresholds observed in simulations may differ somewhat from semi-deterministic predictions (details in the electronic supplementary material, appendix B). Moreover, the semi-deterministic analysis neglects all demographic stochasticity, and thus does not account for the fact that changes in K will affect not just the (systematic) demographic effect of migration but also the (stochastic) effect of demographic fluctuations. Nevertheless, the semi-deterministic analysis is useful as it allows us to explore qualitative dependencies of critical thresholds on the underlying parameters without resorting to time-consuming simulations. Figure 3a,c show the critical selection thresholds Ks c (as in figure 1a) and Ks c,2 versus the scaled mutation target 2LU = 2L(u/r 0 ) for the two values of h, for various carrying capacities (with Ks, Ku, 2LU = 2L(u/r 0 ) held constant, as K is varied). The threshold Ks c , which relates to population outcomes in the zero migration limit, is (by definition) independent of K in the semi-deterministic setting, as increasing K merely weakens the demographic effects of migration.
The threshold Ks c,2 decreases as K decreases, i.e. as the demographic effects of migration become stronger. This means that when alleles are moderately deleterious, i.e. characterized by a certain intermediate value of Ks, populations are stabilized (i.e. rescued from recurrent collapse into the high-load sink state) by increasing migration more easily on smaller islands, where the demographic effects of migration are stronger. Figure 3b,d shows the critical migration thresholds m c,1 , m c,2 and m c,3 versus 2LU for a given Ks value (chosen to be 50). Here, with h = 0.02 for example (figure 3b), increasing migration causes the sink state to vanish for 2LU & 0:5, but degrades the large-population (low-load) state for 2LU * 0:5. The threshold m c,1 (solid lines) is highly sensitive to K: populations can be stabilized at a finite fraction of carrying capacity at much lower levels of migration on smaller islands (given Ks, Ku), suggesting a key role of the demographic effects of migration. We can obtain an explicit expression for the threshold m c,1 in the limit K → ∞ (so that M 0 = m 0 /(r 0 K), which governs the demographic effects of migration, is negligible): m c,1 ¼ LSp ðmÞ À 1 4½1 À LSp ðmÞ ðp ðmÞ þ 2hð1 À p ðmÞ ÞÞ as K ! 1: ð3:1Þ Interestingly, this threshold depends only on the load, LSp (m) ( p (m) + 2h(1 − p (m) )), and breeding value, LSp (m) , among migrants, and is independent of Ks and Ku: this simply reflects the fact that the growth rate of a very small population (just after recolonization) depends primarily on the genetic composition of founders (and not selection on their subsequent descendants). Thus, the critical level of migration required to prevent a genetic Allee effect in the limit of very large carrying capacities is also independent of Ks and Ku. The threshold m c,2 , which signals the collapse of the largepopulation state, is less sensitive to K: this is consistent with our expectation that demographic effects of migration should be less important when numbers are larger. There is, nevertheless, a moderate decrease in m c,2 with K, which can be rationalized as follows: the (detrimental) genetic effects of migration on the large-population state are more effectively compensated by its demographic effects when islands are smaller (lower carrying capacity), allowing the largepopulation state to persist despite higher levels of migration on such islands. Finally, the third threshold m c,3 , which signals the emergence of a migration-dependent low-load state (at large 2LU) is also highly sensitive to the demographic effects of migration, and becomes unrealistically high in the large K limit, where increasing migration can only shift allele frequencies but has little effect on population numbers (relative to carrying capacity).
A comparison of the top and bottom rows in figure 3, which correspond respectively to h = 0.02 and h = 0.1, shows that as deleterious alleles become less recessive (larger h), the parameter regime in which increasing migration eliminates the large-population equilibrium shrinks drastically, and only emerges for very high genome-wide mutation rates-with h = 0.1, the second threshold Ks c,2 exists only for 2LU * 0:8 in figure 3c. This is consistent with the fact that purging of recessive mutations in smaller populations (which allows them to evolve significantly lower load than mainland populations) is only effective for very low h; thus, gene flow from the mainland is less detrimental to the large-population royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210010 equilibrium in island populations for less recessive alleles. In fact, for h * 0:15, increasing migration always causes the sink equilibrium to vanish, irrespective of Ks. Moreover, the migration threshold m c,1 at which the sink state vanishes falls with increasing h (compare figure 3b,d). Thus, we observe a genetic Allee effect only at very low migration rates for moderately recessive alleles. This is consistent with the fact that inbreeding load in small populations (owing to excess homozygosity) becomes lower as alleles becomes less recessive, and is alleviated by even low levels of migration.

Discussion
A key parameter governing the fate of peripheral populations is 2LU = 2L(u/r 0 ), the genome-wide deleterious mutation rate relative to the baseline rate of population growth: low-load (large-population) states are possible only for 2LU & 1, provided selection against deleterious variants is sufficiently strong and/or migration high. Conversely, for 2LU * 1, populations exist only as demographic sinks, irrespective of selection strength. The parameter 2LU is a measure of the 'hardness' of selection and can be small either if the total mutation rate 2Lu (which determines total load in the absence of drift) is small or, more realistically, if the growth rate r 0 , i.e. the logarithm of the baseline fecundity, is high (corresponding to the soft selection limit).
Our analysis highlights qualitatively different effects of migration on population outcomes, depending on fitness effects and the total mutation target of deleterious variants (figures 2 and 3). For example, with 2LU = 0.5 (which corresponds to a % 50% reduction in growth rate owing to genetic load in a deterministic population), typical Ks must be at least approximately 5 for recessive and approximately 10 for additive alleles, if populations are to be metastable in the absence of migration. For weaker selection, gene flow from the mainland is beneficial, and aids population survival by hindering the fixation of deleterious alleles. However, for stronger selection and with recessive alleles, the fitness and size of 'low-load' island populations actually declines with increasing migration from the mainland (owing to higher deleterious allele frequencies in the latter). In the most extreme scenario, where load is primarily owing to segregation of recessive alleles of moderate effect, intermediate levels of migration may actually increase load so much that populations degenerate into high-load demographic sinks (figure 2b; electronic supplementary material, figure S4).
We identify two regimes in which peripheral populations can maintain stable numbers at a substantial fraction of carrying capacity, with qualitatively different roles of migration in the two. When selection is strong, i.e. Ks ≫ Ks c (for additive alleles) or Ks > Ks c,2 (for recessive alleles) for a given 2LU, genetic load is low and populations stable, largely independently of migration. In this case, low levels of migration (typically & 1 migrant per generation; note typical values of m c,1 in figure 3) are sufficient to prevent a genetic Allee effect, should demographic stochasticity or chance fluctuations in load drive population numbers down. On the other hand, with weakly or moderately deleterious alleles, stable populations rely on rather high levels of migration ( ≫1 migrant per generation; note typical m c,3 in figure 3) and are only weakly differentiated with respect to the mainland. Here, migration is essential for maintaining heterozygosity and preventing fixation of deleterious alleles, even though numbers are relatively large.
Both classical quantitative genetics and analyses of allele frequency spectra suggest that most mutation load is owing to weakly deleterious alleles [24,25]. The weak selection on deleterious mutations may be strong relative to random drift in the species as a whole [24], but is likely to be dominated by random drift within local demes (i.e. Ks < 1 in our notation). If this is so, then extinction can be avoided only if many migrants enter demes in each generation. Fortunately, this is generally the case: F ST is typically small, implying that m 0 is between 0.5 and 5 ([26]; note that m 0 is the number of migrant genes, corresponding to 2 Nm in the usual diploid notation). Thus, while selection can act effectively to suppress the mutation load in a well-connected metapopulation, demes that receive few migrants will be vulnerable to the accumulation of load owing to weakly selected mutations. Moreover, when fitness has additional environment-dependent components, local adaptation must depend on alleles of intermediate effect and is hindered by high migration, especially for marginal habitats within the metapopulation [14].
What implications might our work have for the conservation of natural populations? Provided that several migrants are exchanged per generation, selection against deleterious alleles can be effective across the whole population. Indeed, subdivision into small subpopulations can help purge deleterious recessives, making selection more effective than with panmixia. Thus, random drift would lead to a severe load only if local demes are highly isolated-in which case environmental fluctuations are more likely to cause extinction than the gradual accumulation of weakly deleterious mutations [9]. Our work implies that an intermediate rate of migration minimizes mutation load, by preventing extinction of local populations, and yet still allowing some purging. However, the extinction risk arising from environmental fluctuations (which we underestimate by including only demographic stochasticity) favours higher migration. Conversely, local adaptations require selection that is stronger than drift within local demes (Ks > 1) [14]; if this is a concern, then substantial deme sizes are required in the long term.
Our model makes various assumptions: first, we take mutation rates to and from the deleterious state to be equal. Asymmetry in mutation rates would not qualitatively alter our conclusions as long as there is even weak migration (m 0 * 0:1), as load is then alleviated primarily by migration rather than reverse mutation. However, with zero migration and no reverse mutation, selection must be strong enough to prevent long-term ratchet-like accumulation of deleterious variants, resulting in a much higher threshold Ks c for metastability.
Second, we assume deleterious allele frequencies on the mainland to be close to deterministic. This assumption is not crucial: our qualitative conclusions remain unaltered as long as mainland frequencies are much lower than typical island frequencies. However, if mainland populations are small enough to harbour deleterious alleles at high frequencies at a subset of loci (which would, in general, be different from the loci fixed for deleterious alleles on the island), then we expect heterosis and the beneficial effects of migration to be weaker (with one-way migration between the mainland and island). More generally, extending this analysis to the co-evolution of load and population sizes in a metapopulation, where each sub-population may be royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210010 close to fixation for different deleterious alleles, remains an interesting direction for future work.
Third, we assume a rather simple genetic architecture of load: loci are assumed to be unlinked, and load additive across loci. Deviations from additivity, e.g. synergistic epistasis between deleterious variants can lower load [27,28], and may arise, for example, if multiple traits are under stabilizing selection. However, linkage between deleterious variants can inflate load by compromising selection efficacy at individual loci via Hill-Robertson interference [29], making it difficult to arrive at general predictions for the effects of selective interference (owing to linkage and epistasis between selected variants) on the eco-evolutionary dynamics of marginal populations. Fourth, we ignore environmental stochasticity and demographic Allee effects, which may strongly influence outcomes [8,30], especially in parameter regimes where mutation accumulation and demographic stochasticity in themselves are unlikely to cause extinction.
Finally, our analysis relies on a semi-deterministic analysis, which accounts for genetic drift but neglects demographic stochasticity. While this approximation captures qualitatively different population states across parameter space, it gives little insight into the dynamics. In particular, where alternative, i.e. low-load and high-load states are possible, the key assumption underlying the semi-deterministic analysis-namely, that allele frequencies have sufficient time to equilibrate at any given population size-is satisfied only when populations are in one or other state, and not while they transition between states. This makes it challenging to describe the complex coevolution of load and population size during transitions and arrive at a complete understanding of transition timescales.
Our aim has been to base our analysis on as few parameters as possible, in the hope that these can be related to observations from nature. We have reasonably good estimates of the fitness effects of deleterious mutations, and their degree of dominance-albeit largely from Drosophila [24]. We also now have accurate measures of the total mutation rate; the total rate of deleterious mutations is still uncertain [31], but may be substantial in complex organisms [32]. Population structure is less well understood: we have very many estimates of F ST [26], which reflects the numbers of incoming migrants, but local effective deme size is harder to estimate, even if demes can be defined at all. However, the common observation of heterosis implies that different deleterious recessives are common in different populations, suggesting a substantial drift load.
The rather complex effects of migration observed even in this relatively simple model with unconditionally deleterious alleles suggest that a comprehensive understanding of the effects of gene flow on eco-evolutionary dynamics at range limits must account for both environment-dependent (local) and environment-independent (global) components of fitness. These may be influenced by ( partially) overlapping sets of genetic variants, so that genetic load is shaped fundamentally by pleiotropic constraints. Such extended models are key to understanding when, for example, assisted gene flow is beneficial, and whether its mitigatory effect on inbreeding depression may be outweighed by any outbreeding depression that it might generate.
From a conceptual viewpoint, our analysis highlights the importance of considering explicit population dynamics when analysing the influence of gene flow on the efficacy of selection in sub-divided populations. Simple predictions, e.g. that a sub-divided population under hard selection should behave as a single population with an inbreeding coefficient equal to F ST [18,33], may break down when subpopulations can undergo local extinction. In this case, purging may be ineffective and the efficacy of selection reduced relative to undivided populations, in contrast to standard predictions that subdivision always reduces load when selection is hard. We regard the framework proposed here as a starting point for detailed studies in specific metapopulations, which take into account the joint evolution of population size and mutation load.