Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression

Abstract Human‐driven habitat fragmentation and loss have led to a proliferation of small and isolated plant and animal populations with high risk of extinction. One of the main threats to extinction in these populations is inbreeding depression, which is primarily caused by recessive deleterious mutations becoming homozygous due to inbreeding. The typical approach for managing these populations is to maintain high genetic diversity, increasingly by translocating individuals from large populations to initiate a “genetic rescue.” However, the limitations of this approach have recently been highlighted by the demise of the gray wolf population on Isle Royale, which declined to the brink of extinction soon after the arrival of a migrant from the large mainland wolf population. Here, we use a novel population genetic simulation framework to investigate the role of genetic diversity, deleterious variation, and demographic history in mediating extinction risk due to inbreeding depression in small populations. We show that, under realistic models of dominance, large populations harbor high levels of recessive strongly deleterious variation due to these mutations being hidden from selection in the heterozygous state. As a result, when large populations contract, they experience a substantially elevated risk of extinction after these strongly deleterious mutations are exposed by inbreeding. Moreover, we demonstrate that, although genetic rescue is broadly effective as a means to reduce extinction risk, its effectiveness can be greatly increased by drawing migrants from small or moderate‐sized source populations rather than large source populations due to smaller populations harboring lower levels of recessive strongly deleterious variation. Our findings challenge the traditional conservation paradigm that focuses on maximizing genetic diversity in small populations in favor of a view that emphasizes minimizing strongly deleterious variation. These insights have important implications for managing small and isolated populations in the increasingly fragmented landscape of the Anthropocene.


Impact Summary
Numerous threats to extinction exist for small populations, including the detrimental effects of inbreeding. Although much of the focus in reducing these harmful effects in small populations has been on maintaining high genetic diversity, here we use simulations to demonstrate that emphasis should instead be placed on minimizing strongly deleterious variation. More specifically, we show that historically large populations with high levels of genetic diversity also harbor elevated levels of recessive strongly deleterious mutations hidden in the heterozygous state. Thus, when these populations contract, inbreeding can expose these strongly deleterious mutations as homozygous and lead to severe inbreeding depression and rapid extinction. Moreover, we demonstrate that, although translocating individuals to these small populations to perform a "genetic rescue" is broadly beneficial, the long-term effectiveness of this strategy can be greatly increased by targeting moderate-sized source populations where recessive strongly deleterious mutations have been purged. These results challenge long-standing views on how to best conserve small and isolated populations facing the threat of inbreeding depression, and have immediate implications for preserving biodiversity in the increasingly fragmented landscape of the Anthropocene.
The prevailing paradigm in conservation biology prioritizes the maintenance of high genetic diversity in small populations threatened with extinction due to inbreeding depression (Caughley 1994;Spielman et al. 2004). Under this paradigm, genetic diversity is considered one of the primary determinants of fitness (Allendorf and Leary 1986;Reed and Frankham 2003), and the harmful effects of inbreeding are thought to be minimized by maintaining genetic diversity. However, this thinking is challenged by the observation that some species, such as the Channel island fox, can persist at small population size with extremely low genetic diversity and show no signs of inbreeding depression (Robinson et al. 2016(Robinson et al. , 2018. This example and other similar studies suggest that, rather than being mediated by high genetic diversity, the persistence of small populations may instead be enabled by the purging of strongly deleterious mutations (Laws and Jamieson 2011;Xue et al. 2015;Hedrick and Garcia-Dorado 2016;Robinson et al. 2016Robinson et al. , 2018Van Der Valk et al. 2019;Grossen et al. 2020). In this study, we investigate how genetic diversity, deleterious variation, and demographic history influence extinction risk due to inbreeding depression using ecologically realistic population genetic simulations. Our results demonstrate the central role of recessive strongly deleterious variation in determining extinction risk due to inbreeding depression in small and isolated populations, and highlight the counterintuitive effects of strategies aimed at maintaining high genetic diversity. We argue that, in cases where populations are destined to remain small and isolated with high levels of inbreeding, management strategies should aim to minimize strongly deleterious variation rather than maximize genetic diversity.
The motivating example for our simulations is the gray wolf population on Isle Royale, an island in Lake Superior that has long served as a natural laboratory in ecology and conservation biology (Mech 1966;Peterson et al. 1984;Wayne et al. 1991;McLaren and Peterson 1994). Following 70 years of nearcomplete isolation at a size of approximately 10-50 individuals, the population was driven to the brink of extinction by severe inbreeding depression, with just two individuals remaining in 2018 (Hedrick et al. 2019;Robinson et al. 2019). Recent findings have suggested that the collapse of the population may have been prompted by the introduction of high levels of recessive strongly deleterious variation by a migrant individual who arrived from the mainland in 1997 Hedrick et al. 2014Hedrick et al. , 2019Robinson et al. 2019). The high reproductive output of this individual on the island (34 offspring) led to intensive subsequent inbreeding in the population, driving these strongly deleterious mutations to high frequency and leading to severe inbreeding depression and ultimately the collapse of the population.
The example of the Isle Royale wolf population, although extreme, highlighted the potentially negative effects of founding or rescuing small populations with individuals from large and genetically diverse populations. An alternative approach for genetic rescue or reintroduction initiatives could instead target historically smaller populations where strongly deleterious mutations have been purged, or screen populations for individuals with low levels of strongly deleterious variation. The growing evidence for purging in wild populations (e.g., Xue et al. 2015;Robinson et al. 2018;Grossen et al. 2020) suggests that this approach may be effective as a means to decrease the severity of inbreeding depression in small populations at high risk of extinction. Given the ongoing reintroduction of wolves to Isle Royale, and the increasing necessity of reintroduction and genetic rescue more broadly (Whiteley et al. 2015;Frankham et al. 2017;Bell et al. 2019), such an approach could have wide-ranging implications for the conservation of small populations at risk.
The applicability of population genetic models to understanding extinction has historically been limited by unrealistic assumptions that often ignore stochastic ecological factors and rarely consider both weakly and strongly deleterious mutations (Lande 1994;Lynch et al. 1995;O'Grady et al. 2006;Caballero et al. 2017). Here, we use a novel population genetic simulation framework that combines realistic models of population dynamics with exome-scale genetic variation (Haller and Messer 2019) to assess how genetic diversity, deleterious variation, and demographic history influence extinction risk in small populations. Our simulations aim to capture the ecological factors that may contribute to extinction in small populations, such as those observed in the Isle Royale wolf population, by incorporating the effects of demographic and environmental stochasticity, as well as natural catastrophes. Coupled with these stochastic population dynamics, we model a genome with parameters reflecting that of a canine exome, which accumulates neutral and deleterious mutations from an empirically estimated distribution of fitness effects (DFE) ). Although our model is motivated by the Isle Royale wolf population, it is also intended to capture the dynamics of many other classic examples of population decline, inbreeding depression, and genetic rescue, such as Scandinavian wolves (Åkesson et al. 2016), Florida panthers (Johnson et al. 2010), and bighorn sheep (Hogg et al. 2006). Here, we focus on rapid contractions from large historical populations to very small populations, as these populations experience especially high risk of extinction due to inbreeding depression.

MODEL
We conducted non-Wright-Fisher (nonWF) simulations using SLiM 3 (Haller and Messer 2019). The impetus for this model was to allow for more ecologically realistic population genetic simulations by relaxing many of the restrictive assumptions of the Wright-Fisher model (Haller and Messer 2019). Such assumptions include nonoverlapping generations and a fixed population size that is not influenced by fitness, both of which are unrealistic when trying to model the extinction of a population due to genetic factors.
Instead, the SLiM nonWF approach models population size (N) as an emergent property of individual absolute (rather than relative) fitness and a user-defined carrying capacity (K). Thus, if individual fitness declines, a population experiences extinction through a biologically realistic process (a fitnessdriven reduction in population size). Further, the model includes overlapping generations, such that individuals with high fitness can survive and reproduce for multiple generations. At the start of each generation, each individual randomly mates with another individual in the population, with one offspring resulting from each mating. At the end of each generation, individuals die with a probability given by their absolute fitness (ranging from 0 to 1), which is rescaled by the ratio of K/N to model the effects of density dependence. Thus, the carrying capacity here does not directly determine the simulated population size, but rather it indirectly influences it through its impact on individual fitness and viability selection. For the sake of tractability, we assume a hermaphroditic random mating population. A discussion of how the carrying capacity of a SLiM nonWF model is related to its Wright-Fisher effective population size is provided in the Supplementary Methods.

DEMOGRAPHIC SCENARIOS
To obtain a baseline understanding for how ancestral demography influences extinction risk in small populations, we first explored a "population contraction" scenario with our simulations (Fig. 1A). Here, we modeled an instantaneous contraction from a large "ancestral population" (K ancestral ∈ {1000, 5000, 10,000, 15,000}) to a small "endangered population" (K endangered ∈ {25, 50}) (Fig. 1A). This contraction scenario could represent the isolation of a population with historical connectivity (e.g., the Florida panther population) or alternatively the founding of an isolated population through migration or translocation (e.g., the initial founding of the Isle Royale wolf population). For each contraction event, we randomly sampled K endangered number of individuals from the ancestral population to seed the endangered population after a burn-in of 10 * K ancestral generations. All simulations were run until the endangered population went extinct. To examine the effects of a more gradual contraction, we also explored a scenario where an ancestral population with carrying capacity 10,000 first contracted to an intermediate carrying capacity of 1000 for 200 generations, and finally an endangered carrying capacity of 25. We ran 25 simulation replicates for each combination of ancestral and endangered carrying capacities.
We next explored a "genetic rescue" scenario, which similarly consisted of a population contraction from a large ancestral population to a small endangered population (Fig. 1B).
Here, however, we fixed the ancestral carrying capacity to 10,000, and again explored two endangered carrying capacities of 25 and 50. Prior to the contraction, we split off the following source populations for genetic rescue: (1) a large population remaining at the ancestral size (K = 10,000); (2) a moderate-sized population with long-term isolation (K = 1000 for 1000 generations); (3) a small population with relatively recent isolation (K = 100 for 100 generations); and (4) a very small population with very recent isolation (K = 25 for 10 generations). These source population demographic histories were set to reflect a range of biologically relevant scenarios (i.e., large and outbred populations, populations with moderate size and long-term isolation, and smaller populations that have been more recently isolated) as well as provide a range of source population levels of genetic diversity and strongly deleterious variation to examine how these factors influence the efficacy of genetic rescue. Genetic rescue was initiated by translocating five randomly selected individuals from the source population after the endangered population decreased in size to five or fewer individuals for the case when K endangered = 25, and 15 or fewer individuals for the case when K endangered = 50. Importantly, the exact number of generations of isolation for these source populations depended on the number of generations before translocation, which varied for each simulation replicate depending on the stochastic trajectory of the endangered population. For these simulations, we ran 50 replicates for simulations with K endangered = 25 and 25 replicates for simulations with K endangered = 50.
To further explore the dynamics of our genetic rescue scenario, we ran several additional simulations, here with endangered carrying capacity fixed to 25. First, to investigate the impact of selecting individuals for translocation to either maximize genetic diversity or minimize strongly deleterious variation, we ran simulations (50 replicates) where we picked selected five migrants from the K = 10,000 source population with either the highest heterozygosity or fewest number of strongly deleterious alleles (where strongly deleterious mutations were defined to be those mutations with s < −0.01). We also explored scenarios with varying numbers of migrants (1, 5, or 10) with the number of translocations fixed to one, as well as varying the number of translocations (1, 2, or 5) with the number of migrants fixed to five. Here, we ran 25 replicates for each parameter combination. All simulations were run until the endangered population went extinct.

STOCHASTIC POPULATION DYNAMICS
To capture the nongenetic factors that can contribute to extinction in small populations (Caughley 1994), our model includes three sources of ecological stochasticity. First, demographic stochasticity was modeled using the built-in mechanics of the SLiM nonWF model, in which survival from one generation to the next is determined by a Bernoulli trial with the probability of survival determined by the absolute fitness of an individual scaled by K/N (Haller and Messer 2019). Next, we incorporated the effects of environmental stochasticity in our simulations by modeling the carrying capacity of the endangered population as an Ornstein-Uhlenbeck process, in which the carrying capacity in a generation is given by where φ = 0.9, K endangered ∈ {25,50}, and σ = log 10 (1.3) (Fig. S1). The values of φ and σ were set arbitrarily to model environmental stochasticity with a moderate amount of variation and a high degree of autocorrelation. Finally, we modeled the effects of random natural catastrophes in our simulations by drawing a probability of mortality due to a catastrophe each generation from a beta distribution with α = 0.5 and β = 8 (Fig. S1). In each generation, deaths due to a catastrophe are then determined by the outcome of a Bernoulli trial for each individual with the probability given by the beta distribution. Environmental stochasticity and natural catastrophes were only modeled in the small endangered population. Importantly, these stochastic ecological effects rarely lead to extinction in the endangered population in the absence of deleterious variation (8/25 simulation replicates with neutral mutations extinct within 10,000 generations at K endangered = 25, and none extinct for K endangered = 50).

GENOMIC PARAMETERS
We set the genomic parameters in our simulations to model the exome of a wolf-like organism. To do this, each diploid individual in our simulation has 20,000 genes of length 1500 bp, which occur on 38 autosomes, with the number of genes on each chromosome determined by the ratios observed in the dog genome (Lindblad-Toh et al. 2005). Recombination between these genes occurs at a rate of 1 × 10 −3 , with no recombination within genes and free recombination between chromosomes. These genes accumulate neutral and deleterious mutations at a rate of 1 × 10 −8 per site, with the ratio of deleterious to neutral mutations set to 2.31:1 Kim et al. 2017). The selection coefficients for these deleterious mutations were drawn from a DFE estimated from a large sample of humans ; see Methods in the Supporting Information for additional details).
We initially set the dominance coefficients for our simulations according to the following hs relationship from Henn et al. 2016, based on results from Agrawal and Whitlock 2011: This hs relationship intends to capture the pattern evident in empirical data that the dominance coefficient of a mutation tends to be inversely related to its selection coefficient, such that the most strongly deleterious mutations are highly recessive (Simmons and Crow 1977;Agrawal and Whitlock 2011;Huber et al. 2018). However, we found that simulations with this hs relationship in SLiM were extremely computationally intensive, such that we were only able to obtain results for simulations with the smallest ancestral carrying capacities of 1000 and 5000 (see Methods in the Supporting Information for further discussion).
To overcome these computational limitations for realistic models of dominance, we instead implemented an approach assuming that weakly/moderately deleterious mutations (s ≥ −0.01) were partially recessive (h = 0.25) and strongly deleterious mutations (s < −0.01) were fully recessive (h = 0). The aim of this approach (hereafter referred to as our "hmix" model) was to capture the key feature of the hs relationship that more deleterious mutations tend to be more recessive, with the dominance coefficients for these two classes of mutations reflecting their mean dominance coefficient under the hs relationship (Fig. S2). Given our finding that the behavior of this model is extremely similar to that of the hs relationship (see Results), we used the hmix model for all simulations except where otherwise noted. More details on how the hmix model was implemented in SLiM are provided in Methods in the Supporting Information.
To further explore the impact of dominance coefficients, we also ran simulations where we varied the dominance coefficient for all mutations (h ∈ {0, 0.01, 0.05, 0.2, 0.5}). In addition, we explored the impact of decreasing the number of sites in which fully recessive (h = 0) deleterious mutations can occur (i.e., the mutational target size) by varying the number of genes in our simulations (g ∈ {1000, 5000, 10,000, 15,000, 20,000}). These simulations were run solely under the "population contraction" demographic scenario with K ancestral ∈ {1000, 5000, 10,000, 15,000} and K endangered = 25. We ran 25 replicates for each of these simulations, terminating the simulation when the endangered population went extinct.
During the simulations, we kept track of several summaries of the state of the population. These include the population's mean heterozygosity, mean inbreeding coefficient (F ROH ), the mean fitness, and the average number of deleterious alleles per individual binned into weakly (−0.001 < s ≤ −0.00001), moderately (−0.01 < s ≤ −0.001), strongly (s < −0.01), and very strongly (s < −0.05) deleterious classes. Fitness was calculated multiplicatively across sites. Here, we restricted F ROH to include only runs of homozygosity greater than 1 Mb to capture inbreeding due to mating between close relatives. These statistics were calculated from a sample of 30 individuals every 1000 generations during the burn-in and every generation following the contraction.

BURN-IN CONDITIONS FOR THE SIMULATIONS
Our simulations during the burn-in retained fixed mutations and did not model reverse mutation. Retaining fixed mutations during the burn-in was important to ensure that weakly deleterious mutations that drifted to fixation contribute to absolute fitness. However, one consequence of retaining fixed mutations is that there is no mutation-selection-drift equilibrium in the ancestral population because weakly deleterious mutations continue to accumulate as fixations even after heterozygosity of segregating variation has leveled off (Fig. S3). As a result, fitness during the burn-in also never reaches equilibrium, but instead declines gradually as weakly deleterious mutations are fixed. Although fixation of weakly deleterious mutations occurs at a rate that is inverse to population size (i.e., much faster in smaller populations), we found that this effect is cancelled out when the length of the burn-in is proportional to the population's carrying capacity (i.e., 10 * K, leading to a much shorter burn-in for smaller populations), resulting in the same precontraction fitness regardless of population size (Fig. S4A). To explore the impact of these proportional burn-ins, we ran simulations under an hmix model of dominance in which all burn-ins were run for 30,000 generations, regardless of the ancestral size. This led to a notable reduction in precontraction fitness for the K = 1000 population, and a slight increase in fitness for the larger populations (Fig. S4). However, there were no qualitative differences in extinction times relative to when proportional burn-ins were used (Fig. S4C), suggesting that extinction dynamics are governed more by the numbers of recessive strongly deleterious alleles in these populations rather than their prebottleneck fitness. This finding is further supported by our simulation results that demonstrate that strongly deleterious mutations are a far more important determinant of extinction times compared to weakly or moderately deleterious mutations (see Results).

POPULATION CONTRACTION SIMULATIONS UNDER THE hm ix MODEL
Our population contraction simulations demonstrate that, although larger populations have higher genetic diversity ( Fig. 2A), they also harbor higher levels of strongly deleterious variation (Fig. 2B). Consequently, we observe a strong effect of ancestral demography on time to extinction following a population contraction, with populations that were historically large experiencing more rapid extinction (Fig. 2C). For example, given an endangered carrying capacity of 25, a population with an ancestral carrying capacity of 1000 will go extinct in 474 generations on average, whereas a population with an ancestral carrying capacity of 15,000 will do so in 70 generations on average (Fig. 2C). When modeling a more gradual contraction from an ancestral carrying capacity of 10,000, we found that extinction times were slightly increased relative to the instantaneous contraction scenario (Fig. S5), suggesting that more gradual contractions can facilitate purging and ultimately decrease extinction risk.
Examination of individual simulation replicates provides insight into the dynamics of extinction in these populations (Figs. 3A and S6). Endangered populations with an ancestral carrying capacity of 1000 contain fewer strongly deleterious mutations, translating to a decreased severity of inbreeding depression and ultimately longer persistence (Figs. 3A and S6). This decreased severity of inbreeding depression allows these populations to become highly inbred (F ROH ≈ 1) well before going extinct (Figs. 3A and S6). By contrast, endangered populations with a larger ancestral carrying capacity of 15,000 have much higher levels of recessive strongly deleterious variation due to these mutations being hidden from selection in the ancestral population, resulting in a more rapid loss of fitness as these populations become inbred (Figs. 3B and S6). As a result, these populations typically go extinct well before F ROH approaches one (Figs. 3B and S6). These populations also lose fitness due to increased ge-netic drift and inbreeding among more distant relatives, which is not captured by our definition of F ROH .
Following contraction, the ability of the endangered population to purge its load of recessive deleterious mutations also depended on stochastic ecological factors. For example, when the carrying capacity of the endangered population was by chance higher due to environmental stochasticity, natural selection was most effective at removing strongly deleterious alleles, translating to longer persistence (Figs. 3 and S6). By contrast, when the carrying capacity of the endangered population was low soon after contraction, purging tended to be less effective, resulting in more rapid extinction (Figs. 3 and S6). However, in both cases, purging was also counteracted by continual input of strongly deleterious alleles by mutation, which eventually contributed to population extinction. Overall, these various sources of genetic and ecological stochasticity together resulted in highly variable extinction times for any given parameter settings, highlighting the important role of random events in determining whether a small population can persist.
Our simulations also demonstrate the importance of increasing the carrying capacity of the endangered population as a means to ensure population persistence. Larger endangered populations (K = 50) were better able to purge recessive deleterious mutations following the contraction, resulting in much longer persistence (Fig. S7). Moreover, larger populations were less impacted by stochastic ecological factors, which also contributed to increased time to extinction. Nevertheless, extinction times for these larger endangered populations still depended strongly on the ancestral carrying capacity (Fig. S7), demonstrating the importance of considering both recent and historical demography when assessing extinction risk.

RESULTS UNDER VARYING GENOMIC PARAMETERS
We next explored the sensitivity of our results to the genomic parameters assumed in our model. In particular, we focus on the impact of dominance coefficients, as previous work has suggested that the extent to which strongly deleterious mutations accumulate at higher frequencies in larger populations depends strongly on the assumed dominance coefficients (Nei 1968;Hedrick 2002;Hedrick and Garcia-Dorado 2016). Under the most realistic model of an hs relationship, we found that strongly deleterious mutations do accumulate at much higher frequencies in larger populations, leading to much faster extinction following contraction (Fig. 4). However, we also found that simulations with an hs relationship in SLiM were extremely computationally intensive, such that we were unable to obtain results for the larger ancestral carrying capacities of 10,000 and 15,000 (see Supporting Information). Nevertheless, these results demonstrate a strong concordance between the hs relationship and our hmix model, suggesting that our hmix model captures the key features of the hs relationship (Figs. 4 and S2).
When assuming a single dominance coefficient for all mutations, our results demonstrate a range of outcomes that depended on the assumed dominance coefficient. Specifically, when assuming h = 0 or 0.01, we again find that larger populations harbor higher levels of strongly deleterious variation, with similar patterns to the hs relationship or hmix model (Fig. 4). However, this effect is greatly diminished when mutations are only partially recessive (h = 0.05 or 0.2), and is nonexistent when mutations are additive (h = 0.5), as expected (Fig. 4). Notably, although the average dominance coefficient under our assumed hs relationship is approximately 0.2, our results with h = 0.2 are dramatically different from those under the hs relationship (Fig. 4).
To further investigate the importance of strongly versus weakly/moderately deleterious mutations in our hmix model, we ran simulations where we truncated our DFE to only permit strongly deleterious (s < −0.01 and h = 0) or weakly/moderately deleterious mutations (s ≥ −0.01 and h = 0.25) to enter the population. More specifically, in the case of permitting only strongly deleterious mutations, we allowed any mutation with s < −0.01 to enter the population as normal; however, mutations with s ≥ −0.01 instead became neutral mutations. Here, our results for simulations that only included strongly deleterious mutations were notably similar to the full hmix model, with extinction times that depended strongly on ancestral demography (Fig. S8). By contrast, results for simulations that only included weakly/moderately deleterious mutations were dramatically different from the full model, and had greatly increased extinction times (Fig. S8). Given that strongly deleterious muta- tions constitute only approximately 25% of all new deleterious mutations under our assumed DFE , these results demonstrate their disproportionate impact on extinction risk relative to the effects of more weakly or moderately deleterious mutations.
As a final way of exploring the impact of recessive deleterious mutations, we also conducted simulations where we decreased the target size for deleterious mutations assuming h = 0. The motivation for these simulations was to test whether we would still observe an impact of ancestral demography on extinction times when lowering the number of genes that could accumulate fully recessive deleterious mutations. To do this, we decreased the number of genes (g) in our simulations from 20,000 to g ∈ {1000, 5000, 10,000, and 15,000}. Here, we observed that the effect of ancestral demography is still present, although greatly diminished, with as few as 1000 genes, and remains substantial with as few as 5000 genes (Fig. S9).

GENETIC RESCUE SIMULATIONS
We next explored how demography and strongly deleterious variation impact the effectiveness of genetic rescue assuming the hmix model of dominance. Here, we quantify the effectiveness of genetic rescue as the extent to which the introduction of migrants to the endangered population increased extinction times. When translocating five migrants from one of four source populations (Figs. 1B and 5A-B) to an endangered population with K = 25, we found that all source populations led to increases in time to extinction relative to a no-rescue scenario (Fig. 5C). Importantly, however, we found that the magnitude of increase in time to extinction was highly dependent on source population demography and levels of strongly deleterious variation. For example, although genetic rescue from the large source population (K = 10,000) led to a notable increase in mean time to extinction of 16% (one-tailed t-test P = 0.159), rescue from the moderate-sized source population resulted in a much more dramatic increase of 130% (P = 1.75 × 10 -7 ; Fig. 5C). Genetic rescue from small and moderately inbred populations (mean F ROH < 0.1; Fig. S10) also resulted in increases in mean time to extinction that were at least as great as that of the large source population (63% increase for K = 100 [P = 0.00132]; 13% increase for K = 25 [P = 0.224]) (Fig. 5C). We observed similar patterns when conducting simulations with a larger endangered population (K = 50) (Fig. S11), although the beneficial effects of genetic rescue were somewhat diminished, likely due to the larger recipient population size when translocation was initiated (N ≤ 15) as well as the greater efficacy of purging within the larger recipient population.
Examination of individual simulation replicates again offers insight into the factors driving extinction in these populations (Figs. 6 and S12). In all simulated scenarios, we observed a large increase in fitness immediately following the introduction of migrants (Figs. 6, S12, and 13), likely due to the high strength of heterosis (masking of fixed recessive deleterious alleles) following initial admixture. In most cases, these increases in fitness led to population growth, although the extent to which population sizes increased was strongly influenced by environmental stochasticity (Figs. 6 and S12). Soon after genetic rescue, however, the resumption of inbreeding again led to a decline in fitness (Figs. 6 and S12). Here, the rate of fitness decline was determined by the levels of strongly deleterious variation in the migrant individuals and their descendants. For example, when translocation occurred from a large source population (K = 10,000), levels of strongly deleterious variation remained high after genetic rescue, eventually resulting in severe inbreeding depression once inbreeding resumed (Fig. 6). By contrast, when the moderate-sized source population was used (K = 1000), levels of strongly deleterious variation dramatically decreased following genetic rescue, greatly decreasing the severity of inbreeding depression in future generations (Fig. 6).

(B) Example when a moderate-sized (K = 1000) is used. The top of each panel shows population size (N) in black and carrying capacity (K) in gold, the middle shows the average number of strongly deleterious alleles (s < −0.01) per individual, and the bottom shows mean absolute fitness in blue and mean inbreeding coefficient (F ROH ) in green. The dashed gray line indicates the generation in which migration occurred. For both replicates, K endangered = 25 and an hmix model of dominance was assumed. Note the differing x-axis scales for panels A and B.
These results demonstrate that, although the larger source populations have higher fitness and genetic diversity, they also harbor a high number of heterozygous recessive strongly deleterious mutations (Figs. 5 and S10). These mutations are quickly made homozygous by inbreeding in the endangered population following translocation, exacerbating the severity of inbreeding depression and eventually contributing to extinction. In support of this interpretation, we found that time to extinction following genetic rescue is predicted by the average number of strongly deleterious alleles per individual in the source population when examining all source populations simultaneously (Fig. 5D). By contrast, we found that average source population heterozygosity was not very predictive of time to extinction, where in fact we observed a slight negative correlation due to the fact that source populations with higher genetic diversity also tend to harbor higher levels of recessive strongly deleterious variation (Fig. 5E). We obtained similar results when restricting this analysis to the K = 25 source population, the source population for which there was the highest variability in strongly deleterious variation and genetic diversity (Fig. S14). However, in this case, we do see a strong negative correlation between genetic diversity and extinction times, which is likely driven by heterozygosity being correlated with strongly deleterious variation at this intrapopulation scale.
Our finding that source population levels of strongly deleterious variation predict the efficacy of genetic rescue suggests that genome-wide levels of deleterious variation could be used to select individuals for genetic rescue. We explored the efficacy of this strategy by selecting the individuals with the fewest strongly deleterious alleles (s < −0.01) from the large source population (K = 10,000) for rescue. This approach resulted in an increase in mean time to extinction of 74% (P = 1.15 × 10 -4 ) compared to the nonrescue scenario, a 49% increase relative to randomly selecting individuals from the large source population (Fig. 5F). By contrast, when individuals with the highest genome-wide heterozygosity were selected, we observed an increase in time to extinction that was no greater than when individuals were selected at random (P = 0.148; Fig. 5F). These results further support the causal relationship between strongly deleterious variation and extinction risk due to inbreeding depression, as well as the lack of this relationship between heterozygosity and extinction risk.
Last, we explored the effects of varying the number of migrants (1, 5, or 10) as well as the number of translocation events (1, 2, or 5), focusing on the K = 10,000 and K = 1000 source populations. These simulations show persistent increases in time to extinction with increasing number of translocations regardless of the source population (Fig. S15), suggesting that the efficacy of genetic rescue does not diminish with each additional migration event. When varying the number of migrants, we found that extinction times increased slightly in the case where migrants were drawn from the K = 1000 source population, but did not increase when drawing migrants from the K = 10,000 source population (Fig. S15).

Discussion
Our simulations demonstrate the central importance of demographic history and recessive strongly deleterious mutations in determining the risk of extinction due to inbreeding depression in small and isolated populations. We find that populations that were historically large have a much higher risk of extinction following a population contraction compared to historically smaller populations (Fig. 2). These findings may be counterintuitive given the thinking that small populations should be less fit due to an accumulation of weakly deleterious alleles (Kimura et al. 1963;Lynch and Gabriel 1990;Lynch et al. 1995;Battaillon and Kirkpatrick 2000). The key insight that our simulations highlight is that larger ancestral populations harbor more recessive strongly deleterious mutations (Fig. 2), due to these mutations being hidden from purifying selection as heterozygotes. The exposure of these mutations as homozygous by inbreeding in small populations can lead to dramatic reductions in fitness and drive rapid extinction, well before "mutational meltdown" due to an accumulation of weakly deleterious mutations can occur (Lynch and Gabriel 1990;Lynch et al. 1995). By demonstrating that this effect is sufficient to ultimately drive extinction on short timescales, our simulations illustrate the importance of recessive deleterious variation as the primary mechanism of inbreeding depression (Charlesworth and Willis 2009;Hedrick and Garcia-Dorado 2016).
Although the insight that strongly deleterious mutations that are fully recessive can accumulate at greater frequencies in larger populations has been noted elsewhere (Nei 1968;Battaillon and Kirkpatrick 2000;Hedrick 2002;Hedrick and Garcia-Dorado 2016;Szpiech et al. 2019), our results add to this work by demonstrating that these effects persist under realistic genomic parameters and models of dominance, including both an hs relationship and our hmix model (Fig. 4). This result is due to the fact that strongly deleterious mutations under these models tend to be highly recessive (Simmons and Crow 1977;Agrawal and Whitlock 2011;Huber et al. 2018), and are therefore hidden from purifying selection when present as heterozygotes in large populations Moreover, when examining models with a single dominance coefficient for all mutations, we find that our results with fully recessive mutations (h = 0 or 0.01) are highly similar to those under more realistic models of dominance, whereas the impact of ancestral demography is greatly diminished with partially recessive (h = 0.05 or 0.2) mutations and nonexistent with additive (h = 0.5) mutations (Fig. 4). The transition in this behavior has previously been explored using analytical models, where it was similarly shown that the impact of population size on levels of strongly deleterious variation greatly diminishes as h approaches 0.05 (Hedrick 2002;Hedrick and Garcia-Dorado 2016) Although the mean dominance coefficient under our assumed hs relationship is ∼0.2, this model showed very different extinction dynamics compared to the model assuming h = 0.2 for all mutations (Fig. 4). This result suggests that the key factor mediating extinction dynamics is the presence of highly recessive strongly deleterious mutations, rather than the overall mean dominance coefficient. In support of this interpretation, we show that simulations that only include the fraction of mutations that are strongly deleterious and fully recessive (h = 0) exhibit extinction times that are highly similar to those under an hs relationship or hmix model, whereas simulations that only include the fraction of mutations that are weakly or moderately deleterious and partially recessive (h = 0.25) have much longer extinction times that depend only minimally on ancestral demography (Fig. S8). This result is especially striking when considering that strongly deleterious mutations comprise only ∼25% of the overall number of deleterious mutations under our assumed DFE . Taken together, these results further bolster our conclusion that strongly deleterious mutations are a key determinant of extinction risk due to inbreeding depression in small populations. The extent to which these dynamics may be relevant for a given species, however, will depend on their distribution of selection and dominance coefficients, and specifically the extent to which there exists a substantial fraction of mutations that are highly recessive and strongly deleterious. Although these parameters remain poorly known for the vast majority of organisms, highly recessive strongly deleterious mutations are commonly observed across disparate taxa (e.g., Simmons and Crow 1977;McCune et al. 2002), suggesting that their influence on extinction may be widespread. However, these dynamics may prove to be most relevant to mammals and other vertebrates with high levels of "organismal complexity," which has been demonstrated to result in a higher fraction of strongly deleterious mutations .
The considerable influence of ancestral demography that our simulations reveal has important implications for predicting the threat of extinction due to inbreeding depression in wild populations. Quantifying inbreeding depression in the wild and predicting the threat it poses to extinction remains a major challenge in conservation biology, and the reasons why some small populations suffer from inbreeding depression and others do not are often unclear (Hedrick and Garcia-Dorado 2016). Our simulations suggest that these differences may be determined in large part by the ancestral demography of a species. Consequently, we suggest that information on ancestral demography, which is increasingly becoming accessible using genomic data (Beichman et al. 2018), should be more widely incorporated into extinction risk predictions. In particular, given that nearly all threatened populations have recently undergone reductions in size due to anthropogenic pressures, these results suggest that their continued persistence may depend crucially on their ancestral demography and load of recessive strongly deleterious mutations. In addition, our results also suggest that island populations that have historically been small may in fact have reduced extinction risk due to inbreeding depression as a result of to their high isolation and smaller population size, which may have facilitated purging of recessive strongly deleterious mutations (Laws and Jamieson 2011;Robinson et al. 2018). However, our simulations also reveal that the fate of small populations is highly stochastic, and that even under the same ecological and genetic parameters, time to extinction can vary substantially (Figs. 2 and 4). Thus, predictions of extinction risk for any wild population should necessarily be accompanied by a high degree of uncertainty.
Our results also provide insight on how best to conduct genetic rescue, which is becoming increasingly necessary for maintaining small and isolated populations under growing anthropogenic pressures (Whiteley et al. 2015;Hedrick and Garcia-Dorado 2016;Frankham et al. 2017;Bell et al. 2019). Consistent with many other studies (e.g., Hogg et al. 2006;Johnson et al. 2010;Frankham 2015;Åkesson et al. 2016), our simulations highlight the beneficial effects of genetic rescue, supporting recent calls for its more widespread application (Whiteley et al. 2015;Frankham et al. 2017;Bell et al. 2019). However, in stark contrast to existing recommendations (Pickup et al. 2012;Whiteley et al. 2015), we found that targeting large source populations with high genetic diversity was among the least effective genetic rescue strategies (Figs. 5 and S11). Instead, our results demonstrate that the effectiveness of genetic rescue can potentially be maximized by drawing migrants from moderate-sized source populations (Figs. 5 and S11). These source populations are ideal due to being small enough to purge strongly deleterious recessive mutations, but not so small that they accumulate a substantial load of fixed weakly deleterious mutations. Finally, our results demonstrate that even small and somewhat inbred populations (mean F ROH < 0.1) can also serve as effective source populations (Figs. 5C and S10), as other studies have shown (Heber et al. 2012a,b). Our results do not suggest, however, that inbred populations will always serve as more effective source populations than outbred populations, or that highly inbred source populations are likely to be effective, as has been recently claimed (Ralls et al. 2020). Rather, our results demonstrate that source populations that are only slightly inbred (mean F ROH ∼0.06) may be comparable in their effectiveness to large source populations, although both are likely to be less effective than purged source populations with long-term moderate size.
A key factor accounting for the differences between our results and previous experimental and empirical studies examining source population effectiveness (Pickup et al. 2012;Frankham 2015) is a difference of temporal scale. Specifically, existing studies rarely examine the consequences of gene flow beyond the F1 generation, whereas our simulation framework enables tracking these dynamics over the longer term. Previous research has shown that, during the generations immediately following genetic rescue, fitness increases are likely to be greatest when migrants are sourced from larger source populations (Pickup et al. 2012;Frankham 2015). Indeed, our results agree with these findings, and demonstrate greater fitness increases during the F1 generation when drawing migrants from larger source populations (Fig.  S13). Where our results differ, however, is in the generations well after F1, after the initial strong heterosis is diminished and inbreeding resumes. Here, we find that using smaller source populations with decreased levels of strongly deleterious variation can lower the severity of inbreeding depression and increases population persistence (Fig. 6). This result emphasizes the importance of evaluating genetic rescue outcomes over the longer term, because the dynamics in the generations immediately following genetic rescue might differ greatly from those several generations later, as was so dramatically demonstrated by the Isle Royale wolf population.
Although our results suggest that large populations may not be ideal source populations for genetic rescue, we demonstrate that the effectiveness of this strategy can be greatly increased when individuals are screened for low levels of strongly deleterious variation (Fig. 5F). However, applications of such an approach may remain limited by our ability to accurately predict the fitness consequences of individual mutations, which remains a major challenge even in humans (Eilbeck et al. 2017). Our results also demonstrate that repeated genetic rescue from large or moderate-sized source populations can result in persistent beneficial effects (Fig. S15), highlighting the efficacy of this strategy for populations that are destined to remain small and isolated. One inevitable consequence of this or any other genetic rescue strategy, however, is a loss of native ancestry (Johnson et al. 2010;Adams et al. 2011;Harris et al. 2019), which can potentially swamp out locally adapted alleles (although see Fitzpatrick et al. 2020). Although we do not track levels of admixture in our simulations, it is probable that the post-rescue populations are composed of highly admixed individuals (Harris et al. 2019).
In addition to providing guidelines for how to best conduct genetic rescue, our results also have implications for understanding the mechanisms underlying genetic rescue. Two mechanisms have generally been proposed for genetic rescue: heterosis (masking of fixed recessive deleterious mutations) and adaptive evolution (increases in fitness due to selection on new alleles) (Whiteley et al. 2015). By demonstrating that, even in the absence of adaptive variation, migration into small and inbred populations can lead to increases in fitness and population size consistent with those observed in empirical systems (e.g., Hogg et al. 2006;Johnson et al. 2010;Adams et al. 2011;Åkesson et al. 2016), our results suggest that heterosis alone may be able to explain much of the beneficial effects of genetic rescue. However, a more thorough investigation of the relative roles of heterosis and adaptive evolution as mechanisms of genetic rescue could be conducted by incorporating adaptive mutations into our simulation framework. Finally, our finding that increases in fitness following migration do not always lead to increases in population size due to the effects of environmental stochasticity (Fig. S12) suggests that definitions of genetic rescue based on demographic effects alone may be restrictively narrow, as other authors have suggested . In particular, this result has relevance to the debate of whether the Isle Royale wolf population truly experienced genetic rescue, given that the population size did not increase substantially following migration, perhaps due to poor environmental conditions Hedrick et al. 2011) Our simulation framework has several notable limitations. First, due to the high computational load of our simulations (see Methods in the Supporting Information), we were unable to examine ancestral populations with carrying capacities larger than 15,000, limiting the observed impact of ancestral demography. Second, we did not examine heterozygote advantage as a potential mechanism for inbreeding depression, which could influence the negative relationship between extinction times and source population heterozygosity that we observed (Fig. 5E). However, empirical support for heterozygote advantage as a mechanism for inbreeding depression is scarce (Charlesworth and Willis 2009;Hedrick and Garcia-Dorado 2016), suggesting that it may not substantially impact our results. In addition, we did not include adaptive mutations in our model, which may also be relevant to the ability of small populations to persist under environmental change, as well as for exploring the impacts of outbreeding depression due to differential local adaptation. Although recent work has suggested that concerns about outbreeding depression may have been overstated (Whiteley et al. 2015;Frankham et al. 2017), assessing its influence in determining source population selection nevertheless represents an important area for future research. Finally, we emphasize that the demographic scenarios explored here may not be applicable to all small and isolated populations. Specifically, our assumption of an instantaneous contraction from a large ancestral population to a small endangered population may not capture the demographic trajectory of many populations that have experienced more gradual declines. In these cases, a gradual decline may have facilitated purging of strongly deleterious mutations (Fig. S5), decreasing the necessity for genetic rescue (e.g., Xue et al. 2015;Robinson et al. 2018). Overall, we recommend that any conservation actions motivated by our results carefully consider the demography of the population of interest, ideally using simulations.
In conclusion, our results highlight the detrimental effects of strongly deleterious variation in small populations, and suggest that conservation strategies for populations experiencing inbreeding depression may be improved by minimizing strongly deleterious variation rather than maximizing heterozygosity. Heterozygosity at putatively neutral markers (e.g., microsatellites and RAD loci) has long served as an essential tool for assessing population health and risk of extinction due to inbreeding depression (Teixeira and Huber 2020); here, we highlight ways in which this approach may be misleading. These results are especially relevant to the ongoing reintroduction of wolves to Isle Royale, which has been guided by the goal of maximizing the genetic diversity of the new population. Although our simulations are motivated by this example, by examining a wide range of demographic scenarios, our results have implications beyond the Isle Royale wolf population. Future research should explore implementation of these strategies by expanding the use of genomic tools and assessments of deleterious variation in the context of the conserving small and isolated populations . Given the great expense of most translocation programs, incorporating genomic approaches represents a sound investment with the potential to substantially postpone the need for future intervention.
Associate Editor: S. Wright

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Carrying capacity, population size, and effective population size Figure S1. Summary of stochastic population dynamics. Figure S2. Comparison of distribution of dominance coefficients under hs relationship from Henn et al. 2016 and hmix model. Figure S3. Burn-in dynamics for the K=5000 population assuming an hmix model of dominance. Figure S4. Comparison of extinction dynamics between population contraction simulations with burn-in duration of 10 * K ancestral (left panels) and simulations with burn-in duration of 30,000 generations, assuming an hmix model of dominance. Figure S5. Comparison of the time to extinction in an instantaneous vs gradual contraction from K ancestral =10,000 to K endangered =25 under an hmix model of dominance. Figure S6. Population trajectory of individual simulation replicates following contraction to endangered carrying capacity of 25 from ancestral carrying capacity of (A) 1,000, (B) 5,000, (C) 10,000, and (D) 15,000 under an hmix model of dominance. Figure S7. Time to extinction following a population contraction under an hmix model of dominance when K endangered =50. For each ancestral carrying capacity, 25 simulation replicates were run. Figure S8. Time to extinction following a contraction to K endangered = 25 comparing results under the full hmix and hs relationship models to those that only included the fraction of mutations in the hmix model that are strongly deleterious (s < −0.01 and h=0) or weakly/moderately deleterious (s ≥ −0.01 and h=0.25). Figure S9. Population contraction results when reducing the target size for deleterious mutations under a fully recessive model of dominance (h=0). Figure S10. Characteristics of the source populations used for genetic rescue. Figure S11. Time to extinction following genetic rescue from source populations of varying size with recipient population carrying capacity of 50 under an hmix model of dominance. Figure S12. Population trajectory of individual simulation replicates with K endangered = 25 and K ancestral = 10,000 following genetic rescue from source populations of varying demographic history: (A) K source = 10,000, (B) K source = 1,000 for 1,000 generations, (C) K source = 100 for 100 generations, and (D) K source = 25 for 10 generations. Figure S13. Absolute fitness increases during the generation immediately following the introduction of five migrants with K endangered = 25. Figure S14. Relationship between extinction times and source population deleterious variation and genetic diversity for the K=25 source population. Figure S15. Time to extinction following genetic rescue with varying numbers of translocation events and migrants with K endangered = 25.