Population demographic history and evolutionary rescue: Influence of a bottleneck event

Abstract Rapid environmental change presents a significant challenge to the persistence of natural populations. Rapid adaptation that increases population growth, enabling populations that declined following severe environmental change to grow and avoid extinction, is called evolutionary rescue. Numerous studies have shown that evolutionary rescue can indeed prevent extinction. Here, we extend those results by considering the demographic history of populations. To evaluate how demographic history influences evolutionary rescue, we created 80 populations of red flour beetle, Tribolium castaneum, with three classes of demographic history: diverse populations that did not experience a bottleneck, and populations that experienced either an intermediate or a strong bottleneck. We subjected these populations to a new and challenging environment for six discrete generations and tracked extinction and population size. Populations that did not experience a bottleneck in their demographic history avoided extinction entirely, while more than 20% of populations that experienced an intermediate or strong bottleneck went extinct. Similarly, among the extant populations at the end of the experiment, adaptation increased the growth rate in the novel environment the most for populations that had not experienced a bottleneck in their history. Taken together, these results highlight the importance of considering the demographic history of populations to make useful and effective conservation decisions and management strategies for populations experiencing environmental change that pushes them toward extinction.

environment and prevent an otherwise inevitable extinction is called evolutionary rescue (Carlson et al., 2014;Gomulkiewicz & Holt, 1995). This process is characterized by a U-shaped curve in population size through time: initially, population size decreases due to maladaptation to the environment, then population size increases as adaptive alleles spread through the population (Gomulkiewicz & Holt, 1995;Orr & Unckless, 2014).
Theory and experiments with model organisms show that evolutionary rescue can avert extinction. Theory has demonstrated that the probability of evolutionary rescue increases with genetic variation (Barrett & Schluter, 2008;Gomulkiewicz & Holt, 1995) and suggests that it decreases with mutational load and inbreeding depression (Barrett & Kohn, 1991;Ellstrand & Elam, 1993). Despite confirmation of these effects in experiments with model organisms (Agashe et al., 2011;Bijlsma et al., 2010;Hufbauer et al., 2015;Lachapelle & Bell, 2012;Ramsayer et al., 2013), evidence of evolutionary rescue in nature remains scarce, suggesting that the conditions that allow evolutionary rescue to occur may be limited (Carlson et al., 2014). One condition that is particularly important to examine experimentally is the demographic history of populations.
In the wild, populations are frequently subjected to various stressful events such as hunting or habitat destruction that can lead to a temporary reduction in size. These bottlenecks in population size can affect genetic variation, mutational load, and inbreeding levels in the populations that experience them (hereafter bottlenecked populations; Frankham et al., 1999). Such bottleneck events, whether rare extreme events or repeated less-severe restrictions in population size, can reduce the future adaptive potential of populations (Bouzat, 2010;Nei et al., 1975) and thus the probability of evolutionary rescue.
Bottlenecks influence the ability of populations to adapt to a future challenging environment in two main ways (Frankham et al., 1999;O'Brien & Evermann, 1988;Reusch & Wood, 2007).
First, standing genetic variation is lost by genetic drift during bottlenecks due to an increase in stochastic events associated with small population size (Hedrick, 2005). This reduction of genetic diversity decreases the evolutionary potential of populations since potential future adaptive alleles are lost (Futuyma, 1986). In general, adequate genome-wide genetic variation is important for the conservation of natural populations (Kardos et al., 2021). In addition, when populations are small, selection is weak relative to genetic drift. When this occurs, there can be an increase in the frequency and fixation of deleterious alleles (Lynch et al., 1995), leading to populations with lower growth rates (average number of offspring per individual) than populations that have not experienced a bottleneck. However, given the stochastic way in which alleles go to fixation or are lost, by chance, some beneficial alleles may become fixed, potentially leading some bottlenecked populations to have higher growth rates than nonbottlenecked populations (Bouzat, 2010). Thus, in addition to a lower growth rate, greater variation in growth rate is expected across different bottlenecked populations compared to nonbottleneck populations.
The second main way that bottlenecks can influence the ability of populations to adapt is through inbreeding. After bottleneck events, inbreeding increases due to increased homozygosity of alleles common by descent (Charlesworth & Willis, 2009;Fox et al., 2008;Frankham et al., 1999;Hedrick & Fredrickson, 2010;Kirkpatrick & Jarne, 2000;Ralls et al., 1979;Wright, 1977). When homozygosity increases at loci with deleterious alleles, inbreeding leads to inbreeding depression (i.e., a reduction in absolute fitness of inbred individuals). Inbreeding depression increases the risk of extinction for the population by reducing its growth rate (Charlesworth & Charlesworth, 1987;Frankham, 2005aFrankham, , 2005bO'Grady et al., 2006). However, inbreeding in populations can also lead to purging of deleterious alleles, removing them from the population and increasing growth rate of the population (Crnokrak & Barrett, 2002;Facon et al., 2011;Glémin, 2007). Thus, as for genetic drift, the effects of inbreeding due to bottlenecks can increase the variation in growth rate among populations.
Despite the strong evidence for a positive effect of genetic variation (Hedrick, 2005) and a negative effect of inbreeding depression (Frankham et al., 1999) on the rapid adaptation that is central to evolutionary rescue, the consequences of past bottleneck events on adaptation have only recently begun to be explored. Experiments using bacteria demonstrate that the adaptive pathways followed by populations under stress depend on whether bottlenecks are weak or intense (Garoff et al., 2020;Mahrt et al., 2021;Vogwill et al., 2016), supporting the idea that demographic history will influence the probability of evolutionary rescue. In those experiments, with large bacterial populations, de novo mutations generated genetic diversity. However, in natural diploid populations with sexual reproduction, adaptation relies mostly on standing genetic variation, which should result in very different adaptation dynamics. Furthermore, in a recent study, Klerks et al. (2019) evaluated how population demographic history (bottlenecked or not) influenced the response to selection for increased heat tolerance in experimental populations of a diploid, sexually reproducing eucaryote: the least killifish Heterandria formosa. They found that bottlenecked populations had a 50% slower response to selection than normal populations. This suggests that bottlenecks play an important role in the outcomes of evolutionary rescue. However, monitoring population dynamics was outside the scope of Klerks et al. (2019) who studied adaptation in experimental populations of fixed size, and so it is not clear if the U-shaped curve characteristic of evolutionary rescue would have been observed had population size been allowed to vary. Studying population dynamics is, however, necessary to understand the minimum population size that populations will reach without going extinct or how long it will take for populations of a given size to recover. Given the role of bottlenecks in reducing genetic variation and increasing inbreeding depression in diploid eukaryotes that are often of conservation concern (Frankham et al., 1999;Reusch & Wood, 2007), it is important to examine the effects of bottlenecks on evolutionary rescue experimentally.
Here, we studied how bottleneck events influence the process of evolutionary rescue of populations in a severely challenging environment. Because bottlenecks in nature vary in intensity, we created populations that had experienced no bottleneck in their demographic history, an intermediate bottleneck, or a strong bottleneck, and monitored population dynamics in a challenging environment.
By using Tribolium castaneum to model a diploid, semelparous organism, we evaluated how the probability of extinction depended upon population demographic history, and whether surviving populations exhibited the U-shaped population trajectory characteristic of evolutionary rescue, and how that depended upon population demographic history. Finally, we quantified adaptation to the new environment at the end of the experiment by measuring the growth rate of all rescued populations and evaluating the evolution of the rate of development, a life-history trait that could be involved in adaptation (Berner et al., 2004).

| Model system
We used T. castaneum, the red flour beetle, to test the effects of demographic history on evolutionary rescue. Tribolium castaneum populations were kept in 4 cm by 4 cm by 6 cm acrylic containers with 20 g of a standard medium (95% wheat flour and 5% brewer's yeast), which served as both habitat and a food source. We will henceforth refer to a container and medium together as a patch. Multiple patches were used to maintain each distinct laboratory population, and single patches were used in growth assays. During laboratory maintenance, the life cycle of T. castaneum was constrained to mimic that of a seasonally breeding organism with nonoverlapping generations and a discrete dispersal phase (Melbourne & Hastings, 2008), a life history found commonly among plants and animals. Each generation, adult beetles from all the patches were mixed. Then, we placed them in groups of 50 individuals per patch with fresh standard medium for 24 h to reproduce and start the next generation. At the end of the reproduction phase, adults were removed from the patches and the eggs were left to mature into adults for an additional 34 days for a total of 35 days per generation. During laboratory maintenance and the experiments, beetles were kept in incubators at 31°C with 54% ± 14% relative humidity. Several incubators were used, and patches were randomized among and within incubators once each week to prevent systematic effects of incubation conditions.

| Creating populations with different demographic histories
A large diverse population was initiated from the multiple crosses of different stock strains ("GA1", "Lab-S", "Z2", "Z4," and "Estill" from USDA stock maintenance). These strains were originally all from samples from the US southwest, to minimize potential genetic or cytoplasmic incompatibilities between strains. All possible pairwise crosses indicated no reduction in fitness between any given pair . This diverse population was maintained and well mixed for approximately ten generations at a population size of 10,000 individuals. This population was divided into four subsets spaced 1 week apart (hereafter "temporal blocks"), with gene flow between the different temporal blocks. From this large and diverse population, we created 60 independent replicate populations experiencing bottlenecks in population size of two different intensities over two generations. During the first generation of the bottleneck event, all 60 populations were reduced to a size of two individuals.
During the second generation of the bottleneck event, 24 populations were allowed to triple in size to six individuals to initiate the next generation (hereafter called "intermediate bottleneck") and 36 populations were not allowed to grow, and experienced a second bottleneck of two individuals, to initiate the next generation (hereafter called "strong bottleneck"). Contrary to our expectations, the success rate of bottleneck population creation was the same regardless of bottleneck intensity, thus we successfully created more replicates for strong bottleneck populations than for intermediate bottleneck populations. Then, these 60 populations were maintained for six additional generations in the laboratory to increase population size prior to starting the experiment (Figure 1). In contrast to the intermediate and strong bottlenecked populations, the diverse population was maintained as a single large population during those six generations to maintain high genetic diversity and was used to create independent replicate populations at the beginning of the experiment (hereafter called "diverse populations"). During the six generations prior to the initiation of the experiment, the population size of each bottlenecked population was recorded so that the individual population demographic history within each bottleneck category was known.

| Quantifying the demographic history
By tracking the population size of each population each generation, we know their full demographic history since the bottleneck. We calculate the expected heterozygosity remaining in each population relative to the source population as an index that takes into account the variation in population size through time. Following Allendorf (1986), the expected heterozygosity at generation t (He t ), is a function of expected heterozygosity and population size in generation t−1: where He t−1 is the expected heterozygosity at the previous generation t−1 and N t−1 is the population size at the previous generation t−1.
To calculate the proportion of expected heterozygosity remaining, expected heterozygosity of the source population at the time of the bottleneck is defined as 1, that is He 0 = 1. The most accurate way to calculate this is using effective population size, N e (i.e., the size of an ideal population that would lose genetic variation at the same rate as the focal population; Allendorf, 1986). We used census population size (N), as that is what we had available. In the case of T. castaneum these values are correlated (e.g., N e /N around 0.9 for N = 100 individuals, Pray et al., 1996;Wade, 1980). Because N is slightly larger than N e , our calculation is biased in a conservative way, estimating smaller losses

| Evolutionary rescue experiment
To test the effect of demographic history on evolutionary rescue, we compared the probability of rescue of diverse populations to that of bottlenecked populations. In total, we compared the evolution of 36 diverse populations, 24 intermediate bottleneck populations, and 36 strong bottleneck populations distributed over three temporal blocks. To start the experiment, 100 individuals from each population were housed in two patches connected via 2-mm holes drilled in the sides. We challenged these experimental populations over six discrete generations with a stressful medium that consisted of 97% corn flour, 2.85% wheat flour, and 0.15% brewer's yeast, to represent a resource-limited environment similar to Hufbauer et al. (2015). Each generation, the adults in each population were censused completely, then the individuals were mixed and put back evenly in the two connected patches with fresh challenging media for the 5-week life cycle as described above. For each population, the proportion of expected heterozygosity was computed for each generation (Equation 1) until the end of the experiment or the extinction of the population. A population was considered extinct when the population size was 0 or 1 at census.
At the end of the experiment, the growth rate of each population was measured as the number of offspring produced per 30 adults (replicated as much as possible depending on population size). This number of adults was chosen to enable measurements of the growth rate of all the extant populations, even those with small populations (so as not to bias our findings by excluding smaller, potentially F I G U R E 1 Experimental design. The bottlenecked populations experienced a bottleneck event seven and eight generations before the beginning of the experiment (generations represented by negative values G −7 and G −6). During these eight preliminary generations (Part 1, left), the populations were maintained on a rich environment (grey circles). During the experiment (Part 2, right), the populations evolved on a stressful medium mainly composed of corn flour (yellow circles). An extinction event during the experiment is represented with an empty circle.

| Statistical analysis
All statistical models were fitted in R version 3.4.2 (R Core Team, 2014). For all models, we estimated the significance of interactions and main effects by comparing a full model (i.e., with model terms of the same order) to a reduced model without the interaction or main effect of interest using a Likelihood Ratio Test (LRT).
By using likelihood ratio tests to eliminate nonsignificant effects, we identified the most parsimonious model (hereafter "best-fit model") for each dataset.
To understand how the probability of extinction in the challenging environment was influenced by demographic history, we used generalized linear models for the probability of extinction and time to extinction. Since none of the diverse populations went extinct during the experiment, the parameter estimate for this treatment in standard logistic regression and survival analysis would be infinite and estimates of uncertainty and hypothesis tests would be invalid, a phenomenon known as complete separation (Albert & Anderson, 1984). To manage this issue, we separately analyzed (1) the probability of extinction at the last generation, (2) the time to extinction of the populations that went extinct during the experiment, (3) survival of populations over time, in each case using approaches compatible with complete separation. We modeled the probability of extinction as a Bernoulli distribution with logit link fitted with penalized likelihood and Firth's correction (1993) using the package logistf (Heinze et al., 2022). Explanatory variables in the model were demographic history (a categorical variable with three levels: "diverse", "intermediate bottleneck," and "strong bottleneck") and temporal block (a categorical variable with three levels), both as fixed effects (Table 1). We modeled time to extinction (number of generations) as a Poisson distribution with log link after confirming that overdispersion was absent. Explanatory variables in the model were the categorical demographic history and the temporal block, again as fixed effects (Table 1). We modeled survival over time with a survival analysis using a log-rank statistic to compare the survival dynamics across groups using the package survival (Therneau, 2018).
The explanatory variable in the model was demographic history (a categorical variable with three levels: "diverse," "intermediate bottleneck," and "strong bottleneck") as a fixed effect ( Figure S1). Due to the problem of complete separation, a Cox-PH model was not fitted, as it would result in degenerate estimates of the probability of survival or extinction or time to extinction. Note: See Table S1 for the parameter estimates of the probability of extinction model and the time to extinction model, and Table S2 for the parameter estimates of the population size model.

TA B L E 1 Likelihood ratio test
Abbreviations: D, deviance computed from the likelihood or the penalized likelihood (see Section 2); df, degrees of freedom; p-value, p-value of the LRT; ΔD, change in deviance between the two models (i.e., statistic of the LRT following a χ 2 distribution); Δdf, change in degrees of freedom between the two models.
A U-shaped curve in population size characterizes evolutionary rescue. To evaluate if and how the shape of rescue curves is influenced by demographic history, we focused on populations extant at the end of the experiment and analyzed size through time using a Generalized Additive Mixed Model (GAMM; Wood, 2017) with a negative binomial distribution to account for overdispersion using the package mgcv (Wood, 2022). To allow for the nonlinear relationship, we modeled population size (starting at generation 2, the peak in abundance before the decline) as a smooth function of generation (a continuous variable from 2 to 6) using a thin plate regression spline (Wood, 2017), implemented with the number of "knots" set to 5 using the mgcv package (Wood, 2022). Explanatory variables included categorical demographic history, the interaction between generation and demographic history (thus allowing the shape of the curve to vary by demographic history), and temporal block (three levels) as fixed effects and population identity as a random effect. This model was compared with a reduced model (Table 1)  in our dataset, this direct transformation was not possible. To consider appropriate transformations, we checked the distribution of residuals from growth rate analyses with either a log transformation of growth rate + 1, a square root transformation, or without transformation. This preliminary analysis showed that not transforming the data provided the best approximation of normally distributed residuals. Thus, growth rate was modeled with a normal distribution and identity link, using the package lme4 (Bates et al., 2015). The development rate was evaluated by comparing the presence of the different life stages (i.e., larva, pupa, or adult) weighted by the number in each patch, using a cumulative linkage model for ordinal data to consider the order of appearance of each stage (larva, then pupa, then adult), using the package ordinal (Christensen, 2015). For the growth rate model and for the development rate model, explanatory variables in the model were the categorical demographic history and temporal blocks as fixed effects, with population identity included as a random effect to account for replication within populations (Table 2). In addition, the observation time (4 weeks or 5 weeks) was also included as a fixed effect in the developmental rate model, as well as the replicate identity as a random effect (  Table S3 for the parameter estimates of the growth rate models and the development time model. Abbreviations: D, deviance computed from the likelihood; df, degrees of freedom; p-value, p-value of the LRT; ΔD, change in deviance between the two models (i.e., statistic of the LRT following a χ 2 distribution); Δdf, change in degrees of freedom between the two models. assuming the fitted model was the data-generating process (5000 samples, percentile method; Davison & Hinkley, 1997). Knowing the confidence interval of the difference allows us to differentiate between the case where the populations are biologically similar (small mean difference with small confidence interval reflecting a narrow range of differences compatible with the data) or if we lack evidence to precisely determine the difference between these two types of bottlenecked populations (small mean difference but large confidence interval reflecting a wide range of differences compatible with the data).

| Extinction
There was a significant effect of demographic history on the probability of extinction (p = 0.009, Table 1  For populations that went extinct, there was not necessarily a gradual decrease in population size that led up to extinction ( Figure S2). Instead, population size mostly declined markedly one generation before extinction (Figure 2, red-barred circles in Figure S2), reflecting the sharp population declines.     Figure 3a). These estimates remain similar when using the full model that best represents the experimental design (population size at the last generation for di-  The best-fit GAMM model does not include an interaction between the demographic history and the time, which means that the shape of the predictions cannot differ by demographic history. The full GAMM model includes an interaction between demographic history and time, which means that the shape of the curve can differ by demographic history. The shaded and colored area around each line is the 95% confidence interval. Points in different treatments are staggered slightly and jittered on the generation axis to reduce overlap (generation is discrete).

| Growth rate and development time
Most of the extant populations for the three demographic histories had on average at the end of the experiment a growth rate higher than 1, which means that more descendants were produced than there were parents ( Figure 5). At the end of the experiment, the growth rate was significantly different by demographic history (demography history effect: p < 0.001, Table 2 p < 0.001, Table 2).
The demographic history of rescued populations had little to no influence on developmental timing. The proportions of larvae, pupae, and adults were similar among the three demographic histories ( Figure S3) and there was insufficient evidence for a difference (p = 0.94, Table 2). These proportions varied significantly only with respect to timestep, that is, between 4 weeks and 5 weeks (p < 0.001, Table 2).

| DISCUSS ION
In this study, we evaluated how the demographic history of populations influences the probability of extinction and adaptation to a challenging environment. We found that populations that experienced a bottleneck event in their demographic history were more likely to go extinct, and, if they survived, had smaller population sizes, and were less well adapted to the challenging environment than populations that had not experienced bottleneck events.

| Extinction associated with large initial drop in population size for bottlenecked populations
While all diverse populations persisted throughout the six generations of the experiment, multiple bottlenecked populations went extinct. These extinctions were probably due to the high inbreeding depression and loss of genetic diversity present in bottlenecked populations (Frankham, 2019). Notably, except for one population, none of the populations that experienced an initial drastic drop in population size associated with a 55% reduction in heterozygosity, were able to avoid extinction. This may be due to an extinction vortex, a process in which small populations remain small due to a reinforcing negative feedback loop that makes them vulnerable to extinction due to demographic stochasticity (Gilpin & Soulé, 1986).
The diverse populations did not decline to as small a size as the bottlenecked populations, and perhaps this enabled them to avoid falling into an extinction vortex. In the future, it would be interesting to know if diverse populations that decline more sharply, due to stochasticity or experimental manipulation, can escape more easily via rapid adaptation than populations that have experienced bottleneck events.

| Dynamics of rescue for extant populations
Extant populations from all treatments exhibited the U-shaped population size trajectory characteristic of evolutionary rescue. Diverse populations showed a higher level of adaptation, quantified by population growth rate, than the bottlenecked populations.
These results confirm that rapid adaptation is compromised F I G U R E 4 Coefficient of variation of population size through time. The mean coefficient of variation for each demographic history is represented by a point, with the 95% confidence interval represented by an error bar. Extinct populations were excluded but results were similar when they are considered. Points in different treatments are staggered slightly on the generation axis to reduce overlap (generation is discrete).   (Bouzat, 2010).

| Similar patterns of evolutionary rescue for intermediate and strong bottleneck populations
The lack of support for differing coefficients of variation at the beginning of the experiment might be because all populations were introduced into a new environment and experience it in a similar way. In addition, contrary to theoretical predictions and experimental results (Vogwill et al., 2016; but see Windels et al., 2021), we did not observe significantly higher coefficients of variation for strong bottlenecks than for intermediate bottlenecks. As described above, this lack of evidence for a difference could be due to a lack of statistical power or a lack of difference due to purging or a lack of phenotypic difference (see the section on similar responses between bottlenecked populations for more details).

| Life-history traits involved in the adaptation process
To better understand the different levels of adaptation between treatments, we examined the proportions of individuals in different life stages at two time points at the end of the experiment, to search for potential differences in development time. We found that these proportions were on average very similar among demographic histories. Life-history traits other than development time may thus be responsible for adaptation to the challenging environment. For instance, Agashe et al. (2011) demonstrated an increase in fecundity and egg survival, as well as a higher rate of egg cannibalism, associated with adaptation to a new corn environment for populations of T. castaneum.

| Contributions and ideas for future genomic research projects
While genomic analyses were not within the scope of this study, this project opens up many fruitful avenues of research, as mentioned above, that could be addressed with genomics. We propose three concrete ways in which sequencing these types of experimental or natural populations would address the different perspectives mentioned above. First, sequencing these types of populations would allow us to measure the components of genetic load (i.e., "the proportion by which the population fitness is decreased in comparison with an optimum genotype") as described in Bertorelle et al. (2022) to determine whether genetic load can be used as a good indicator for predicting extinction probabilities. Secondly, sequencing populations whose population dynamics are known would allow us to directly measure and calculate heterozygosity to know if this measure allows us to correctly describe the demographic history of populations (Frankham, 2010;Frankham et al., 2002). Third, by looking at populations adapted to different stressors, it would be interesting to know if the pathways used are more diverse involving a large range of polygenic adaptation, when the populations have or have not undergone complex demographic histories (Gamblin et al., 2023;Garoff et al., 2020).

| Conclusion and perspectives
Our findings have implications for conservation in today's changing world. Managers often have been concerned with maintaining adaptation to local environments for populations of conservation concern (Dalrymple et al., 2021). In the Anthropocene, where environments are changing extremely rapidly, a focus on supporting sufficient genetic variation for ongoing and future adaptation to shifting conditions is important. The reduced ability of bottlenecked populations in this study to adapt to altered conditions suggests that more active management strategies may be needed for populations that have experienced bottlenecks. Therefore, our work supports the substantial literature on genetic rescue, the process of alleviating inbreeding depression and restoring genetic diversity of inbred populations through gene flow (Ralls et al., 2020;Whiteley et al., 2015).

ACK N OWLED G M ENTS
We thank the excellent team of undergraduate students who participated with data collection: D.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R scripts for our analyses are available at: https:// github.com/olazlaure/AnalysisEvolRescueExp and in the Dryad Repository: https://doi.org/10.5061/dryad.xsj3tx9mx.